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 288000 : 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 288000 : if (present(timestep)) then
290 126720 : nstep = timestep
291 : else
292 161280 : nstep = get_nstep()
293 : end if
294 :
295 207360 : 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 207360 : .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 80640 : .or. nstep <= irad_always
306 :
307 : case default
308 288000 : call endrun('radiation_do: unknown operation:'//op)
309 :
310 : end select
311 1536 : end function radiation_do
312 :
313 : !================================================================================================
314 :
315 80640 : 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 80640 : radiation_nextsw_cday = -1._r8
329 80640 : dosw = .false.
330 80640 : nstep = get_nstep()
331 80640 : dtime = get_step_size()
332 80640 : offset = 0
333 : do while (.not. dosw)
334 126720 : nstep = nstep + 1
335 126720 : offset = offset + dtime
336 126720 : if (radiation_do('sw', nstep)) then
337 80640 : radiation_nextsw_cday = get_curr_calday(offset=offset)
338 : dosw = .true.
339 : end if
340 : end do
341 80640 : 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 80640 : if (get_nstep() >= 1) then
347 76800 : caldayp1 = get_curr_calday(offset=int(dtime))
348 76800 : if (caldayp1 /= radiation_nextsw_cday) radiation_nextsw_cday = -1._r8
349 : end if
350 :
351 80640 : 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 4608 : cosp_cnt(begchunk:endchunk) = cosp_cnt_init
432 : else
433 4608 : 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 31691520 : 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_cam, 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 80640 : real(r8), pointer :: cld(:,:) ! cloud fraction
785 80640 : real(r8), pointer :: cldfsnow(:,:) ! cloud fraction of just "snow clouds- whatever they are"
786 80640 : real(r8), pointer :: cldfgrau(:,:) ! cloud fraction of just "snow clouds- whatever they are"
787 80640 : real(r8), pointer :: qrs(:,:) ! shortwave radiative heating rate
788 80640 : real(r8), pointer :: qrl(:,:) ! longwave radiative heating rate
789 80640 : real(r8), pointer :: fsds(:) ! Surface solar down flux
790 80640 : real(r8), pointer :: fsns(:) ! Surface solar absorbed flux
791 80640 : real(r8), pointer :: fsnt(:) ! Net column abs solar flux at model top
792 80640 : real(r8), pointer :: flns(:) ! Srf longwave cooling (up-down) flux
793 80640 : 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 80640 : 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 80640 : lchnk = state%lchnk
881 80640 : ncol = state%ncol
882 :
883 80640 : if (present(rd_out)) then
884 : rd => rd_out
885 : write_output = .false.
886 : else
887 80640 : allocate(rd)
888 : write_output=.true.
889 : end if
890 :
891 80640 : dosw = radiation_do('sw') ! do shortwave heating calc this timestep?
892 80640 : dolw = radiation_do('lw') ! do longwave heating calc this timestep?
893 :
894 : ! Cosine solar zenith angle for current time step
895 80640 : calday = get_curr_calday()
896 80640 : call get_rlat_all_p(lchnk, ncol, clat)
897 80640 : call get_rlon_all_p(lchnk, ncol, clon)
898 :
899 : call shr_orb_decl(calday, eccen, mvelpp, lambm0, obliqr, &
900 80640 : delta, eccf)
901 :
902 80640 : 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 1241856 : do i = 1, ncol
908 1241856 : 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 80640 : Nday = 0
914 80640 : Nnite = 0
915 1241856 : do i = 1, ncol
916 1241856 : if ( coszrs(i) > 0.0_r8 ) then
917 580608 : Nday = Nday + 1
918 580608 : IdxDay(Nday) = i
919 : else
920 580608 : Nnite = Nnite + 1
921 580608 : IdxNite(Nnite) = i
922 : end if
923 : end do
924 :
925 : ! Associate pointers to physics buffer fields
926 80640 : itim_old = pbuf_old_tim_idx()
927 80640 : if (cldfsnow_idx > 0) then
928 322560 : call pbuf_get_field(pbuf, cldfsnow_idx, cldfsnow, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
929 : endif
930 80640 : 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 322560 : call pbuf_get_field(pbuf, cld_idx, cld, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
935 :
936 80640 : call pbuf_get_field(pbuf, qrs_idx, qrs)
937 80640 : call pbuf_get_field(pbuf, qrl_idx, qrl)
938 :
939 80640 : call pbuf_get_field(pbuf, fsnt_idx, fsnt)
940 80640 : call pbuf_get_field(pbuf, fsds_idx, fsds)
941 80640 : call pbuf_get_field(pbuf, fsns_idx, fsns)
942 80640 : call pbuf_get_field(pbuf, flns_idx, flns)
943 80640 : call pbuf_get_field(pbuf, flnt_idx, flnt)
944 :
945 80640 : 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 80640 : 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 80640 : if (hist_fld_active('FSNR') .or. hist_fld_active('FLNR')) then
961 : !REMOVECAM - no longer need this when CAM is retired and pcols no longer exists
962 0 : troplev(:) = 0
963 0 : p_trop(:) = 0._r8
964 : !REMOVECAM_END
965 0 : call tropopause_find_cam(state, troplev, tropP=p_trop, primary=TROP_ALG_HYBSTOB, backup=TROP_ALG_CLIMATE)
966 : endif
967 :
968 : ! Get time of next radiation calculation - albedos will need to be
969 : ! calculated by each surface model at this time
970 80640 : nextsw_cday = radiation_nextsw_cday()
971 :
972 80640 : if (dosw .or. dolw) then
973 :
974 : ! construct an RRTMG state object
975 38400 : r_state = rrtmg_state_create( state, cam_in )
976 :
977 38400 : call t_startf('cldoptics')
978 :
979 38400 : if (cldfsnow_idx > 0) then
980 1267200 : do k = 1, pver
981 18961920 : do i = 1, ncol
982 18923520 : cldfprime(i,k) = max(cld(i,k), cldfsnow(i,k))
983 : end do
984 : end do
985 : else
986 0 : cldfprime(:ncol,:) = cld(:ncol,:)
987 : end if
988 :
989 38400 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
990 0 : do k = 1, pver
991 0 : do i = 1, ncol
992 0 : cldfprime(i,k) = max(cld(i,k), cldfgrau(i,k))
993 : end do
994 : end do
995 : end if
996 :
997 38400 : if (dosw) then
998 :
999 38400 : if (oldcldoptics) then
1000 0 : call ec_ice_optics_sw(state, pbuf, ice_tau, ice_tau_w, ice_tau_w_g, ice_tau_w_f, oldicewp=.false.)
1001 0 : call slingo_liq_optics_sw(state, pbuf, liq_tau, liq_tau_w, liq_tau_w_g, liq_tau_w_f, oldliqwp=.false.)
1002 : else
1003 0 : select case (icecldoptics)
1004 : case ('ebertcurry')
1005 0 : call ec_ice_optics_sw(state, pbuf, ice_tau, ice_tau_w, ice_tau_w_g, ice_tau_w_f, oldicewp=.true.)
1006 : case ('mitchell')
1007 38400 : call get_ice_optics_sw(state, pbuf, ice_tau, ice_tau_w, ice_tau_w_g, ice_tau_w_f)
1008 : case default
1009 38400 : call endrun('iccldoptics must be one either ebertcurry or mitchell')
1010 : end select
1011 :
1012 0 : select case (liqcldoptics)
1013 : case ('slingo')
1014 0 : call slingo_liq_optics_sw(state, pbuf, liq_tau, liq_tau_w, liq_tau_w_g, liq_tau_w_f, oldliqwp=.true.)
1015 : case ('gammadist')
1016 38400 : call get_liquid_optics_sw(state, pbuf, liq_tau, liq_tau_w, liq_tau_w_g, liq_tau_w_f)
1017 : case default
1018 38400 : call endrun('liqcldoptics must be either slingo or gammadist')
1019 : end select
1020 : end if
1021 :
1022 266688000 : cld_tau(:,:ncol,:) = liq_tau(:,:ncol,:) + ice_tau(:,:ncol,:)
1023 266688000 : cld_tau_w(:,:ncol,:) = liq_tau_w(:,:ncol,:) + ice_tau_w(:,:ncol,:)
1024 266688000 : cld_tau_w_g(:,:ncol,:) = liq_tau_w_g(:,:ncol,:) + ice_tau_w_g(:,:ncol,:)
1025 266688000 : cld_tau_w_f(:,:ncol,:) = liq_tau_w_f(:,:ncol,:) + ice_tau_w_f(:,:ncol,:)
1026 :
1027 38400 : if (cldfsnow_idx > 0) then
1028 : ! add in snow
1029 38400 : call get_snow_optics_sw(state, pbuf, snow_tau, snow_tau_w, snow_tau_w_g, snow_tau_w_f)
1030 591360 : do i = 1, ncol
1031 18286080 : do k = 1, pver
1032 :
1033 18247680 : if (cldfprime(i,k) > 0._r8) then
1034 :
1035 0 : c_cld_tau(:,i,k) = ( cldfsnow(i,k)*snow_tau(:,i,k) &
1036 87645285 : + cld(i,k)*cld_tau(:,i,k) )/cldfprime(i,k)
1037 :
1038 : c_cld_tau_w(:,i,k) = ( cldfsnow(i,k)*snow_tau_w(:,i,k) &
1039 87645285 : + cld(i,k)*cld_tau_w(:,i,k) )/cldfprime(i,k)
1040 :
1041 : c_cld_tau_w_g(:,i,k) = ( cldfsnow(i,k)*snow_tau_w_g(:,i,k) &
1042 87645285 : + cld(i,k)*cld_tau_w_g(:,i,k) )/cldfprime(i,k)
1043 :
1044 : c_cld_tau_w_f(:,i,k) = ( cldfsnow(i,k)*snow_tau_w_f(:,i,k) &
1045 87645285 : + cld(i,k)*cld_tau_w_f(:,i,k) )/cldfprime(i,k)
1046 : else
1047 177775515 : c_cld_tau(:,i,k) = 0._r8
1048 177775515 : c_cld_tau_w(:,i,k) = 0._r8
1049 177775515 : c_cld_tau_w_g(:,i,k) = 0._r8
1050 177775515 : c_cld_tau_w_f(:,i,k) = 0._r8
1051 : end if
1052 : end do
1053 : end do
1054 : else
1055 0 : c_cld_tau(:,:ncol,:) = cld_tau(:,:ncol,:)
1056 0 : c_cld_tau_w(:,:ncol,:) = cld_tau_w(:,:ncol,:)
1057 0 : c_cld_tau_w_g(:,:ncol,:) = cld_tau_w_g(:,:ncol,:)
1058 0 : c_cld_tau_w_f(:,:ncol,:) = cld_tau_w_f(:,:ncol,:)
1059 : end if
1060 :
1061 38400 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
1062 : ! add in graupel
1063 0 : call get_grau_optics_sw(state, pbuf, grau_tau, grau_tau_w, grau_tau_w_g, grau_tau_w_f)
1064 0 : do i = 1, ncol
1065 0 : do k = 1, pver
1066 :
1067 0 : if (cldfprime(i,k) > 0._r8) then
1068 :
1069 0 : c_cld_tau(:,i,k) = ( cldfgrau(i,k)*grau_tau(:,i,k) &
1070 0 : + cld(i,k)*c_cld_tau(:,i,k) )/cldfprime(i,k)
1071 :
1072 : c_cld_tau_w(:,i,k) = ( cldfgrau(i,k)*grau_tau_w(:,i,k) &
1073 0 : + cld(i,k)*c_cld_tau_w(:,i,k) )/cldfprime(i,k)
1074 :
1075 : c_cld_tau_w_g(:,i,k) = ( cldfgrau(i,k)*grau_tau_w_g(:,i,k) &
1076 0 : + cld(i,k)*c_cld_tau_w_g(:,i,k) )/cldfprime(i,k)
1077 :
1078 : c_cld_tau_w_f(:,i,k) = ( cldfgrau(i,k)*grau_tau_w_f(:,i,k) &
1079 0 : + cld(i,k)*c_cld_tau_w_f(:,i,k) )/cldfprime(i,k)
1080 : else
1081 0 : c_cld_tau(:,i,k) = 0._r8
1082 0 : c_cld_tau_w(:,i,k) = 0._r8
1083 0 : c_cld_tau_w_g(:,i,k) = 0._r8
1084 0 : c_cld_tau_w_f(:,i,k) = 0._r8
1085 : end if
1086 : end do
1087 : end do
1088 : end if
1089 :
1090 : ! Output cloud optical depth fields for the visible band
1091 18961920 : rd%tot_icld_vistau(:ncol,:) = c_cld_tau(idx_sw_diag,:ncol,:)
1092 18961920 : rd%liq_icld_vistau(:ncol,:) = liq_tau(idx_sw_diag,:ncol,:)
1093 18961920 : rd%ice_icld_vistau(:ncol,:) = ice_tau(idx_sw_diag,:ncol,:)
1094 :
1095 38400 : if (cldfsnow_idx > 0) then
1096 18961920 : rd%snow_icld_vistau(:ncol,:) = snow_tau(idx_sw_diag,:ncol,:)
1097 : endif
1098 :
1099 38400 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
1100 0 : rd%grau_icld_vistau(:ncol,:) = grau_tau(idx_sw_diag,:ncol,:)
1101 : endif
1102 :
1103 : ! multiply by total cloud fraction to get gridbox value
1104 18961920 : rd%tot_cld_vistau(:ncol,:) = c_cld_tau(idx_sw_diag,:ncol,:)*cldfprime(:ncol,:)
1105 :
1106 : ! add fillvalue for night columns
1107 314880 : do i = 1, Nnite
1108 9123840 : rd%tot_cld_vistau(IdxNite(i),:) = fillvalue
1109 9123840 : rd%tot_icld_vistau(IdxNite(i),:) = fillvalue
1110 9123840 : rd%liq_icld_vistau(IdxNite(i),:) = fillvalue
1111 9123840 : rd%ice_icld_vistau(IdxNite(i),:) = fillvalue
1112 276480 : if (cldfsnow_idx > 0) then
1113 9123840 : rd%snow_icld_vistau(IdxNite(i),:) = fillvalue
1114 : end if
1115 314880 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
1116 0 : rd%grau_icld_vistau(IdxNite(i),:) = fillvalue
1117 : end if
1118 : end do
1119 :
1120 38400 : if (write_output) call radiation_output_cld(lchnk, ncol, rd)
1121 :
1122 : end if ! if (dosw)
1123 :
1124 38400 : if (dolw) then
1125 :
1126 38400 : if (oldcldoptics) then
1127 0 : call cloud_rad_props_get_lw(state, pbuf, cld_lw_abs, oldcloud=.true.)
1128 : else
1129 0 : select case (icecldoptics)
1130 : case ('ebertcurry')
1131 0 : call ec_ice_get_rad_props_lw(state, pbuf, ice_lw_abs, oldicewp=.true.)
1132 : case ('mitchell')
1133 38400 : call ice_cloud_get_rad_props_lw(state, pbuf, ice_lw_abs)
1134 : case default
1135 38400 : call endrun('iccldoptics must be one either ebertcurry or mitchell')
1136 : end select
1137 :
1138 0 : select case (liqcldoptics)
1139 : case ('slingo')
1140 0 : call slingo_liq_get_rad_props_lw(state, pbuf, liq_lw_abs, oldliqwp=.true.)
1141 : case ('gammadist')
1142 38400 : call liquid_cloud_get_rad_props_lw(state, pbuf, liq_lw_abs)
1143 : case default
1144 38400 : call endrun('liqcldoptics must be either slingo or gammadist')
1145 : end select
1146 :
1147 302077440 : cld_lw_abs(:,:ncol,:) = liq_lw_abs(:,:ncol,:) + ice_lw_abs(:,:ncol,:)
1148 :
1149 : end if
1150 :
1151 38400 : if (cldfsnow_idx > 0) then
1152 :
1153 : ! add in snow
1154 38400 : call snow_cloud_get_rad_props_lw(state, pbuf, snow_lw_abs)
1155 :
1156 591360 : do i = 1, ncol
1157 18286080 : do k = 1, pver
1158 18247680 : if (cldfprime(i,k) > 0._r8) then
1159 0 : c_cld_lw_abs(:,i,k) = ( cldfsnow(i,k)*snow_lw_abs(:,i,k) &
1160 99331323 : + cld(i,k)*cld_lw_abs(:,i,k) )/cldfprime(i,k)
1161 : else
1162 201478917 : c_cld_lw_abs(:,i,k) = 0._r8
1163 : end if
1164 : end do
1165 : end do
1166 : else
1167 0 : c_cld_lw_abs(:,:ncol,:) = cld_lw_abs(:,:ncol,:)
1168 : end if
1169 :
1170 38400 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
1171 :
1172 : ! add in graupel
1173 0 : call grau_cloud_get_rad_props_lw(state, pbuf, grau_lw_abs)
1174 :
1175 0 : do i = 1, ncol
1176 0 : do k = 1, pver
1177 0 : if (cldfprime(i,k) > 0._r8) then
1178 0 : c_cld_lw_abs(:,i,k) = ( cldfgrau(i,k)*grau_lw_abs(:,i,k) &
1179 0 : + cld(i,k)*c_cld_lw_abs(:,i,k) )/cldfprime(i,k)
1180 : else
1181 0 : c_cld_lw_abs(:,i,k) = 0._r8
1182 : end if
1183 : end do
1184 : end do
1185 : end if
1186 :
1187 : end if ! if (dolw)
1188 :
1189 38400 : call t_stopf('cldoptics')
1190 :
1191 : ! Solar radiation computation
1192 :
1193 38400 : if (dosw) then
1194 :
1195 38400 : call get_variability(sfac)
1196 :
1197 : ! Get the active climate/diagnostic shortwave calculations
1198 38400 : call rad_cnst_get_call_list(active_calls)
1199 :
1200 : ! The climate (icall==0) calculation must occur last.
1201 460800 : do icall = N_DIAG, 0, -1
1202 :
1203 460800 : if (active_calls(icall)) then
1204 :
1205 : ! update the concentrations in the RRTMG state object
1206 38400 : call rrtmg_state_update(state, pbuf, icall, r_state)
1207 :
1208 : call aer_rad_props_sw(icall, state, pbuf, nnite, idxnite, &
1209 38400 : aer_tau, aer_tau_w, aer_tau_w_g, aer_tau_w_f)
1210 :
1211 18961920 : rd%cld_tau_cloudsim(:ncol,:) = cld_tau(rrtmg_sw_cloudsim_band,:ncol,:)
1212 19553280 : rd%aer_tau550(:ncol,:) = aer_tau(:ncol,:,idx_sw_diag)
1213 19553280 : rd%aer_tau400(:ncol,:) = aer_tau(:ncol,:,idx_sw_diag+1)
1214 19553280 : rd%aer_tau700(:ncol,:) = aer_tau(:ncol,:,idx_sw_diag-1)
1215 :
1216 : call rad_rrtmg_sw( &
1217 : lchnk, ncol, num_rrtmg_levs, r_state, state%pmid, &
1218 : cldfprime, aer_tau, aer_tau_w, aer_tau_w_g, aer_tau_w_f, &
1219 : eccf, coszrs, rd%solin, sfac, cam_in%asdir, &
1220 : cam_in%asdif, cam_in%aldir, cam_in%aldif, qrs, rd%qrsc, &
1221 : fsnt, rd%fsntc, rd%fsntoa, rd%fsutoa, rd%fsntoac, &
1222 : rd%fsnirt, rd%fsnrtc, rd%fsnirtsq, fsns, rd%fsnsc, &
1223 : rd%fsdsc, fsds, cam_out%sols, cam_out%soll, cam_out%solsd, &
1224 : cam_out%solld, fns, fcns, Nday, Nnite, &
1225 : IdxDay, IdxNite, su, sd, E_cld_tau=c_cld_tau, &
1226 : E_cld_tau_w=c_cld_tau_w, E_cld_tau_w_g=c_cld_tau_w_g, &
1227 38400 : E_cld_tau_w_f=c_cld_tau_w_f, old_convert=.false.)
1228 :
1229 : ! Output net fluxes at 200 mb
1230 38400 : call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fcns, rd%fsn200c)
1231 38400 : call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fns, rd%fsn200)
1232 38400 : if (hist_fld_active('FSNR')) then
1233 0 : do i = 1,ncol
1234 0 : call vertinterp(1, 1, pverp, state%pint(i,:), p_trop(i), fns(i,:), rd%fsnr(i))
1235 : end do
1236 : end if
1237 :
1238 38400 : if (write_output) call radiation_output_sw(lchnk, ncol, icall, rd, pbuf, cam_out)
1239 :
1240 : end if
1241 : end do
1242 :
1243 : end if
1244 :
1245 : ! Output aerosol mmr
1246 38400 : call rad_cnst_out(0, state, pbuf)
1247 :
1248 : ! Longwave radiation computation
1249 :
1250 38400 : if (dolw) then
1251 :
1252 38400 : call rad_cnst_get_call_list(active_calls)
1253 :
1254 : ! The climate (icall==0) calculation must occur last.
1255 460800 : do icall = N_DIAG, 0, -1
1256 :
1257 460800 : if (active_calls(icall)) then
1258 :
1259 : ! update the conctrations in the RRTMG state object
1260 38400 : call rrtmg_state_update( state, pbuf, icall, r_state)
1261 :
1262 38400 : call aer_rad_props_lw(icall, state, pbuf, aer_lw_abs)
1263 :
1264 : call rad_rrtmg_lw( &
1265 : lchnk, ncol, num_rrtmg_levs, r_state, state%pmid, &
1266 : aer_lw_abs, cldfprime, c_cld_lw_abs, qrl, rd%qrlc, &
1267 : flns, flnt, rd%flnsc, rd%flntc, cam_out%flwds, &
1268 : rd%flut, rd%flutc, fnl, fcnl, rd%fldsc, &
1269 38400 : lu, ld)
1270 :
1271 : ! Output fluxes at 200 mb
1272 38400 : call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fnl, rd%fln200)
1273 38400 : call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fcnl, rd%fln200c)
1274 38400 : if (hist_fld_active('FLNR')) then
1275 0 : do i = 1,ncol
1276 0 : call vertinterp(1, 1, pverp, state%pint(i,:), p_trop(i), fnl(i,:), rd%flnr(i))
1277 : end do
1278 : end if
1279 :
1280 38400 : flntclr(:) = 0._r8
1281 38400 : freqclr(:) = 0._r8
1282 591360 : do i = 1, ncol
1283 18839040 : if (maxval(cldfprime(i,:)) <= 0.1_r8) then
1284 73307 : freqclr(i) = 1._r8
1285 73307 : flntclr(i) = rd%flntc(i)
1286 : end if
1287 : end do
1288 :
1289 38400 : if (write_output) call radiation_output_lw(lchnk, ncol, icall, rd, pbuf, cam_out, freqclr, flntclr)
1290 :
1291 : end if
1292 : end do
1293 :
1294 : end if
1295 :
1296 : ! deconstruct the RRTMG state object
1297 38400 : call rrtmg_state_destroy(r_state)
1298 :
1299 38400 : if (docosp) then
1300 :
1301 : ! initialize and calculate emis
1302 0 : emis(:,:) = 0._r8
1303 0 : emis(:ncol,:) = 1._r8 - exp(-cld_lw_abs(rrtmg_lw_cloudsim_band,:ncol,:))
1304 0 : call outfld('EMIS', emis, pcols, lchnk)
1305 :
1306 : ! compute grid-box mean SW and LW snow optical depth for use by COSP
1307 0 : gb_snow_tau(:,:) = 0._r8
1308 0 : gb_snow_lw(:,:) = 0._r8
1309 0 : if (cldfsnow_idx > 0) then
1310 0 : do i = 1, ncol
1311 0 : do k = 1, pver
1312 0 : if (cldfsnow(i,k) > 0._r8) then
1313 :
1314 : ! Add graupel to snow tau for cosp
1315 0 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
1316 0 : gb_snow_tau(i,k) = snow_tau(rrtmg_sw_cloudsim_band,i,k)*cldfsnow(i,k) + &
1317 0 : grau_tau(rrtmg_sw_cloudsim_band,i,k)*cldfgrau(i,k)
1318 : gb_snow_lw(i,k) = snow_lw_abs(rrtmg_lw_cloudsim_band,i,k)*cldfsnow(i,k) + &
1319 0 : grau_lw_abs(rrtmg_lw_cloudsim_band,i,k)*cldfgrau(i,k)
1320 : else
1321 0 : gb_snow_tau(i,k) = snow_tau(rrtmg_sw_cloudsim_band,i,k)*cldfsnow(i,k)
1322 0 : gb_snow_lw(i,k) = snow_lw_abs(rrtmg_lw_cloudsim_band,i,k)*cldfsnow(i,k)
1323 : end if
1324 : end if
1325 : end do
1326 : end do
1327 : end if
1328 :
1329 : ! advance counter for this timestep (chunk dimension required for thread safety)
1330 0 : cosp_cnt(lchnk) = cosp_cnt(lchnk) + 1
1331 :
1332 : ! if counter is the same as cosp_nradsteps, run cosp and reset counter
1333 0 : if (cosp_nradsteps .eq. cosp_cnt(lchnk)) then
1334 :
1335 : ! N.B.: For snow optical properties, the GRID-BOX MEAN shortwave and longwave
1336 : ! optical depths are passed.
1337 : call cospsimulator_intr_run(state, pbuf, cam_in, emis, coszrs, &
1338 : cld_swtau_in=cld_tau(rrtmg_sw_cloudsim_band,:,:),&
1339 0 : snow_tau_in=gb_snow_tau, snow_emis_in=gb_snow_lw)
1340 0 : cosp_cnt(lchnk) = 0
1341 : end if
1342 : end if
1343 :
1344 : else ! if (dosw .or. dolw) then
1345 :
1346 : ! convert radiative heating rates from Q*dp to Q for energy conservation
1347 1393920 : do k =1 , pver
1348 20858112 : do i = 1, ncol
1349 19464192 : qrs(i,k) = qrs(i,k)/state%pdel(i,k)
1350 20815872 : qrl(i,k) = qrl(i,k)/state%pdel(i,k)
1351 : end do
1352 : end do
1353 :
1354 : end if ! if (dosw .or. dolw) then
1355 :
1356 : ! output rad inputs and resulting heating rates
1357 80640 : call rad_data_write( pbuf, state, cam_in, coszrs )
1358 :
1359 : ! Compute net radiative heating tendency
1360 : call radheat_tend(state, pbuf, ptend, qrl, qrs, fsns, &
1361 80640 : fsnt, flns, flnt, cam_in%asdir, net_flx)
1362 :
1363 80640 : if (write_output) then
1364 : ! Compute heating rate for dtheta/dt
1365 2661120 : do k = 1, pver
1366 39820032 : do i = 1, ncol
1367 39739392 : ftem(i,k) = (qrs(i,k) + qrl(i,k))/cpair * (1.e5_r8/state%pmid(i,k))**cappa
1368 : end do
1369 : end do
1370 80640 : call outfld('HR', ftem, pcols, lchnk)
1371 : end if
1372 :
1373 : ! convert radiative heating rates to Q*dp for energy conservation
1374 2661120 : do k = 1, pver
1375 39820032 : do i = 1, ncol
1376 37158912 : qrs(i,k) = qrs(i,k)*state%pdel(i,k)
1377 39739392 : qrl(i,k) = qrl(i,k)*state%pdel(i,k)
1378 : end do
1379 : end do
1380 :
1381 1241856 : cam_out%netsw(:ncol) = fsns(:ncol)
1382 :
1383 80640 : if (.not. present(rd_out)) then
1384 80640 : deallocate(rd)
1385 : end if
1386 :
1387 161280 : end subroutine radiation_tend
1388 :
1389 : !===============================================================================
1390 :
1391 38400 : subroutine radiation_output_sw(lchnk, ncol, icall, rd, pbuf, cam_out)
1392 :
1393 : ! Dump shortwave radiation information to history buffer.
1394 :
1395 : integer , intent(in) :: lchnk
1396 : integer, intent(in) :: ncol
1397 : integer, intent(in) :: icall
1398 : type(rad_out_t), intent(in) :: rd
1399 : type(physics_buffer_desc), pointer :: pbuf(:)
1400 : type(cam_out_t), intent(in) :: cam_out
1401 :
1402 : ! local variables
1403 38400 : real(r8), pointer :: qrs(:,:)
1404 38400 : real(r8), pointer :: fsnt(:)
1405 38400 : real(r8), pointer :: fsns(:)
1406 38400 : real(r8), pointer :: fsds(:)
1407 :
1408 : real(r8) :: ftem(pcols)
1409 : !----------------------------------------------------------------------------
1410 :
1411 38400 : call pbuf_get_field(pbuf, qrs_idx, qrs)
1412 38400 : call pbuf_get_field(pbuf, fsnt_idx, fsnt)
1413 38400 : call pbuf_get_field(pbuf, fsns_idx, fsns)
1414 38400 : call pbuf_get_field(pbuf, fsds_idx, fsds)
1415 :
1416 38400 : call outfld('SOLIN'//diag(icall), rd%solin, pcols, lchnk)
1417 :
1418 18961920 : call outfld('QRS'//diag(icall), qrs(:ncol,:)/cpair, ncol, lchnk)
1419 18961920 : call outfld('QRSC'//diag(icall), rd%qrsc(:ncol,:)/cpair, ncol, lchnk)
1420 :
1421 38400 : call outfld('FSNT'//diag(icall), fsnt, pcols, lchnk)
1422 38400 : call outfld('FSNTC'//diag(icall), rd%fsntc, pcols, lchnk)
1423 38400 : call outfld('FSNTOA'//diag(icall), rd%fsntoa, pcols, lchnk)
1424 38400 : call outfld('FSNTOAC'//diag(icall), rd%fsntoac, pcols, lchnk)
1425 :
1426 591360 : ftem(:ncol) = rd%fsntoa(:ncol) - rd%fsntoac(:ncol)
1427 38400 : call outfld('SWCF'//diag(icall), ftem, pcols, lchnk)
1428 :
1429 38400 : call outfld('FSUTOA'//diag(icall), rd%fsutoa, pcols, lchnk)
1430 :
1431 38400 : call outfld('FSNIRTOA'//diag(icall), rd%fsnirt, pcols, lchnk)
1432 38400 : call outfld('FSNRTOAC'//diag(icall), rd%fsnrtc, pcols, lchnk)
1433 38400 : call outfld('FSNRTOAS'//diag(icall), rd%fsnirtsq, pcols, lchnk)
1434 :
1435 38400 : call outfld('FSN200'//diag(icall), rd%fsn200, pcols, lchnk)
1436 38400 : call outfld('FSN200C'//diag(icall), rd%fsn200c, pcols, lchnk)
1437 :
1438 38400 : call outfld('FSNR'//diag(icall), rd%fsnr, pcols, lchnk)
1439 :
1440 38400 : call outfld('SOLS'//diag(icall), cam_out%sols, pcols, lchnk)
1441 38400 : call outfld('SOLL'//diag(icall), cam_out%soll, pcols, lchnk)
1442 38400 : call outfld('SOLSD'//diag(icall), cam_out%solsd, pcols, lchnk)
1443 38400 : call outfld('SOLLD'//diag(icall), cam_out%solld, pcols, lchnk)
1444 :
1445 38400 : call outfld('FSNS'//diag(icall), fsns, pcols, lchnk)
1446 38400 : call outfld('FSNSC'//diag(icall), rd%fsnsc, pcols, lchnk)
1447 :
1448 38400 : call outfld('FSDS'//diag(icall), fsds, pcols, lchnk)
1449 38400 : call outfld('FSDSC'//diag(icall), rd%fsdsc, pcols, lchnk)
1450 :
1451 119040 : end subroutine radiation_output_sw
1452 :
1453 :
1454 : !===============================================================================
1455 :
1456 38400 : subroutine radiation_output_cld(lchnk, ncol, rd)
1457 :
1458 : ! Dump shortwave cloud optics information to history buffer.
1459 :
1460 : integer , intent(in) :: lchnk
1461 : integer, intent(in) :: ncol
1462 : type(rad_out_t), intent(in) :: rd
1463 : !----------------------------------------------------------------------------
1464 :
1465 38400 : call outfld('TOT_CLD_VISTAU', rd%tot_cld_vistau, pcols, lchnk)
1466 38400 : call outfld('TOT_ICLD_VISTAU', rd%tot_icld_vistau, pcols, lchnk)
1467 38400 : call outfld('LIQ_ICLD_VISTAU', rd%liq_icld_vistau, pcols, lchnk)
1468 38400 : call outfld('ICE_ICLD_VISTAU', rd%ice_icld_vistau, pcols, lchnk)
1469 38400 : if (cldfsnow_idx > 0) then
1470 38400 : call outfld('SNOW_ICLD_VISTAU', rd%snow_icld_vistau, pcols, lchnk)
1471 : endif
1472 38400 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
1473 0 : call outfld('GRAU_ICLD_VISTAU', rd%grau_icld_vistau, pcols, lchnk)
1474 : endif
1475 38400 : end subroutine radiation_output_cld
1476 :
1477 : !===============================================================================
1478 :
1479 38400 : subroutine radiation_output_lw(lchnk, ncol, icall, rd, pbuf, cam_out, freqclr, flntclr)
1480 :
1481 : ! Dump longwave radiation information to history buffer
1482 :
1483 : integer, intent(in) :: lchnk
1484 : integer, intent(in) :: ncol
1485 : integer, intent(in) :: icall ! icall=0 for climate diagnostics
1486 : type(rad_out_t), intent(in) :: rd
1487 : type(physics_buffer_desc), pointer :: pbuf(:)
1488 : type(cam_out_t), intent(in) :: cam_out
1489 : real(r8), intent(in) :: freqclr(pcols)
1490 : real(r8), intent(in) :: flntclr(pcols)
1491 :
1492 : ! local variables
1493 38400 : real(r8), pointer :: qrl(:,:)
1494 38400 : real(r8), pointer :: flnt(:)
1495 38400 : real(r8), pointer :: flns(:)
1496 :
1497 : real(r8) :: ftem(pcols)
1498 : !----------------------------------------------------------------------------
1499 :
1500 38400 : call pbuf_get_field(pbuf, qrl_idx, qrl)
1501 38400 : call pbuf_get_field(pbuf, flnt_idx, flnt)
1502 38400 : call pbuf_get_field(pbuf, flns_idx, flns)
1503 :
1504 18961920 : call outfld('QRL'//diag(icall), qrl(:ncol,:)/cpair, ncol, lchnk)
1505 18961920 : call outfld('QRLC'//diag(icall), rd%qrlc(:ncol,:)/cpair, ncol, lchnk)
1506 :
1507 38400 : call outfld('FLNT'//diag(icall), flnt, pcols, lchnk)
1508 38400 : call outfld('FLNTC'//diag(icall), rd%flntc, pcols, lchnk)
1509 :
1510 38400 : call outfld('FREQCLR'//diag(icall), freqclr, pcols, lchnk)
1511 38400 : call outfld('FLNTCLR'//diag(icall), flntclr, pcols, lchnk)
1512 :
1513 38400 : call outfld('FLUT'//diag(icall), rd%flut, pcols, lchnk)
1514 38400 : call outfld('FLUTC'//diag(icall), rd%flutc, pcols, lchnk)
1515 :
1516 591360 : ftem(:ncol) = rd%flutc(:ncol) - rd%flut(:ncol)
1517 38400 : call outfld('LWCF'//diag(icall), ftem, pcols, lchnk)
1518 :
1519 38400 : call outfld('FLN200'//diag(icall), rd%fln200, pcols, lchnk)
1520 38400 : call outfld('FLN200C'//diag(icall), rd%fln200c, pcols, lchnk)
1521 :
1522 38400 : call outfld('FLNR'//diag(icall), rd%flnr, pcols, lchnk)
1523 :
1524 38400 : call outfld('FLNS'//diag(icall), flns, pcols, lchnk)
1525 38400 : call outfld('FLNSC'//diag(icall), rd%flnsc, pcols, lchnk)
1526 :
1527 38400 : call outfld('FLDS'//diag(icall), cam_out%flwds, pcols, lchnk)
1528 38400 : call outfld('FLDSC'//diag(icall), rd%fldsc, pcols, lchnk)
1529 :
1530 38400 : end subroutine radiation_output_lw
1531 :
1532 : !===============================================================================
1533 :
1534 : subroutine calc_col_mean(state, mmr_pointer, mean_value)
1535 :
1536 : ! Compute the column mean mass mixing ratio.
1537 :
1538 : type(physics_state), intent(in) :: state
1539 : real(r8), dimension(:,:), pointer :: mmr_pointer ! mass mixing ratio (lev)
1540 : real(r8), dimension(pcols), intent(out) :: mean_value ! column mean mmr
1541 :
1542 : integer :: i, k, ncol
1543 : real(r8) :: ptot(pcols)
1544 : !-----------------------------------------------------------------------
1545 :
1546 : ncol = state%ncol
1547 : mean_value = 0.0_r8
1548 : ptot = 0.0_r8
1549 :
1550 : do k=1,pver
1551 : do i=1,ncol
1552 : mean_value(i) = mean_value(i) + mmr_pointer(i,k)*state%pdeldry(i,k)
1553 : ptot(i) = ptot(i) + state%pdeldry(i,k)
1554 : end do
1555 : end do
1556 : do i=1,ncol
1557 : mean_value(i) = mean_value(i) / ptot(i)
1558 : end do
1559 :
1560 : end subroutine calc_col_mean
1561 :
1562 : !===============================================================================
1563 :
1564 0 : end module radiation
|
