Line data Source code
1 : module radiation
2 :
3 : !---------------------------------------------------------------------------------
4 : !
5 : ! CAM interface to RRTMG radiation parameterization
6 : !
7 : !---------------------------------------------------------------------------------
8 :
9 : use shr_kind_mod, only: r8=>shr_kind_r8
10 : use spmd_utils, only: masterproc
11 : use ppgrid, only: pcols, pver, pverp, begchunk, endchunk
12 : use physics_types, only: physics_state, physics_ptend
13 : use physics_buffer, only: physics_buffer_desc, pbuf_get_field, pbuf_old_tim_idx
14 : use camsrfexch, only: cam_out_t, cam_in_t
15 : use physconst, only: cappa, cpair
16 :
17 : use time_manager, only: get_nstep, is_first_restart_step, &
18 : get_curr_calday, get_step_size
19 :
20 : use rad_constituents, only: N_DIAG, rad_cnst_get_call_list, &
21 : rad_cnst_get_gas, rad_cnst_out, oldcldoptics, &
22 : liqcldoptics, icecldoptics
23 :
24 : use radconstants, only: nswbands, nlwbands, rrtmg_sw_cloudsim_band, rrtmg_lw_cloudsim_band, &
25 : idx_sw_diag
26 :
27 : use cospsimulator_intr, only: docosp, cospsimulator_intr_init, &
28 : cospsimulator_intr_run, cosp_nradsteps
29 :
30 : use scamMod, only: scm_crm_mode, single_column, have_cld, cldobs
31 :
32 : use cam_history, only: addfld, add_default, horiz_only, outfld, hist_fld_active
33 : use cam_history_support, only: fillvalue
34 :
35 : use pio, only: file_desc_t, var_desc_t, &
36 : pio_int, pio_double, pio_noerr, &
37 : pio_seterrorhandling, pio_bcast_error, &
38 : pio_inq_varid, pio_def_var, &
39 : pio_put_var, pio_get_var, pio_put_att
40 :
41 : use cam_abortutils, only: endrun
42 : use error_messages, only: handle_err
43 : use perf_mod, only: t_startf, t_stopf
44 : use cam_logfile, only: iulog
45 :
46 : implicit none
47 : private
48 : save
49 :
50 : public :: &
51 : radiation_readnl, &! read namelist variables
52 : radiation_register, &! registers radiation physics buffer fields
53 : radiation_nextsw_cday, &! calendar day of next radiation calculation
54 : radiation_do, &! query which radiation calcs are done this timestep
55 : radiation_init, &! initialization
56 : radiation_define_restart, &! define variables for restart
57 : radiation_write_restart, &! write variables to restart
58 : radiation_read_restart, &! read variables from restart
59 : radiation_tend, &! compute heating rates and fluxes
60 : rad_out_t ! type for diagnostic outputs
61 :
62 : integer,public, allocatable :: cosp_cnt(:) ! counter for cosp
63 : integer,public :: cosp_cnt_init = 0 !initial value for cosp counter
64 : real(r8), public, protected :: nextsw_cday ! future radiation calday for surface models
65 :
66 : type rad_out_t
67 : real(r8) :: solin(pcols) ! Solar incident flux
68 :
69 : real(r8) :: qrsc(pcols,pver)
70 :
71 : real(r8) :: fsntc(pcols) ! Clear sky total column abs solar flux
72 : real(r8) :: fsntoa(pcols) ! Net solar flux at TOA
73 : real(r8) :: fsntoac(pcols) ! Clear sky net solar flux at TOA
74 : real(r8) :: fsutoa(pcols) ! upwelling solar flux at TOA
75 :
76 : real(r8) :: fsnirt(pcols) ! Near-IR flux absorbed at toa
77 : real(r8) :: fsnrtc(pcols) ! Clear sky near-IR flux absorbed at toa
78 : real(r8) :: fsnirtsq(pcols) ! Near-IR flux absorbed at toa >= 0.7 microns
79 :
80 : real(r8) :: fsn200(pcols) ! fns interpolated to 200 mb
81 : real(r8) :: fsn200c(pcols) ! fcns interpolated to 200 mb
82 : real(r8) :: fsnr(pcols) ! fns interpolated to tropopause
83 :
84 : real(r8) :: fsnsc(pcols) ! Clear sky surface abs solar flux
85 : real(r8) :: fsdsc(pcols) ! Clear sky surface downwelling solar flux
86 :
87 : real(r8) :: qrlc(pcols,pver)
88 :
89 : real(r8) :: flntc(pcols) ! Clear sky lw flux at model top
90 : real(r8) :: flut(pcols) ! Upward flux at top of model
91 : real(r8) :: flutc(pcols) ! Upward Clear Sky flux at top of model
92 : real(r8) :: lwcf(pcols) ! longwave cloud forcing
93 :
94 : real(r8) :: fln200(pcols) ! net longwave flux interpolated to 200 mb
95 : real(r8) :: fln200c(pcols) ! net clearsky longwave flux interpolated to 200 mb
96 : real(r8) :: flnr(pcols) ! net longwave flux interpolated to tropopause
97 :
98 : real(r8) :: flnsc(pcols) ! Clear sky lw flux at srf (up-down)
99 : real(r8) :: fldsc(pcols) ! Clear sky lw flux at srf (down)
100 :
101 : real(r8) :: tot_cld_vistau(pcols,pver) ! gbx water+ice cloud optical depth (only during day, night = fillvalue)
102 : real(r8) :: tot_icld_vistau(pcols,pver) ! in-cld water+ice cloud optical depth (only during day, night = fillvalue)
103 : real(r8) :: liq_icld_vistau(pcols,pver) ! in-cld liq cloud optical depth (only during day, night = fillvalue)
104 : real(r8) :: ice_icld_vistau(pcols,pver) ! in-cld ice cloud optical depth (only during day, night = fillvalue)
105 : real(r8) :: snow_icld_vistau(pcols,pver) ! snow in-cloud visible sw optical depth for output on history files
106 : real(r8) :: grau_icld_vistau(pcols,pver) ! Graupel in-cloud visible sw optical depth for output on history files
107 :
108 : real(r8) :: cld_tau_cloudsim(pcols,pver)
109 : real(r8) :: aer_tau400(pcols,0:pver)
110 : real(r8) :: aer_tau550(pcols,0:pver)
111 : real(r8) :: aer_tau700(pcols,0:pver)
112 :
113 : end type rad_out_t
114 :
115 : ! Namelist variables
116 :
117 : integer :: iradsw = -1 ! freq. of shortwave radiation calc in time steps (positive)
118 : ! or hours (negative).
119 : integer :: iradlw = -1 ! frequency of longwave rad. calc. in time steps (positive)
120 : ! or hours (negative).
121 :
122 : integer :: irad_always = 0 ! Specifies length of time in timesteps (positive)
123 : ! or hours (negative) SW/LW radiation will be
124 : ! run continuously from the start of an
125 : ! initial or restart run
126 : logical :: use_rad_dt_cosz = .false. ! if true, use radiation dt for all cosz calculations
127 : logical :: spectralflux = .false. ! calculate fluxes (up and down) per band.
128 : logical :: graupel_in_rad = .false. ! graupel in radiation code
129 : logical :: use_rad_uniform_angle = .false. ! if true, use the namelist rad_uniform_angle for the coszrs calculation
130 :
131 :
132 : ! Physics buffer indices
133 : integer :: qrs_idx = 0
134 : integer :: qrl_idx = 0
135 : integer :: su_idx = 0
136 : integer :: sd_idx = 0
137 : integer :: lu_idx = 0
138 : integer :: ld_idx = 0
139 : integer :: fsds_idx = 0
140 : integer :: fsns_idx = 0
141 : integer :: fsnt_idx = 0
142 : integer :: flns_idx = 0
143 : integer :: flnt_idx = 0
144 : integer :: cldfsnow_idx = 0
145 : integer :: cld_idx = 0
146 : integer :: cldfgrau_idx = 0
147 :
148 : character(len=4) :: diag(0:N_DIAG) =(/' ','_d1 ','_d2 ','_d3 ','_d4 ','_d5 ','_d6 ','_d7 ','_d8 ','_d9 ','_d10'/)
149 :
150 : ! averaging time interval for zenith angle
151 : real(r8) :: dt_avg = 0._r8
152 :
153 : real(r8) :: rad_uniform_angle = -99._r8
154 :
155 : ! PIO descriptors (for restarts)
156 : type(var_desc_t) :: cospcnt_desc
157 : type(var_desc_t) :: nextsw_cday_desc
158 : !===============================================================================
159 : contains
160 : !===============================================================================
161 :
162 1536 : subroutine radiation_readnl(nlfile)
163 :
164 : ! Read radiation_nl namelist group.
165 :
166 : use namelist_utils, only: find_group_name
167 : use units, only: getunit, freeunit
168 : use spmd_utils, only: mpicom, mstrid=>masterprocid, mpi_integer, mpi_logical, mpi_real8
169 :
170 :
171 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
172 :
173 : ! Local variables
174 : integer :: unitn, ierr
175 : integer :: dtime ! timestep size
176 : character(len=*), parameter :: sub = 'radiation_readnl'
177 :
178 : namelist /radiation_nl/ iradsw, iradlw, irad_always, &
179 : use_rad_dt_cosz, spectralflux,graupel_in_rad, use_rad_uniform_angle, rad_uniform_angle
180 :
181 : !-----------------------------------------------------------------------------
182 :
183 1536 : if (masterproc) then
184 2 : unitn = getunit()
185 2 : open( unitn, file=trim(nlfile), status='old' )
186 2 : call find_group_name(unitn, 'radiation_nl', status=ierr)
187 2 : if (ierr == 0) then
188 2 : read(unitn, radiation_nl, iostat=ierr)
189 2 : if (ierr /= 0) then
190 0 : call endrun(sub // ':: ERROR reading namelist')
191 : end if
192 : end if
193 2 : close(unitn)
194 2 : call freeunit(unitn)
195 : end if
196 :
197 : ! Broadcast namelist variables
198 1536 : call mpi_bcast(iradsw, 1, mpi_integer, mstrid, mpicom, ierr)
199 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: iradsw")
200 1536 : call mpi_bcast(iradlw, 1, mpi_integer, mstrid, mpicom, ierr)
201 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: iradlw")
202 1536 : call mpi_bcast(irad_always, 1, mpi_integer, mstrid, mpicom, ierr)
203 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: irad_always")
204 1536 : call mpi_bcast(use_rad_dt_cosz, 1, mpi_logical, mstrid, mpicom, ierr)
205 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: use_rad_dt_cosz")
206 1536 : call mpi_bcast(spectralflux, 1, mpi_logical, mstrid, mpicom, ierr)
207 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: spectralflux")
208 1536 : call mpi_bcast(graupel_in_rad, 1, mpi_logical, mstrid, mpicom, ierr)
209 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: graupel_in_rad")
210 1536 : call mpi_bcast(use_rad_uniform_angle, 1, mpi_logical, mstrid, mpicom, ierr)
211 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: use_rad_uniform_angle")
212 1536 : call mpi_bcast(rad_uniform_angle, 1, mpi_real8, mstrid, mpicom, ierr)
213 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: rad_uniform_angle")
214 :
215 1536 : if (use_rad_uniform_angle .and. rad_uniform_angle == -99._r8) then
216 0 : call endrun(sub // ' ERROR - use_rad_uniform_angle is set to .true, but rad_uniform_angle is not set ')
217 : end if
218 :
219 : ! Convert iradsw, iradlw and irad_always from hours to timesteps if necessary
220 1536 : dtime = get_step_size()
221 1536 : if (iradsw < 0) iradsw = nint((-iradsw *3600._r8)/dtime)
222 1536 : if (iradlw < 0) iradlw = nint((-iradlw *3600._r8)/dtime)
223 1536 : if (irad_always < 0) irad_always = nint((-irad_always*3600._r8)/dtime)
224 :
225 : !-----------------------------------------------------------------------
226 : ! Print runtime options to log.
227 : !-----------------------------------------------------------------------
228 :
229 1536 : if (masterproc) then
230 2 : write(iulog,*) 'RRTMG radiation scheme parameters:'
231 2 : write(iulog,10) iradsw, iradlw, irad_always, use_rad_dt_cosz, spectralflux, graupel_in_rad
232 : end if
233 :
234 : 10 format(' Frequency (timesteps) of Shortwave Radiation calc: ',i5/, &
235 : ' Frequency (timesteps) of Longwave Radiation calc: ',i5/, &
236 : ' SW/LW calc done every timestep for first N steps. N=',i5/, &
237 : ' Use average zenith angle: ',l5/, &
238 : ' Output spectrally resolved fluxes: ',l5/, &
239 : ' Graupel in Radiation Code: ',l5/)
240 :
241 1536 : end subroutine radiation_readnl
242 :
243 : !================================================================================================
244 :
245 1536 : subroutine radiation_register
246 :
247 : ! Register radiation fields in the physics buffer
248 :
249 : use physics_buffer, only: pbuf_add_field, dtype_r8
250 : use radiation_data, only: rad_data_register
251 :
252 1536 : call pbuf_add_field('QRS' , 'global',dtype_r8,(/pcols,pver/), qrs_idx) ! shortwave radiative heating rate
253 1536 : call pbuf_add_field('QRL' , 'global',dtype_r8,(/pcols,pver/), qrl_idx) ! longwave radiative heating rate
254 :
255 1536 : call pbuf_add_field('FSDS' , 'global',dtype_r8,(/pcols/), fsds_idx) ! Surface solar downward flux
256 1536 : call pbuf_add_field('FSNS' , 'global',dtype_r8,(/pcols/), fsns_idx) ! Surface net shortwave flux
257 1536 : call pbuf_add_field('FSNT' , 'global',dtype_r8,(/pcols/), fsnt_idx) ! Top-of-model net shortwave flux
258 :
259 1536 : call pbuf_add_field('FLNS' , 'global',dtype_r8,(/pcols/), flns_idx) ! Surface net longwave flux
260 1536 : call pbuf_add_field('FLNT' , 'global',dtype_r8,(/pcols/), flnt_idx) ! Top-of-model net longwave flux
261 :
262 : ! If the namelist has been configured for preserving the spectral fluxes, then create
263 : ! physics buffer variables to store the results.
264 1536 : if (spectralflux) then
265 0 : call pbuf_add_field('SU' , 'global',dtype_r8,(/pcols,pverp,nswbands/), su_idx) ! shortwave upward flux (per band)
266 0 : call pbuf_add_field('SD' , 'global',dtype_r8,(/pcols,pverp,nswbands/), sd_idx) ! shortwave downward flux (per band)
267 0 : call pbuf_add_field('LU' , 'global',dtype_r8,(/pcols,pverp,nlwbands/), lu_idx) ! longwave upward flux (per band)
268 0 : call pbuf_add_field('LD' , 'global',dtype_r8,(/pcols,pverp,nlwbands/), ld_idx) ! longwave downward flux (per band)
269 : end if
270 :
271 1536 : call rad_data_register()
272 :
273 1536 : end subroutine radiation_register
274 :
275 : !================================================================================================
276 :
277 226008 : function radiation_do(op, timestep)
278 :
279 : ! Return true if the specified operation is done this timestep.
280 :
281 : character(len=*), intent(in) :: op ! name of operation
282 : integer, intent(in), optional:: timestep
283 : logical :: radiation_do ! return value
284 :
285 : ! Local variables
286 : integer :: nstep ! current timestep number
287 : !-----------------------------------------------------------------------
288 :
289 226008 : if (present(timestep)) then
290 102168 : nstep = timestep
291 : else
292 123840 : nstep = get_nstep()
293 : end if
294 :
295 164088 : select case (op)
296 :
297 : case ('sw') ! do a shortwave heating calc this timestep?
298 : radiation_do = nstep == 0 .or. iradsw == 1 &
299 : .or. (mod(nstep-1,iradsw) == 0 .and. nstep /= 1) &
300 164088 : .or. nstep <= irad_always
301 :
302 : case ('lw') ! do a longwave heating calc this timestep?
303 : radiation_do = nstep == 0 .or. iradlw == 1 &
304 : .or. (mod(nstep-1,iradlw) == 0 .and. nstep /= 1) &
305 61920 : .or. nstep <= irad_always
306 :
307 : case default
308 226008 : call endrun('radiation_do: unknown operation:'//op)
309 :
310 : end select
311 1536 : end function radiation_do
312 :
313 : !================================================================================================
314 :
315 61920 : real(r8) function radiation_nextsw_cday()
316 :
317 : ! Return calendar day of next sw radiation calculation
318 :
319 : ! Local variables
320 : integer :: nstep ! timestep counter
321 : logical :: dosw ! true => do shosrtwave calc
322 : integer :: offset ! offset for calendar day calculation
323 : integer :: dtime ! integer timestep size
324 : real(r8):: calday ! calendar day of
325 : real(r8):: caldayp1 ! calendar day of next time-step
326 : !-----------------------------------------------------------------------
327 :
328 61920 : radiation_nextsw_cday = -1._r8
329 61920 : dosw = .false.
330 61920 : nstep = get_nstep()
331 61920 : dtime = get_step_size()
332 61920 : offset = 0
333 : do while (.not. dosw)
334 102168 : nstep = nstep + 1
335 102168 : offset = offset + dtime
336 102168 : if (radiation_do('sw', nstep)) then
337 61920 : radiation_nextsw_cday = get_curr_calday(offset=offset)
338 : dosw = .true.
339 : end if
340 : end do
341 61920 : if(radiation_nextsw_cday == -1._r8) then
342 0 : call endrun('error in radiation_nextsw_cday')
343 : end if
344 :
345 : ! determine if next radiation time-step not equal to next time-step
346 61920 : if (get_nstep() >= 1) then
347 55728 : caldayp1 = get_curr_calday(offset=int(dtime))
348 55728 : if (caldayp1 /= radiation_nextsw_cday) radiation_nextsw_cday = -1._r8
349 : end if
350 :
351 61920 : end function radiation_nextsw_cday
352 :
353 : !================================================================================================
354 :
355 1536 : subroutine radiation_init(pbuf2d)
356 :
357 : ! Initialize the radiation parameterization, add fields to the history buffer
358 :
359 : use physics_buffer, only: pbuf_get_index, pbuf_set_field
360 : use phys_control, only: phys_getopts
361 : use radsw, only: radsw_init
362 : use radlw, only: radlw_init
363 : use rad_solar_var, only: rad_solar_var_init
364 : use radiation_data, only: rad_data_init
365 : use cloud_rad_props, only: cloud_rad_props_init
366 : use rrtmg_state, only: rrtmg_state_init
367 : use time_manager, only: is_first_step
368 :
369 :
370 : ! arguments
371 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
372 :
373 : ! local variables
374 : integer :: icall
375 : logical :: active_calls(0:N_DIAG)
376 : integer :: nstep ! current timestep number
377 : logical :: history_amwg ! output the variables used by the AMWG diag package
378 : logical :: history_vdiag ! output the variables used by the AMWG variability diag package
379 : logical :: history_budget ! output tendencies and state variables for CAM4
380 : ! temperature, water vapor, cloud ice and cloud
381 : ! liquid budgets.
382 : integer :: history_budget_histfile_num ! output history file number for budget fields
383 : integer :: err
384 :
385 : integer :: dtime
386 : !-----------------------------------------------------------------------
387 :
388 1536 : call rad_solar_var_init()
389 1536 : call rrtmg_state_init()
390 1536 : call rad_data_init(pbuf2d) ! initialize output fields for offline driver
391 1536 : call radsw_init()
392 1536 : call radlw_init()
393 1536 : call cloud_rad_props_init()
394 :
395 1536 : cld_idx = pbuf_get_index('CLD')
396 1536 : cldfsnow_idx = pbuf_get_index('CLDFSNOW',errcode=err)
397 1536 : cldfgrau_idx = pbuf_get_index('CLDFGRAU',errcode=err)
398 :
399 1536 : if (is_first_step()) then
400 768 : call pbuf_set_field(pbuf2d, qrl_idx, 0._r8)
401 : end if
402 :
403 : ! Set the radiation timestep for cosz calculations if requested using the adjusted iradsw value from radiation
404 1536 : if (use_rad_dt_cosz) then
405 0 : dtime = get_step_size()
406 0 : dt_avg = iradsw*dtime
407 : end if
408 :
409 : ! Surface components to get radiation computed today
410 1536 : if (.not. is_first_restart_step()) then
411 768 : nextsw_cday = get_curr_calday()
412 : end if
413 :
414 : call phys_getopts(history_amwg_out = history_amwg, &
415 : history_vdiag_out = history_vdiag, &
416 : history_budget_out = history_budget, &
417 1536 : history_budget_histfile_num_out = history_budget_histfile_num)
418 :
419 : ! "irad_always" is number of time steps to execute radiation continuously from start of
420 : ! initial OR restart run
421 1536 : nstep = get_nstep()
422 1536 : if (irad_always > 0) then
423 0 : nstep = get_nstep()
424 0 : irad_always = irad_always + nstep
425 : end if
426 :
427 1536 : if (docosp) call cospsimulator_intr_init
428 :
429 4608 : allocate(cosp_cnt(begchunk:endchunk))
430 1536 : if (is_first_restart_step()) then
431 3864 : cosp_cnt(begchunk:endchunk) = cosp_cnt_init
432 : else
433 3864 : cosp_cnt(begchunk:endchunk) = 0
434 : end if
435 :
436 1536 : call addfld('O3colAbove', horiz_only, 'A', 'DU', 'Column O3 above model top', sampling_seq='rad_lwsw')
437 :
438 : call addfld('TOT_CLD_VISTAU', (/ 'lev' /), 'A', '1', 'Total gbx cloud extinction visible sw optical depth', &
439 3072 : sampling_seq='rad_lwsw', flag_xyfill=.true.)
440 : call addfld('TOT_ICLD_VISTAU', (/ 'lev' /), 'A', '1', 'Total in-cloud extinction visible sw optical depth', &
441 3072 : sampling_seq='rad_lwsw', flag_xyfill=.true.)
442 : call addfld('LIQ_ICLD_VISTAU', (/ 'lev' /), 'A', '1', 'Liquid in-cloud extinction visible sw optical depth', &
443 3072 : sampling_seq='rad_lwsw', flag_xyfill=.true.)
444 : call addfld('ICE_ICLD_VISTAU', (/ 'lev' /), 'A', '1', 'Ice in-cloud extinction visible sw optical depth', &
445 3072 : sampling_seq='rad_lwsw', flag_xyfill=.true.)
446 :
447 1536 : if (cldfsnow_idx > 0) then
448 : call addfld('SNOW_ICLD_VISTAU', (/ 'lev' /), 'A', '1', 'Snow in-cloud extinction visible sw optical depth', &
449 3072 : sampling_seq='rad_lwsw', flag_xyfill=.true.)
450 : endif
451 1536 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
452 : call addfld('GRAU_ICLD_VISTAU', (/ 'lev' /), 'A', '1', 'Graupel in-cloud extinction visible sw optical depth', &
453 0 : sampling_seq='rad_lwsw', flag_xyfill=.true.)
454 : endif
455 :
456 : ! get list of active radiation calls
457 1536 : call rad_cnst_get_call_list(active_calls)
458 :
459 : ! Add shortwave radiation fields to history master field list.
460 :
461 18432 : do icall = 0, N_DIAG
462 :
463 18432 : if (active_calls(icall)) then
464 :
465 1536 : call addfld('SOLIN'//diag(icall), horiz_only, 'A', 'W/m2', 'Solar insolation', sampling_seq='rad_lwsw')
466 :
467 3072 : call addfld('QRS'//diag(icall), (/ 'lev' /), 'A', 'K/s', 'Solar heating rate', sampling_seq='rad_lwsw')
468 : call addfld('QRSC'//diag(icall), (/ 'lev' /), 'A', 'K/s', 'Clearsky solar heating rate', &
469 3072 : sampling_seq='rad_lwsw')
470 : call addfld('FSNT'//diag(icall), horiz_only, 'A', 'W/m2', 'Net solar flux at top of model', &
471 1536 : sampling_seq='rad_lwsw')
472 : call addfld('FSNTC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky net solar flux at top of model', &
473 1536 : sampling_seq='rad_lwsw')
474 : call addfld('FSNTOA'//diag(icall), horiz_only, 'A', 'W/m2', 'Net solar flux at top of atmosphere', &
475 1536 : sampling_seq='rad_lwsw')
476 : call addfld('FSNTOAC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky net solar flux at top of atmosphere', &
477 1536 : sampling_seq='rad_lwsw')
478 : call addfld('SWCF'//diag(icall), horiz_only, 'A', 'W/m2', 'Shortwave cloud forcing', &
479 1536 : sampling_seq='rad_lwsw')
480 : call addfld('FSUTOA'//diag(icall), horiz_only, 'A', 'W/m2', 'Upwelling solar flux at top of atmosphere', &
481 1536 : sampling_seq='rad_lwsw')
482 : call addfld('FSNIRTOA'//diag(icall), horiz_only, 'A', 'W/m2', &
483 1536 : 'Net near-infrared flux (Nimbus-7 WFOV) at top of atmosphere', sampling_seq='rad_lwsw')
484 : call addfld('FSNRTOAC'//diag(icall), horiz_only, 'A', 'W/m2', &
485 1536 : 'Clearsky net near-infrared flux (Nimbus-7 WFOV) at top of atmosphere', sampling_seq='rad_lwsw')
486 : call addfld('FSNRTOAS'//diag(icall), horiz_only, 'A', 'W/m2', &
487 1536 : 'Net near-infrared flux (>= 0.7 microns) at top of atmosphere', sampling_seq='rad_lwsw')
488 :
489 : call addfld('FSN200'//diag(icall), horiz_only, 'A', 'W/m2', 'Net shortwave flux at 200 mb', &
490 1536 : sampling_seq='rad_lwsw')
491 : call addfld('FSN200C'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky net shortwave flux at 200 mb', &
492 1536 : sampling_seq='rad_lwsw')
493 :
494 : call addfld('FSNR'//diag(icall), horiz_only, 'A', 'W/m2', 'Net solar flux at tropopause', &
495 1536 : sampling_seq='rad_lwsw')
496 :
497 : call addfld('SOLL'//diag(icall), horiz_only, 'A', 'W/m2', 'Solar downward near infrared direct to surface', &
498 1536 : sampling_seq='rad_lwsw')
499 : call addfld('SOLS'//diag(icall), horiz_only, 'A', 'W/m2', 'Solar downward visible direct to surface', &
500 1536 : sampling_seq='rad_lwsw')
501 : call addfld('SOLLD'//diag(icall), horiz_only, 'A', 'W/m2', 'Solar downward near infrared diffuse to surface', &
502 1536 : sampling_seq='rad_lwsw')
503 : call addfld('SOLSD'//diag(icall), horiz_only, 'A', 'W/m2', 'Solar downward visible diffuse to surface', &
504 1536 : sampling_seq='rad_lwsw')
505 : call addfld('FSNS'//diag(icall), horiz_only, 'A', 'W/m2', 'Net solar flux at surface', &
506 1536 : sampling_seq='rad_lwsw')
507 : call addfld('FSNSC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky net solar flux at surface', &
508 1536 : sampling_seq='rad_lwsw')
509 :
510 : call addfld('FSDS'//diag(icall), horiz_only, 'A', 'W/m2', 'Downwelling solar flux at surface', &
511 1536 : sampling_seq='rad_lwsw')
512 : call addfld('FSDSC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky downwelling solar flux at surface', &
513 1536 : sampling_seq='rad_lwsw')
514 :
515 3072 : call addfld('FUS'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', 'Shortwave upward flux')
516 3072 : call addfld('FDS'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', 'Shortwave downward flux')
517 3072 : call addfld('FUSC'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', 'Shortwave clear-sky upward flux')
518 3072 : call addfld('FDSC'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', 'Shortwave clear-sky downward flux')
519 :
520 1536 : if (history_amwg) then
521 1536 : call add_default('SOLIN'//diag(icall), 1, ' ')
522 1536 : call add_default('QRS'//diag(icall), 1, ' ')
523 1536 : call add_default('FSNT'//diag(icall), 1, ' ')
524 1536 : call add_default('FSNTC'//diag(icall), 1, ' ')
525 1536 : call add_default('FSNTOA'//diag(icall), 1, ' ')
526 1536 : call add_default('FSNTOAC'//diag(icall), 1, ' ')
527 1536 : call add_default('SWCF'//diag(icall), 1, ' ')
528 1536 : call add_default('FSNS'//diag(icall), 1, ' ')
529 1536 : call add_default('FSNSC'//diag(icall), 1, ' ')
530 1536 : call add_default('FSUTOA'//diag(icall), 1, ' ')
531 1536 : call add_default('FSDSC'//diag(icall), 1, ' ')
532 1536 : call add_default('FSDS'//diag(icall), 1, ' ')
533 : endif
534 :
535 : end if
536 : end do
537 :
538 1536 : if (scm_crm_mode) then
539 0 : call add_default('FUS ', 1, ' ')
540 0 : call add_default('FUSC ', 1, ' ')
541 0 : call add_default('FDS ', 1, ' ')
542 0 : call add_default('FDSC ', 1, ' ')
543 : endif
544 :
545 : ! Add longwave radiation fields to history master field list.
546 :
547 18432 : do icall = 0, N_DIAG
548 :
549 18432 : if (active_calls(icall)) then
550 :
551 3072 : call addfld('QRL'//diag(icall), (/ 'lev' /), 'A', 'K/s', 'Longwave heating rate', sampling_seq='rad_lwsw')
552 : call addfld('QRLC'//diag(icall), (/ 'lev' /), 'A', 'K/s', 'Clearsky longwave heating rate', &
553 3072 : sampling_seq='rad_lwsw')
554 : call addfld('FLNT'//diag(icall), horiz_only, 'A', 'W/m2', 'Net longwave flux at top of model', &
555 1536 : sampling_seq='rad_lwsw')
556 : call addfld('FLNTC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky net longwave flux at top of model', &
557 1536 : sampling_seq='rad_lwsw')
558 : call addfld('FLNTCLR'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky ONLY points net longwave flux at top of model',&
559 1536 : sampling_seq='rad_lwsw')
560 : call addfld('FREQCLR'//diag(icall), horiz_only, 'A', 'Frac', 'Frequency of Occurrence of Clearsky', &
561 1536 : sampling_seq='rad_lwsw')
562 : call addfld('FLUT'//diag(icall), horiz_only, 'A', 'W/m2', 'Upwelling longwave flux at top of model', &
563 1536 : sampling_seq='rad_lwsw')
564 : call addfld('FLUTC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky upwelling longwave flux at top of model', &
565 1536 : sampling_seq='rad_lwsw')
566 1536 : call addfld('LWCF'//diag(icall), horiz_only, 'A', 'W/m2', 'Longwave cloud forcing', sampling_seq='rad_lwsw')
567 :
568 : call addfld('FLN200'//diag(icall), horiz_only, 'A', 'W/m2', 'Net longwave flux at 200 mb', &
569 1536 : sampling_seq='rad_lwsw')
570 : call addfld('FLN200C'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky net longwave flux at 200 mb', &
571 1536 : sampling_seq='rad_lwsw')
572 : call addfld('FLNR'//diag(icall), horiz_only, 'A', 'W/m2', 'Net longwave flux at tropopause', &
573 1536 : sampling_seq='rad_lwsw')
574 :
575 : call addfld('FLNS'//diag(icall), horiz_only, 'A', 'W/m2', 'Net longwave flux at surface', &
576 1536 : sampling_seq='rad_lwsw')
577 : call addfld('FLNSC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky net longwave flux at surface', &
578 1536 : sampling_seq='rad_lwsw')
579 : call addfld('FLDS'//diag(icall), horiz_only, 'A', 'W/m2', 'Downwelling longwave flux at surface', &
580 1536 : sampling_seq='rad_lwsw')
581 : call addfld('FLDSC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky Downwelling longwave flux at surface', &
582 1536 : sampling_seq='rad_lwsw')
583 3072 : call addfld('FUL'//diag(icall), (/ 'ilev' /),'I', 'W/m2', 'Longwave upward flux')
584 3072 : call addfld('FDL'//diag(icall), (/ 'ilev' /),'I', 'W/m2', 'Longwave downward flux')
585 3072 : call addfld('FULC'//diag(icall), (/ 'ilev' /),'I', 'W/m2', 'Longwave clear-sky upward flux')
586 3072 : call addfld('FDLC'//diag(icall), (/ 'ilev' /),'I', 'W/m2', 'Longwave clear-sky downward flux')
587 :
588 1536 : if (history_amwg) then
589 1536 : call add_default('QRL'//diag(icall), 1, ' ')
590 1536 : call add_default('FLNT'//diag(icall), 1, ' ')
591 1536 : call add_default('FLNTC'//diag(icall), 1, ' ')
592 1536 : call add_default('FLNTCLR'//diag(icall), 1, ' ')
593 1536 : call add_default('FREQCLR'//diag(icall), 1, ' ')
594 1536 : call add_default('FLUT'//diag(icall), 1, ' ')
595 1536 : call add_default('FLUTC'//diag(icall), 1, ' ')
596 1536 : call add_default('LWCF'//diag(icall), 1, ' ')
597 1536 : call add_default('FLNS'//diag(icall), 1, ' ')
598 1536 : call add_default('FLNSC'//diag(icall), 1, ' ')
599 1536 : call add_default('FLDS'//diag(icall), 1, ' ')
600 : endif
601 :
602 : end if
603 : end do
604 :
605 3072 : call addfld('EMIS', (/ 'lev' /), 'A', '1', 'Cloud longwave emissivity')
606 :
607 1536 : if (scm_crm_mode) then
608 0 : call add_default ('FUL ', 1, ' ')
609 0 : call add_default ('FULC ', 1, ' ')
610 0 : call add_default ('FDL ', 1, ' ')
611 0 : call add_default ('FDLC ', 1, ' ')
612 : endif
613 :
614 : ! Heating rate needed for d(theta)/dt computation
615 3072 : call addfld ('HR',(/ 'lev' /), 'A','K/s','Heating rate needed for d(theta)/dt computation')
616 :
617 1536 : if ( history_budget .and. history_budget_histfile_num > 1 ) then
618 0 : call add_default ('QRL ', history_budget_histfile_num, ' ')
619 0 : call add_default ('QRS ', history_budget_histfile_num, ' ')
620 : end if
621 :
622 1536 : if (history_vdiag) then
623 0 : call add_default('FLUT', 2, ' ')
624 0 : call add_default('FLUT', 3, ' ')
625 : end if
626 :
627 3072 : end subroutine radiation_init
628 :
629 : !===============================================================================
630 :
631 1536 : subroutine radiation_define_restart(file)
632 :
633 : ! define variables to be written to restart file
634 :
635 : ! arguments
636 : type(file_desc_t), intent(inout) :: file
637 :
638 : ! local variables
639 : integer :: ierr
640 : !----------------------------------------------------------------------------
641 :
642 1536 : call pio_seterrorhandling(File, PIO_BCAST_ERROR)
643 :
644 1536 : ierr = pio_def_var(File, 'nextsw_cday', pio_double, nextsw_cday_desc)
645 1536 : ierr = pio_put_att(File, nextsw_cday_desc, 'long_name', 'future radiation calday for surface models')
646 1536 : if (docosp) then
647 0 : ierr = pio_def_var(File, 'cosp_cnt_init', pio_int, cospcnt_desc)
648 : end if
649 :
650 1536 : end subroutine radiation_define_restart
651 :
652 : !===============================================================================
653 :
654 1536 : subroutine radiation_write_restart(file)
655 :
656 : ! write variables to restart file
657 :
658 : ! arguments
659 : type(file_desc_t), intent(inout) :: file
660 :
661 : ! local variables
662 : integer :: ierr
663 : !----------------------------------------------------------------------------
664 :
665 3072 : ierr = pio_put_var(File, nextsw_cday_desc, (/ nextsw_cday /))
666 1536 : if (docosp) then
667 0 : ierr = pio_put_var(File, cospcnt_desc, (/cosp_cnt(begchunk)/))
668 : end if
669 :
670 1536 : end subroutine radiation_write_restart
671 :
672 : !===============================================================================
673 :
674 768 : subroutine radiation_read_restart(file)
675 :
676 : ! read variables from restart file
677 :
678 : ! arguments
679 : type(file_desc_t), intent(inout) :: file
680 :
681 : ! local variables
682 :
683 : integer :: err_handling
684 : integer :: ierr
685 : real(r8) :: temp_var
686 :
687 : type(var_desc_t) :: vardesc
688 : !----------------------------------------------------------------------------
689 :
690 768 : if (docosp) then
691 0 : call pio_seterrorhandling(File, PIO_BCAST_ERROR, err_handling)
692 0 : ierr = pio_inq_varid(File, 'cosp_cnt_init', vardesc)
693 0 : call pio_seterrorhandling(File, err_handling)
694 0 : if (ierr /= PIO_NOERR) then
695 0 : cosp_cnt_init = 0
696 : else
697 0 : ierr = pio_get_var(File, vardesc, cosp_cnt_init)
698 : end if
699 : end if
700 :
701 768 : ierr = pio_inq_varid(File, 'nextsw_cday', vardesc)
702 768 : ierr = pio_get_var(File, vardesc, temp_var)
703 768 : nextsw_cday = temp_var
704 :
705 768 : end subroutine radiation_read_restart
706 :
707 : !===============================================================================
708 :
709 2600640 : subroutine radiation_tend( &
710 : state, ptend, pbuf, cam_out, cam_in, net_flx, rd_out)
711 :
712 : !-----------------------------------------------------------------------
713 : !
714 : ! Driver for radiation computation.
715 : !
716 : ! Revision history:
717 : ! 2007-11-05 M. Iacono Install rrtmg_lw and sw as radiation model.
718 : ! 2007-12-27 M. Iacono Modify to use CAM cloud optical properties with rrtmg.
719 : !-----------------------------------------------------------------------
720 :
721 : use phys_grid, only: get_rlat_all_p, get_rlon_all_p
722 : use cam_control_mod, only: eccen, mvelpp, lambm0, obliqr
723 : use shr_orb_mod, only: shr_orb_decl, shr_orb_cosz
724 :
725 : use aer_rad_props, only: aer_rad_props_sw, aer_rad_props_lw
726 :
727 : use cloud_rad_props, only: get_ice_optics_sw, get_liquid_optics_sw, liquid_cloud_get_rad_props_lw, &
728 : ice_cloud_get_rad_props_lw, cloud_rad_props_get_lw, &
729 : grau_cloud_get_rad_props_lw, get_grau_optics_sw, &
730 : snow_cloud_get_rad_props_lw, get_snow_optics_sw
731 : use slingo_liq_optics, only: slingo_liq_get_rad_props_lw, slingo_liq_optics_sw
732 : use ebert_curry_ice_optics, only: ec_ice_optics_sw, ec_ice_get_rad_props_lw
733 :
734 : use rad_solar_var, only: get_variability
735 : use radsw, only: rad_rrtmg_sw
736 : use radlw, only: rad_rrtmg_lw
737 : use radheat, only: radheat_tend
738 :
739 : use radiation_data, only: rad_data_write
740 : use rrtmg_state, only: rrtmg_state_create, rrtmg_state_update, rrtmg_state_destroy, rrtmg_state_t, &
741 : num_rrtmg_levs
742 :
743 : use interpolate_data, only: vertinterp
744 : use tropopause, only: tropopause_find, TROP_ALG_HYBSTOB, TROP_ALG_CLIMATE
745 :
746 : use cospsimulator_intr, only: docosp, cospsimulator_intr_run, cosp_nradsteps
747 :
748 : ! Arguments
749 : type(physics_state), intent(in), target :: state
750 : type(physics_ptend), intent(out) :: ptend
751 :
752 : type(physics_buffer_desc), pointer :: pbuf(:)
753 : type(cam_out_t), intent(inout) :: cam_out
754 : type(cam_in_t), intent(in) :: cam_in
755 : real(r8), intent(out) :: net_flx(pcols)
756 :
757 : type(rad_out_t), target, optional, intent(out) :: rd_out
758 :
759 :
760 : ! Local variables
761 : type(rad_out_t), pointer :: rd ! allow rd_out to be optional by allocating a local object
762 : ! if the argument is not present
763 : logical :: write_output
764 :
765 : integer :: i, k
766 : integer :: lchnk, ncol
767 : logical :: dosw, dolw
768 : real(r8) :: calday ! current calendar day
769 : real(r8) :: delta ! Solar declination angle in radians
770 : real(r8) :: eccf ! Earth orbit eccentricity factor
771 : real(r8) :: clat(pcols) ! current latitudes(radians)
772 : real(r8) :: clon(pcols) ! current longitudes(radians)
773 : real(r8) :: coszrs(pcols) ! Cosine solar zenith angle
774 :
775 : ! Gathered indices of day and night columns
776 : ! chunk_column_index = IdxDay(daylight_column_index)
777 : integer :: Nday ! Number of daylight columns
778 : integer :: Nnite ! Number of night columns
779 : integer :: IdxDay(pcols) ! Indices of daylight columns
780 : integer :: IdxNite(pcols) ! Indices of night columns
781 :
782 : integer :: itim_old
783 :
784 61920 : real(r8), pointer :: cld(:,:) ! cloud fraction
785 61920 : real(r8), pointer :: cldfsnow(:,:) ! cloud fraction of just "snow clouds- whatever they are"
786 61920 : real(r8), pointer :: cldfgrau(:,:) ! cloud fraction of just "snow clouds- whatever they are"
787 61920 : real(r8), pointer :: qrs(:,:) ! shortwave radiative heating rate
788 61920 : real(r8), pointer :: qrl(:,:) ! longwave radiative heating rate
789 61920 : real(r8), pointer :: fsds(:) ! Surface solar down flux
790 61920 : real(r8), pointer :: fsns(:) ! Surface solar absorbed flux
791 61920 : real(r8), pointer :: fsnt(:) ! Net column abs solar flux at model top
792 61920 : real(r8), pointer :: flns(:) ! Srf longwave cooling (up-down) flux
793 61920 : real(r8), pointer :: flnt(:) ! Net outgoing lw flux at model top
794 :
795 : real(r8), pointer, dimension(:,:,:) :: su => NULL() ! shortwave spectral flux up
796 : real(r8), pointer, dimension(:,:,:) :: sd => NULL() ! shortwave spectral flux down
797 : real(r8), pointer, dimension(:,:,:) :: lu => NULL() ! longwave spectral flux up
798 : real(r8), pointer, dimension(:,:,:) :: ld => NULL() ! longwave spectral flux down
799 :
800 : ! tropopause diagnostic
801 : integer :: troplev(pcols)
802 : real(r8) :: p_trop(pcols)
803 :
804 61920 : type(rrtmg_state_t) :: r_state ! contains the atm concentrations in layers needed for RRTMG
805 :
806 : ! cloud radiative parameters are "in cloud" not "in cell"
807 : real(r8) :: ice_tau (nswbands,pcols,pver) ! ice extinction optical depth
808 : real(r8) :: ice_tau_w (nswbands,pcols,pver) ! ice single scattering albedo * tau
809 : real(r8) :: ice_tau_w_g(nswbands,pcols,pver) ! ice asymmetry parameter * tau * w
810 : real(r8) :: ice_tau_w_f(nswbands,pcols,pver) ! ice forward scattered fraction * tau * w
811 : real(r8) :: ice_lw_abs (nlwbands,pcols,pver) ! ice absorption optics depth (LW)
812 :
813 : ! cloud radiative parameters are "in cloud" not "in cell"
814 : real(r8) :: liq_tau (nswbands,pcols,pver) ! liquid extinction optical depth
815 : real(r8) :: liq_tau_w (nswbands,pcols,pver) ! liquid single scattering albedo * tau
816 : real(r8) :: liq_tau_w_g(nswbands,pcols,pver) ! liquid asymmetry parameter * tau * w
817 : real(r8) :: liq_tau_w_f(nswbands,pcols,pver) ! liquid forward scattered fraction * tau * w
818 : real(r8) :: liq_lw_abs (nlwbands,pcols,pver) ! liquid absorption optics depth (LW)
819 :
820 : ! cloud radiative parameters are "in cloud" not "in cell"
821 : real(r8) :: cld_tau (nswbands,pcols,pver) ! cloud extinction optical depth
822 : real(r8) :: cld_tau_w (nswbands,pcols,pver) ! cloud single scattering albedo * tau
823 : real(r8) :: cld_tau_w_g(nswbands,pcols,pver) ! cloud asymmetry parameter * w * tau
824 : real(r8) :: cld_tau_w_f(nswbands,pcols,pver) ! cloud forward scattered fraction * w * tau
825 : real(r8) :: cld_lw_abs (nlwbands,pcols,pver) ! cloud absorption optics depth (LW)
826 :
827 : ! cloud radiative parameters are "in cloud" not "in cell"
828 : real(r8) :: snow_tau (nswbands,pcols,pver) ! snow extinction optical depth
829 : real(r8) :: snow_tau_w (nswbands,pcols,pver) ! snow single scattering albedo * tau
830 : real(r8) :: snow_tau_w_g(nswbands,pcols,pver) ! snow asymmetry parameter * tau * w
831 : real(r8) :: snow_tau_w_f(nswbands,pcols,pver) ! snow forward scattered fraction * tau * w
832 : real(r8) :: snow_lw_abs (nlwbands,pcols,pver)! snow absorption optics depth (LW)
833 :
834 : ! Add graupel as another snow species.
835 : ! cloud radiative parameters are "in cloud" not "in cell"
836 : real(r8) :: grau_tau (nswbands,pcols,pver) ! graupel extinction optical depth
837 : real(r8) :: grau_tau_w (nswbands,pcols,pver) ! graupel single scattering albedo * tau
838 : real(r8) :: grau_tau_w_g(nswbands,pcols,pver) ! graupel asymmetry parameter * tau * w
839 : real(r8) :: grau_tau_w_f(nswbands,pcols,pver) ! graupel forward scattered fraction * tau * w
840 : real(r8) :: grau_lw_abs (nlwbands,pcols,pver)! graupel absorption optics depth (LW)
841 :
842 : ! combined cloud radiative parameters are "in cloud" not "in cell"
843 : real(r8) :: cldfprime(pcols,pver) ! combined cloud fraction (snow plus regular)
844 : real(r8) :: c_cld_tau (nswbands,pcols,pver) ! combined cloud extinction optical depth
845 : real(r8) :: c_cld_tau_w (nswbands,pcols,pver) ! combined cloud single scattering albedo * tau
846 : real(r8) :: c_cld_tau_w_g(nswbands,pcols,pver) ! combined cloud asymmetry parameter * w * tau
847 : real(r8) :: c_cld_tau_w_f(nswbands,pcols,pver) ! combined cloud forward scattered fraction * w * tau
848 : real(r8) :: c_cld_lw_abs (nlwbands,pcols,pver) ! combined cloud absorption optics depth (LW)
849 :
850 : real(r8) :: sfac(1:nswbands) ! time varying scaling factors due to Solar Spectral Irrad at 1 A.U. per band
851 :
852 : integer :: icall ! index through climate/diagnostic radiation calls
853 : logical :: active_calls(0:N_DIAG)
854 :
855 : ! Aerosol radiative properties
856 : real(r8) :: aer_tau (pcols,0:pver,nswbands) ! aerosol extinction optical depth
857 : real(r8) :: aer_tau_w (pcols,0:pver,nswbands) ! aerosol single scattering albedo * tau
858 : real(r8) :: aer_tau_w_g(pcols,0:pver,nswbands) ! aerosol asymmetry parameter * w * tau
859 : real(r8) :: aer_tau_w_f(pcols,0:pver,nswbands) ! aerosol forward scattered fraction * w * tau
860 : real(r8) :: aer_lw_abs (pcols,pver,nlwbands) ! aerosol absorption optics depth (LW)
861 :
862 : real(r8) :: fns(pcols,pverp) ! net shortwave flux
863 : real(r8) :: fcns(pcols,pverp) ! net clear-sky shortwave flux
864 : real(r8) :: fnl(pcols,pverp) ! net longwave flux
865 : real(r8) :: fcnl(pcols,pverp) ! net clear-sky longwave flux
866 :
867 : ! for COSP
868 : real(r8) :: emis(pcols,pver) ! Cloud longwave emissivity
869 : real(r8) :: gb_snow_tau(pcols,pver) ! grid-box mean snow_tau
870 : real(r8) :: gb_snow_lw(pcols,pver) ! grid-box mean LW snow optical depth
871 :
872 : real(r8) :: ftem(pcols,pver) ! Temporary workspace for outfld variables
873 :
874 : real(r8) :: freqclr(pcols) ! Frequency of occurrence of clear sky columns
875 : real(r8) :: flntclr(pcols) ! Clearsky only columns (zero if cloudy)
876 :
877 : character(*), parameter :: name = 'radiation_tend'
878 : !--------------------------------------------------------------------------------------
879 :
880 61920 : lchnk = state%lchnk
881 61920 : ncol = state%ncol
882 :
883 61920 : if (present(rd_out)) then
884 : rd => rd_out
885 : write_output = .false.
886 : else
887 61920 : allocate(rd)
888 : write_output=.true.
889 : end if
890 :
891 61920 : dosw = radiation_do('sw') ! do shortwave heating calc this timestep?
892 61920 : dolw = radiation_do('lw') ! do longwave heating calc this timestep?
893 :
894 : ! Cosine solar zenith angle for current time step
895 61920 : calday = get_curr_calday()
896 61920 : call get_rlat_all_p(lchnk, ncol, clat)
897 61920 : call get_rlon_all_p(lchnk, ncol, clon)
898 :
899 : call shr_orb_decl(calday, eccen, mvelpp, lambm0, obliqr, &
900 61920 : delta, eccf)
901 :
902 61920 : if (use_rad_uniform_angle) then
903 0 : do i = 1, ncol
904 0 : coszrs(i) = shr_orb_cosz(calday, clat(i), clon(i), delta, dt_avg, uniform_angle=rad_uniform_angle)
905 : end do
906 : else
907 1033920 : do i = 1, ncol
908 1033920 : coszrs(i) = shr_orb_cosz(calday, clat(i), clon(i), delta, dt_avg)
909 : end do
910 : end if
911 :
912 : ! Gather night/day column indices.
913 61920 : Nday = 0
914 61920 : Nnite = 0
915 1033920 : do i = 1, ncol
916 1033920 : if ( coszrs(i) > 0.0_r8 ) then
917 486000 : Nday = Nday + 1
918 486000 : IdxDay(Nday) = i
919 : else
920 486000 : Nnite = Nnite + 1
921 486000 : IdxNite(Nnite) = i
922 : end if
923 : end do
924 :
925 : ! Associate pointers to physics buffer fields
926 61920 : itim_old = pbuf_old_tim_idx()
927 61920 : if (cldfsnow_idx > 0) then
928 247680 : call pbuf_get_field(pbuf, cldfsnow_idx, cldfsnow, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
929 : endif
930 61920 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
931 0 : call pbuf_get_field(pbuf, cldfgrau_idx, cldfgrau, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
932 : endif
933 :
934 247680 : call pbuf_get_field(pbuf, cld_idx, cld, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
935 :
936 61920 : call pbuf_get_field(pbuf, qrs_idx, qrs)
937 61920 : call pbuf_get_field(pbuf, qrl_idx, qrl)
938 :
939 61920 : call pbuf_get_field(pbuf, fsnt_idx, fsnt)
940 61920 : call pbuf_get_field(pbuf, fsds_idx, fsds)
941 61920 : call pbuf_get_field(pbuf, fsns_idx, fsns)
942 61920 : call pbuf_get_field(pbuf, flns_idx, flns)
943 61920 : call pbuf_get_field(pbuf, flnt_idx, flnt)
944 :
945 61920 : if (spectralflux) then
946 0 : call pbuf_get_field(pbuf, su_idx, su)
947 0 : call pbuf_get_field(pbuf, sd_idx, sd)
948 0 : call pbuf_get_field(pbuf, lu_idx, lu)
949 0 : call pbuf_get_field(pbuf, ld_idx, ld)
950 : end if
951 :
952 : ! For CRM, make cloud equal to input observations:
953 61920 : if (scm_crm_mode .and. have_cld) then
954 0 : do k = 1, pver
955 0 : cld(:ncol,k)= cldobs(k)
956 : end do
957 : end if
958 :
959 : ! Find tropopause height if needed for diagnostic output
960 61920 : if (hist_fld_active('FSNR') .or. hist_fld_active('FLNR')) then
961 0 : call tropopause_find(state, troplev, tropP=p_trop, primary=TROP_ALG_HYBSTOB, backup=TROP_ALG_CLIMATE)
962 : endif
963 :
964 : ! Get time of next radiation calculation - albedos will need to be
965 : ! calculated by each surface model at this time
966 61920 : nextsw_cday = radiation_nextsw_cday()
967 :
968 61920 : if (dosw .or. dolw) then
969 :
970 : ! construct an RRTMG state object
971 30960 : r_state = rrtmg_state_create( state, cam_in )
972 :
973 30960 : call t_startf('cldoptics')
974 :
975 30960 : if (cldfsnow_idx > 0) then
976 2910240 : do k = 1, pver
977 48108240 : do i = 1, ncol
978 48077280 : cldfprime(i,k) = max(cld(i,k), cldfsnow(i,k))
979 : end do
980 : end do
981 : else
982 0 : cldfprime(:ncol,:) = cld(:ncol,:)
983 : end if
984 :
985 30960 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
986 0 : do k = 1, pver
987 0 : do i = 1, ncol
988 0 : cldfprime(i,k) = max(cld(i,k), cldfgrau(i,k))
989 : end do
990 : end do
991 : end if
992 :
993 30960 : if (dosw) then
994 :
995 30960 : if (oldcldoptics) then
996 0 : call ec_ice_optics_sw(state, pbuf, ice_tau, ice_tau_w, ice_tau_w_g, ice_tau_w_f, oldicewp=.false.)
997 0 : call slingo_liq_optics_sw(state, pbuf, liq_tau, liq_tau_w, liq_tau_w_g, liq_tau_w_f, oldliqwp=.false.)
998 : else
999 0 : select case (icecldoptics)
1000 : case ('ebertcurry')
1001 0 : call ec_ice_optics_sw(state, pbuf, ice_tau, ice_tau_w, ice_tau_w_g, ice_tau_w_f, oldicewp=.true.)
1002 : case ('mitchell')
1003 30960 : call get_ice_optics_sw(state, pbuf, ice_tau, ice_tau_w, ice_tau_w_g, ice_tau_w_f)
1004 : case default
1005 30960 : call endrun('iccldoptics must be one either ebertcurry or mitchell')
1006 : end select
1007 :
1008 0 : select case (liqcldoptics)
1009 : case ('slingo')
1010 0 : call slingo_liq_optics_sw(state, pbuf, liq_tau, liq_tau_w, liq_tau_w_g, liq_tau_w_f, oldliqwp=.true.)
1011 : case ('gammadist')
1012 30960 : call get_liquid_optics_sw(state, pbuf, liq_tau, liq_tau_w, liq_tau_w_g, liq_tau_w_f)
1013 : case default
1014 30960 : call endrun('liqcldoptics must be either slingo or gammadist')
1015 : end select
1016 : end if
1017 :
1018 680880240 : cld_tau(:,:ncol,:) = liq_tau(:,:ncol,:) + ice_tau(:,:ncol,:)
1019 680880240 : cld_tau_w(:,:ncol,:) = liq_tau_w(:,:ncol,:) + ice_tau_w(:,:ncol,:)
1020 680880240 : cld_tau_w_g(:,:ncol,:) = liq_tau_w_g(:,:ncol,:) + ice_tau_w_g(:,:ncol,:)
1021 680880240 : cld_tau_w_f(:,:ncol,:) = liq_tau_w_f(:,:ncol,:) + ice_tau_w_f(:,:ncol,:)
1022 :
1023 30960 : if (cldfsnow_idx > 0) then
1024 : ! add in snow
1025 30960 : call get_snow_optics_sw(state, pbuf, snow_tau, snow_tau_w, snow_tau_w_g, snow_tau_w_f)
1026 516960 : do i = 1, ncol
1027 45714960 : do k = 1, pver
1028 :
1029 45684000 : if (cldfprime(i,k) > 0._r8) then
1030 :
1031 0 : c_cld_tau(:,i,k) = ( cldfsnow(i,k)*snow_tau(:,i,k) &
1032 124975605 : + cld(i,k)*cld_tau(:,i,k) )/cldfprime(i,k)
1033 :
1034 : c_cld_tau_w(:,i,k) = ( cldfsnow(i,k)*snow_tau_w(:,i,k) &
1035 124975605 : + cld(i,k)*cld_tau_w(:,i,k) )/cldfprime(i,k)
1036 :
1037 : c_cld_tau_w_g(:,i,k) = ( cldfsnow(i,k)*snow_tau_w_g(:,i,k) &
1038 124975605 : + cld(i,k)*cld_tau_w_g(:,i,k) )/cldfprime(i,k)
1039 :
1040 : c_cld_tau_w_f(:,i,k) = ( cldfsnow(i,k)*snow_tau_w_f(:,i,k) &
1041 124975605 : + cld(i,k)*cld_tau_w_f(:,i,k) )/cldfprime(i,k)
1042 : else
1043 552994395 : c_cld_tau(:,i,k) = 0._r8
1044 552994395 : c_cld_tau_w(:,i,k) = 0._r8
1045 552994395 : c_cld_tau_w_g(:,i,k) = 0._r8
1046 552994395 : c_cld_tau_w_f(:,i,k) = 0._r8
1047 : end if
1048 : end do
1049 : end do
1050 : else
1051 0 : c_cld_tau(:,:ncol,:) = cld_tau(:,:ncol,:)
1052 0 : c_cld_tau_w(:,:ncol,:) = cld_tau_w(:,:ncol,:)
1053 0 : c_cld_tau_w_g(:,:ncol,:) = cld_tau_w_g(:,:ncol,:)
1054 0 : c_cld_tau_w_f(:,:ncol,:) = cld_tau_w_f(:,:ncol,:)
1055 : end if
1056 :
1057 30960 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
1058 : ! add in graupel
1059 0 : call get_grau_optics_sw(state, pbuf, grau_tau, grau_tau_w, grau_tau_w_g, grau_tau_w_f)
1060 0 : do i = 1, ncol
1061 0 : do k = 1, pver
1062 :
1063 0 : if (cldfprime(i,k) > 0._r8) then
1064 :
1065 0 : c_cld_tau(:,i,k) = ( cldfgrau(i,k)*grau_tau(:,i,k) &
1066 0 : + cld(i,k)*c_cld_tau(:,i,k) )/cldfprime(i,k)
1067 :
1068 : c_cld_tau_w(:,i,k) = ( cldfgrau(i,k)*grau_tau_w(:,i,k) &
1069 0 : + cld(i,k)*c_cld_tau_w(:,i,k) )/cldfprime(i,k)
1070 :
1071 : c_cld_tau_w_g(:,i,k) = ( cldfgrau(i,k)*grau_tau_w_g(:,i,k) &
1072 0 : + cld(i,k)*c_cld_tau_w_g(:,i,k) )/cldfprime(i,k)
1073 :
1074 : c_cld_tau_w_f(:,i,k) = ( cldfgrau(i,k)*grau_tau_w_f(:,i,k) &
1075 0 : + cld(i,k)*c_cld_tau_w_f(:,i,k) )/cldfprime(i,k)
1076 : else
1077 0 : c_cld_tau(:,i,k) = 0._r8
1078 0 : c_cld_tau_w(:,i,k) = 0._r8
1079 0 : c_cld_tau_w_g(:,i,k) = 0._r8
1080 0 : c_cld_tau_w_f(:,i,k) = 0._r8
1081 : end if
1082 : end do
1083 : end do
1084 : end if
1085 :
1086 : ! Output cloud optical depth fields for the visible band
1087 48108240 : rd%tot_icld_vistau(:ncol,:) = c_cld_tau(idx_sw_diag,:ncol,:)
1088 48108240 : rd%liq_icld_vistau(:ncol,:) = liq_tau(idx_sw_diag,:ncol,:)
1089 48108240 : rd%ice_icld_vistau(:ncol,:) = ice_tau(idx_sw_diag,:ncol,:)
1090 :
1091 30960 : if (cldfsnow_idx > 0) then
1092 48108240 : rd%snow_icld_vistau(:ncol,:) = snow_tau(idx_sw_diag,:ncol,:)
1093 : endif
1094 :
1095 30960 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
1096 0 : rd%grau_icld_vistau(:ncol,:) = grau_tau(idx_sw_diag,:ncol,:)
1097 : endif
1098 :
1099 : ! multiply by total cloud fraction to get gridbox value
1100 48108240 : rd%tot_cld_vistau(:ncol,:) = c_cld_tau(idx_sw_diag,:ncol,:)*cldfprime(:ncol,:)
1101 :
1102 : ! add fillvalue for night columns
1103 273960 : do i = 1, Nnite
1104 22842000 : rd%tot_cld_vistau(IdxNite(i),:) = fillvalue
1105 22842000 : rd%tot_icld_vistau(IdxNite(i),:) = fillvalue
1106 22842000 : rd%liq_icld_vistau(IdxNite(i),:) = fillvalue
1107 22842000 : rd%ice_icld_vistau(IdxNite(i),:) = fillvalue
1108 243000 : if (cldfsnow_idx > 0) then
1109 22842000 : rd%snow_icld_vistau(IdxNite(i),:) = fillvalue
1110 : end if
1111 273960 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
1112 0 : rd%grau_icld_vistau(IdxNite(i),:) = fillvalue
1113 : end if
1114 : end do
1115 :
1116 30960 : if (write_output) call radiation_output_cld(lchnk, ncol, rd)
1117 :
1118 : end if ! if (dosw)
1119 :
1120 30960 : if (dolw) then
1121 :
1122 30960 : if (oldcldoptics) then
1123 0 : call cloud_rad_props_get_lw(state, pbuf, cld_lw_abs, oldcloud=.true.)
1124 : else
1125 0 : select case (icecldoptics)
1126 : case ('ebertcurry')
1127 0 : call ec_ice_get_rad_props_lw(state, pbuf, ice_lw_abs, oldicewp=.true.)
1128 : case ('mitchell')
1129 30960 : call ice_cloud_get_rad_props_lw(state, pbuf, ice_lw_abs)
1130 : case default
1131 30960 : call endrun('iccldoptics must be one either ebertcurry or mitchell')
1132 : end select
1133 :
1134 0 : select case (liqcldoptics)
1135 : case ('slingo')
1136 0 : call slingo_liq_get_rad_props_lw(state, pbuf, liq_lw_abs, oldliqwp=.true.)
1137 : case ('gammadist')
1138 30960 : call liquid_cloud_get_rad_props_lw(state, pbuf, liq_lw_abs)
1139 : case default
1140 30960 : call endrun('liqcldoptics must be either slingo or gammadist')
1141 : end select
1142 :
1143 771276240 : cld_lw_abs(:,:ncol,:) = liq_lw_abs(:,:ncol,:) + ice_lw_abs(:,:ncol,:)
1144 :
1145 : end if
1146 :
1147 30960 : if (cldfsnow_idx > 0) then
1148 :
1149 : ! add in snow
1150 30960 : call snow_cloud_get_rad_props_lw(state, pbuf, snow_lw_abs)
1151 :
1152 516960 : do i = 1, ncol
1153 45714960 : do k = 1, pver
1154 45684000 : if (cldfprime(i,k) > 0._r8) then
1155 0 : c_cld_lw_abs(:,i,k) = ( cldfsnow(i,k)*snow_lw_abs(:,i,k) &
1156 141639019 : + cld(i,k)*cld_lw_abs(:,i,k) )/cldfprime(i,k)
1157 : else
1158 626726981 : c_cld_lw_abs(:,i,k) = 0._r8
1159 : end if
1160 : end do
1161 : end do
1162 : else
1163 0 : c_cld_lw_abs(:,:ncol,:) = cld_lw_abs(:,:ncol,:)
1164 : end if
1165 :
1166 30960 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
1167 :
1168 : ! add in graupel
1169 0 : call grau_cloud_get_rad_props_lw(state, pbuf, grau_lw_abs)
1170 :
1171 0 : do i = 1, ncol
1172 0 : do k = 1, pver
1173 0 : if (cldfprime(i,k) > 0._r8) then
1174 0 : c_cld_lw_abs(:,i,k) = ( cldfgrau(i,k)*grau_lw_abs(:,i,k) &
1175 0 : + cld(i,k)*c_cld_lw_abs(:,i,k) )/cldfprime(i,k)
1176 : else
1177 0 : c_cld_lw_abs(:,i,k) = 0._r8
1178 : end if
1179 : end do
1180 : end do
1181 : end if
1182 :
1183 : end if ! if (dolw)
1184 :
1185 30960 : call t_stopf('cldoptics')
1186 :
1187 : ! Solar radiation computation
1188 :
1189 30960 : if (dosw) then
1190 :
1191 30960 : call get_variability(sfac)
1192 :
1193 : ! Get the active climate/diagnostic shortwave calculations
1194 30960 : call rad_cnst_get_call_list(active_calls)
1195 :
1196 : ! The climate (icall==0) calculation must occur last.
1197 371520 : do icall = N_DIAG, 0, -1
1198 :
1199 371520 : if (active_calls(icall)) then
1200 :
1201 : ! update the concentrations in the RRTMG state object
1202 30960 : call rrtmg_state_update(state, pbuf, icall, r_state)
1203 :
1204 : call aer_rad_props_sw(icall, state, pbuf, nnite, idxnite, &
1205 30960 : aer_tau, aer_tau_w, aer_tau_w_g, aer_tau_w_f)
1206 :
1207 48108240 : rd%cld_tau_cloudsim(:ncol,:) = cld_tau(rrtmg_sw_cloudsim_band,:ncol,:)
1208 48625200 : rd%aer_tau550(:ncol,:) = aer_tau(:ncol,:,idx_sw_diag)
1209 48625200 : rd%aer_tau400(:ncol,:) = aer_tau(:ncol,:,idx_sw_diag+1)
1210 48625200 : rd%aer_tau700(:ncol,:) = aer_tau(:ncol,:,idx_sw_diag-1)
1211 :
1212 : call rad_rrtmg_sw( &
1213 : lchnk, ncol, num_rrtmg_levs, r_state, state%pmid, &
1214 : cldfprime, aer_tau, aer_tau_w, aer_tau_w_g, aer_tau_w_f, &
1215 : eccf, coszrs, rd%solin, sfac, cam_in%asdir, &
1216 : cam_in%asdif, cam_in%aldir, cam_in%aldif, qrs, rd%qrsc, &
1217 : fsnt, rd%fsntc, rd%fsntoa, rd%fsutoa, rd%fsntoac, &
1218 : rd%fsnirt, rd%fsnrtc, rd%fsnirtsq, fsns, rd%fsnsc, &
1219 : rd%fsdsc, fsds, cam_out%sols, cam_out%soll, cam_out%solsd, &
1220 : cam_out%solld, fns, fcns, Nday, Nnite, &
1221 : IdxDay, IdxNite, su, sd, E_cld_tau=c_cld_tau, &
1222 : E_cld_tau_w=c_cld_tau_w, E_cld_tau_w_g=c_cld_tau_w_g, &
1223 30960 : E_cld_tau_w_f=c_cld_tau_w_f, old_convert=.false.)
1224 :
1225 : ! Output net fluxes at 200 mb
1226 30960 : call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fcns, rd%fsn200c)
1227 30960 : call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fns, rd%fsn200)
1228 30960 : if (hist_fld_active('FSNR')) then
1229 0 : do i = 1,ncol
1230 0 : call vertinterp(1, 1, pverp, state%pint(i,:), p_trop(i), fns(i,:), rd%fsnr(i))
1231 : end do
1232 : end if
1233 :
1234 30960 : if (write_output) call radiation_output_sw(lchnk, ncol, icall, rd, pbuf, cam_out)
1235 :
1236 : end if
1237 : end do
1238 :
1239 : end if
1240 :
1241 : ! Output aerosol mmr
1242 30960 : call rad_cnst_out(0, state, pbuf)
1243 :
1244 : ! Longwave radiation computation
1245 :
1246 30960 : if (dolw) then
1247 :
1248 30960 : call rad_cnst_get_call_list(active_calls)
1249 :
1250 : ! The climate (icall==0) calculation must occur last.
1251 371520 : do icall = N_DIAG, 0, -1
1252 :
1253 371520 : if (active_calls(icall)) then
1254 :
1255 : ! update the conctrations in the RRTMG state object
1256 30960 : call rrtmg_state_update( state, pbuf, icall, r_state)
1257 :
1258 30960 : call aer_rad_props_lw(icall, state, pbuf, aer_lw_abs)
1259 :
1260 : call rad_rrtmg_lw( &
1261 : lchnk, ncol, num_rrtmg_levs, r_state, state%pmid, &
1262 : aer_lw_abs, cldfprime, c_cld_lw_abs, qrl, rd%qrlc, &
1263 : flns, flnt, rd%flnsc, rd%flntc, cam_out%flwds, &
1264 : rd%flut, rd%flutc, fnl, fcnl, rd%fldsc, &
1265 30960 : lu, ld)
1266 :
1267 : ! Output fluxes at 200 mb
1268 30960 : call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fnl, rd%fln200)
1269 30960 : call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fcnl, rd%fln200c)
1270 30960 : if (hist_fld_active('FLNR')) then
1271 0 : do i = 1,ncol
1272 0 : call vertinterp(1, 1, pverp, state%pint(i,:), p_trop(i), fnl(i,:), rd%flnr(i))
1273 : end do
1274 : end if
1275 :
1276 30960 : flntclr(:) = 0._r8
1277 30960 : freqclr(:) = 0._r8
1278 516960 : do i = 1, ncol
1279 46200960 : if (maxval(cldfprime(i,:)) <= 0.1_r8) then
1280 135602 : freqclr(i) = 1._r8
1281 135602 : flntclr(i) = rd%flntc(i)
1282 : end if
1283 : end do
1284 :
1285 30960 : if (write_output) call radiation_output_lw(lchnk, ncol, icall, rd, pbuf, cam_out, freqclr, flntclr)
1286 :
1287 : end if
1288 : end do
1289 :
1290 : end if
1291 :
1292 : ! deconstruct the RRTMG state object
1293 30960 : call rrtmg_state_destroy(r_state)
1294 :
1295 30960 : if (docosp) then
1296 :
1297 : ! initialize and calculate emis
1298 0 : emis(:,:) = 0._r8
1299 0 : emis(:ncol,:) = 1._r8 - exp(-cld_lw_abs(rrtmg_lw_cloudsim_band,:ncol,:))
1300 0 : call outfld('EMIS', emis, pcols, lchnk)
1301 :
1302 : ! compute grid-box mean SW and LW snow optical depth for use by COSP
1303 0 : gb_snow_tau(:,:) = 0._r8
1304 0 : gb_snow_lw(:,:) = 0._r8
1305 0 : if (cldfsnow_idx > 0) then
1306 0 : do i = 1, ncol
1307 0 : do k = 1, pver
1308 0 : if (cldfsnow(i,k) > 0._r8) then
1309 :
1310 : ! Add graupel to snow tau for cosp
1311 0 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
1312 0 : gb_snow_tau(i,k) = snow_tau(rrtmg_sw_cloudsim_band,i,k)*cldfsnow(i,k) + &
1313 0 : grau_tau(rrtmg_sw_cloudsim_band,i,k)*cldfgrau(i,k)
1314 : gb_snow_lw(i,k) = snow_lw_abs(rrtmg_lw_cloudsim_band,i,k)*cldfsnow(i,k) + &
1315 0 : grau_lw_abs(rrtmg_lw_cloudsim_band,i,k)*cldfgrau(i,k)
1316 : else
1317 0 : gb_snow_tau(i,k) = snow_tau(rrtmg_sw_cloudsim_band,i,k)*cldfsnow(i,k)
1318 0 : gb_snow_lw(i,k) = snow_lw_abs(rrtmg_lw_cloudsim_band,i,k)*cldfsnow(i,k)
1319 : end if
1320 : end if
1321 : end do
1322 : end do
1323 : end if
1324 :
1325 : ! advance counter for this timestep (chunk dimension required for thread safety)
1326 0 : cosp_cnt(lchnk) = cosp_cnt(lchnk) + 1
1327 :
1328 : ! if counter is the same as cosp_nradsteps, run cosp and reset counter
1329 0 : if (cosp_nradsteps .eq. cosp_cnt(lchnk)) then
1330 :
1331 : ! N.B.: For snow optical properties, the GRID-BOX MEAN shortwave and longwave
1332 : ! optical depths are passed.
1333 : call cospsimulator_intr_run(state, pbuf, cam_in, emis, coszrs, &
1334 : cld_swtau_in=cld_tau(rrtmg_sw_cloudsim_band,:,:),&
1335 0 : snow_tau_in=gb_snow_tau, snow_emis_in=gb_snow_lw)
1336 0 : cosp_cnt(lchnk) = 0
1337 : end if
1338 : end if
1339 :
1340 : else ! if (dosw .or. dolw) then
1341 :
1342 : ! convert radiative heating rates from Q*dp to Q for energy conservation
1343 2910240 : do k =1 , pver
1344 48108240 : do i = 1, ncol
1345 45198000 : qrs(i,k) = qrs(i,k)/state%pdel(i,k)
1346 48077280 : qrl(i,k) = qrl(i,k)/state%pdel(i,k)
1347 : end do
1348 : end do
1349 :
1350 : end if ! if (dosw .or. dolw) then
1351 :
1352 : ! output rad inputs and resulting heating rates
1353 61920 : call rad_data_write( pbuf, state, cam_in, coszrs )
1354 :
1355 : ! Compute net radiative heating tendency
1356 : call radheat_tend(state, pbuf, ptend, qrl, qrs, fsns, &
1357 61920 : fsnt, flns, flnt, cam_in%asdir, net_flx)
1358 :
1359 61920 : if (write_output) then
1360 : ! Compute heating rate for dtheta/dt
1361 5820480 : do k = 1, pver
1362 96216480 : do i = 1, ncol
1363 96154560 : ftem(i,k) = (qrs(i,k) + qrl(i,k))/cpair * (1.e5_r8/state%pmid(i,k))**cappa
1364 : end do
1365 : end do
1366 61920 : call outfld('HR', ftem, pcols, lchnk)
1367 : end if
1368 :
1369 : ! convert radiative heating rates to Q*dp for energy conservation
1370 5820480 : do k = 1, pver
1371 96216480 : do i = 1, ncol
1372 90396000 : qrs(i,k) = qrs(i,k)*state%pdel(i,k)
1373 96154560 : qrl(i,k) = qrl(i,k)*state%pdel(i,k)
1374 : end do
1375 : end do
1376 :
1377 1033920 : cam_out%netsw(:ncol) = fsns(:ncol)
1378 :
1379 61920 : if (.not. present(rd_out)) then
1380 61920 : deallocate(rd)
1381 : end if
1382 :
1383 123840 : end subroutine radiation_tend
1384 :
1385 : !===============================================================================
1386 :
1387 30960 : subroutine radiation_output_sw(lchnk, ncol, icall, rd, pbuf, cam_out)
1388 :
1389 : ! Dump shortwave radiation information to history buffer.
1390 :
1391 : integer , intent(in) :: lchnk
1392 : integer, intent(in) :: ncol
1393 : integer, intent(in) :: icall
1394 : type(rad_out_t), intent(in) :: rd
1395 : type(physics_buffer_desc), pointer :: pbuf(:)
1396 : type(cam_out_t), intent(in) :: cam_out
1397 :
1398 : ! local variables
1399 30960 : real(r8), pointer :: qrs(:,:)
1400 30960 : real(r8), pointer :: fsnt(:)
1401 30960 : real(r8), pointer :: fsns(:)
1402 30960 : real(r8), pointer :: fsds(:)
1403 :
1404 : real(r8) :: ftem(pcols)
1405 : !----------------------------------------------------------------------------
1406 :
1407 30960 : call pbuf_get_field(pbuf, qrs_idx, qrs)
1408 30960 : call pbuf_get_field(pbuf, fsnt_idx, fsnt)
1409 30960 : call pbuf_get_field(pbuf, fsns_idx, fsns)
1410 30960 : call pbuf_get_field(pbuf, fsds_idx, fsds)
1411 :
1412 30960 : call outfld('SOLIN'//diag(icall), rd%solin, pcols, lchnk)
1413 :
1414 48108240 : call outfld('QRS'//diag(icall), qrs(:ncol,:)/cpair, ncol, lchnk)
1415 48108240 : call outfld('QRSC'//diag(icall), rd%qrsc(:ncol,:)/cpair, ncol, lchnk)
1416 :
1417 30960 : call outfld('FSNT'//diag(icall), fsnt, pcols, lchnk)
1418 30960 : call outfld('FSNTC'//diag(icall), rd%fsntc, pcols, lchnk)
1419 30960 : call outfld('FSNTOA'//diag(icall), rd%fsntoa, pcols, lchnk)
1420 30960 : call outfld('FSNTOAC'//diag(icall), rd%fsntoac, pcols, lchnk)
1421 :
1422 516960 : ftem(:ncol) = rd%fsntoa(:ncol) - rd%fsntoac(:ncol)
1423 30960 : call outfld('SWCF'//diag(icall), ftem, pcols, lchnk)
1424 :
1425 30960 : call outfld('FSUTOA'//diag(icall), rd%fsutoa, pcols, lchnk)
1426 :
1427 30960 : call outfld('FSNIRTOA'//diag(icall), rd%fsnirt, pcols, lchnk)
1428 30960 : call outfld('FSNRTOAC'//diag(icall), rd%fsnrtc, pcols, lchnk)
1429 30960 : call outfld('FSNRTOAS'//diag(icall), rd%fsnirtsq, pcols, lchnk)
1430 :
1431 30960 : call outfld('FSN200'//diag(icall), rd%fsn200, pcols, lchnk)
1432 30960 : call outfld('FSN200C'//diag(icall), rd%fsn200c, pcols, lchnk)
1433 :
1434 30960 : call outfld('FSNR'//diag(icall), rd%fsnr, pcols, lchnk)
1435 :
1436 30960 : call outfld('SOLS'//diag(icall), cam_out%sols, pcols, lchnk)
1437 30960 : call outfld('SOLL'//diag(icall), cam_out%soll, pcols, lchnk)
1438 30960 : call outfld('SOLSD'//diag(icall), cam_out%solsd, pcols, lchnk)
1439 30960 : call outfld('SOLLD'//diag(icall), cam_out%solld, pcols, lchnk)
1440 :
1441 30960 : call outfld('FSNS'//diag(icall), fsns, pcols, lchnk)
1442 30960 : call outfld('FSNSC'//diag(icall), rd%fsnsc, pcols, lchnk)
1443 :
1444 30960 : call outfld('FSDS'//diag(icall), fsds, pcols, lchnk)
1445 30960 : call outfld('FSDSC'//diag(icall), rd%fsdsc, pcols, lchnk)
1446 :
1447 92880 : end subroutine radiation_output_sw
1448 :
1449 :
1450 : !===============================================================================
1451 :
1452 30960 : subroutine radiation_output_cld(lchnk, ncol, rd)
1453 :
1454 : ! Dump shortwave cloud optics information to history buffer.
1455 :
1456 : integer , intent(in) :: lchnk
1457 : integer, intent(in) :: ncol
1458 : type(rad_out_t), intent(in) :: rd
1459 : !----------------------------------------------------------------------------
1460 :
1461 30960 : call outfld('TOT_CLD_VISTAU', rd%tot_cld_vistau, pcols, lchnk)
1462 30960 : call outfld('TOT_ICLD_VISTAU', rd%tot_icld_vistau, pcols, lchnk)
1463 30960 : call outfld('LIQ_ICLD_VISTAU', rd%liq_icld_vistau, pcols, lchnk)
1464 30960 : call outfld('ICE_ICLD_VISTAU', rd%ice_icld_vistau, pcols, lchnk)
1465 30960 : if (cldfsnow_idx > 0) then
1466 30960 : call outfld('SNOW_ICLD_VISTAU', rd%snow_icld_vistau, pcols, lchnk)
1467 : endif
1468 30960 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
1469 0 : call outfld('GRAU_ICLD_VISTAU', rd%grau_icld_vistau, pcols, lchnk)
1470 : endif
1471 30960 : end subroutine radiation_output_cld
1472 :
1473 : !===============================================================================
1474 :
1475 30960 : subroutine radiation_output_lw(lchnk, ncol, icall, rd, pbuf, cam_out, freqclr, flntclr)
1476 :
1477 : ! Dump longwave radiation information to history buffer
1478 :
1479 : integer, intent(in) :: lchnk
1480 : integer, intent(in) :: ncol
1481 : integer, intent(in) :: icall ! icall=0 for climate diagnostics
1482 : type(rad_out_t), intent(in) :: rd
1483 : type(physics_buffer_desc), pointer :: pbuf(:)
1484 : type(cam_out_t), intent(in) :: cam_out
1485 : real(r8), intent(in) :: freqclr(pcols)
1486 : real(r8), intent(in) :: flntclr(pcols)
1487 :
1488 : ! local variables
1489 30960 : real(r8), pointer :: qrl(:,:)
1490 30960 : real(r8), pointer :: flnt(:)
1491 30960 : real(r8), pointer :: flns(:)
1492 :
1493 : real(r8) :: ftem(pcols)
1494 : !----------------------------------------------------------------------------
1495 :
1496 30960 : call pbuf_get_field(pbuf, qrl_idx, qrl)
1497 30960 : call pbuf_get_field(pbuf, flnt_idx, flnt)
1498 30960 : call pbuf_get_field(pbuf, flns_idx, flns)
1499 :
1500 48108240 : call outfld('QRL'//diag(icall), qrl(:ncol,:)/cpair, ncol, lchnk)
1501 48108240 : call outfld('QRLC'//diag(icall), rd%qrlc(:ncol,:)/cpair, ncol, lchnk)
1502 :
1503 30960 : call outfld('FLNT'//diag(icall), flnt, pcols, lchnk)
1504 30960 : call outfld('FLNTC'//diag(icall), rd%flntc, pcols, lchnk)
1505 :
1506 30960 : call outfld('FREQCLR'//diag(icall), freqclr, pcols, lchnk)
1507 30960 : call outfld('FLNTCLR'//diag(icall), flntclr, pcols, lchnk)
1508 :
1509 30960 : call outfld('FLUT'//diag(icall), rd%flut, pcols, lchnk)
1510 30960 : call outfld('FLUTC'//diag(icall), rd%flutc, pcols, lchnk)
1511 :
1512 516960 : ftem(:ncol) = rd%flutc(:ncol) - rd%flut(:ncol)
1513 30960 : call outfld('LWCF'//diag(icall), ftem, pcols, lchnk)
1514 :
1515 30960 : call outfld('FLN200'//diag(icall), rd%fln200, pcols, lchnk)
1516 30960 : call outfld('FLN200C'//diag(icall), rd%fln200c, pcols, lchnk)
1517 :
1518 30960 : call outfld('FLNR'//diag(icall), rd%flnr, pcols, lchnk)
1519 :
1520 30960 : call outfld('FLNS'//diag(icall), flns, pcols, lchnk)
1521 30960 : call outfld('FLNSC'//diag(icall), rd%flnsc, pcols, lchnk)
1522 :
1523 30960 : call outfld('FLDS'//diag(icall), cam_out%flwds, pcols, lchnk)
1524 30960 : call outfld('FLDSC'//diag(icall), rd%fldsc, pcols, lchnk)
1525 :
1526 30960 : end subroutine radiation_output_lw
1527 :
1528 : !===============================================================================
1529 :
1530 : subroutine calc_col_mean(state, mmr_pointer, mean_value)
1531 :
1532 : ! Compute the column mean mass mixing ratio.
1533 :
1534 : type(physics_state), intent(in) :: state
1535 : real(r8), dimension(:,:), pointer :: mmr_pointer ! mass mixing ratio (lev)
1536 : real(r8), dimension(pcols), intent(out) :: mean_value ! column mean mmr
1537 :
1538 : integer :: i, k, ncol
1539 : real(r8) :: ptot(pcols)
1540 : !-----------------------------------------------------------------------
1541 :
1542 : ncol = state%ncol
1543 : mean_value = 0.0_r8
1544 : ptot = 0.0_r8
1545 :
1546 : do k=1,pver
1547 : do i=1,ncol
1548 : mean_value(i) = mean_value(i) + mmr_pointer(i,k)*state%pdeldry(i,k)
1549 : ptot(i) = ptot(i) + state%pdeldry(i,k)
1550 : end do
1551 : end do
1552 : do i=1,ncol
1553 : mean_value(i) = mean_value(i) / ptot(i)
1554 : end do
1555 :
1556 : end subroutine calc_col_mean
1557 :
1558 : !===============================================================================
1559 :
1560 0 : end module radiation
|
