Line data Source code
1 : module radiation
2 :
3 : !---------------------------------------------------------------------------------
4 : !
5 : ! CAM interface to RRTMGP radiation parameterization.
6 : !
7 : !---------------------------------------------------------------------------------
8 :
9 : use shr_kind_mod, only: r8=>shr_kind_r8, cl=>shr_kind_cl
10 : use spmd_utils, only: masterproc
11 : use ppgrid, only: pcols, pver, pverp, begchunk, endchunk
12 : use ref_pres, only: pref_edge
13 : use physics_types, only: physics_state, physics_ptend
14 : use phys_control, only: phys_getopts
15 : use physics_buffer, only: physics_buffer_desc, pbuf_add_field, dtype_r8, pbuf_get_index, &
16 : pbuf_set_field, pbuf_get_field, pbuf_old_tim_idx
17 : use camsrfexch, only: cam_out_t, cam_in_t
18 : use physconst, only: cappa, cpair, gravit
19 :
20 : use time_manager, only: get_nstep, is_first_step, is_first_restart_step, &
21 : get_curr_calday, get_step_size
22 :
23 : use rad_constituents, only: N_DIAG, rad_cnst_get_call_list, rad_cnst_get_gas, rad_cnst_out
24 :
25 : use rrtmgp_inputs, only: rrtmgp_inputs_init
26 :
27 : use radconstants, only: nradgas, gasnamelength, gaslist, nswbands, nlwbands, &
28 : nswgpts, set_wavenumber_bands
29 : use rad_solar_var, only: rad_solar_var_init, get_variability
30 :
31 : use cloud_rad_props, only: cloud_rad_props_init
32 :
33 : use cospsimulator_intr, only: docosp, cospsimulator_intr_init, &
34 : cospsimulator_intr_run, cosp_nradsteps
35 :
36 : use scamMod, only: scm_crm_mode, single_column, have_cld, cldobs
37 :
38 : use cam_history, only: addfld, add_default, horiz_only, outfld, hist_fld_active
39 :
40 : use radiation_data, only: rad_data_register, rad_data_init
41 :
42 : use ioFileMod, only: getfil
43 : use cam_pio_utils, only: cam_pio_openfile
44 : use pio, only: file_desc_t, var_desc_t, &
45 : pio_int, pio_double, PIO_NOERR, &
46 : pio_seterrorhandling, PIO_BCAST_ERROR, &
47 : pio_inq_dimlen, pio_inq_dimid, pio_inq_varid, &
48 : pio_def_var, pio_put_var, pio_get_var, &
49 : pio_put_att, PIO_NOWRITE, pio_closefile
50 :
51 : use mo_gas_concentrations, only: ty_gas_concs
52 : use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp
53 : use mo_optical_props, only: ty_optical_props_1scl, ty_optical_props_2str
54 : use mo_source_functions, only: ty_source_func_lw
55 : use mo_fluxes, only: ty_fluxes_broadband
56 : use mo_fluxes_byband, only: ty_fluxes_byband
57 :
58 : use string_utils, only: to_lower
59 : use cam_abortutils, only: endrun, handle_allocate_error
60 : use cam_logfile, only: iulog
61 :
62 :
63 : implicit none
64 : private
65 : save
66 :
67 : public :: &
68 : radiation_readnl, &! read namelist variables
69 : radiation_register, &! registers radiation physics buffer fields
70 : radiation_do, &! query which radiation calcs are done this timestep
71 : radiation_init, &! initialization
72 : radiation_define_restart, &! define variables for restart
73 : radiation_write_restart, &! write variables to restart
74 : radiation_read_restart, &! read variables from restart
75 : radiation_tend, &! compute heating rates and fluxes
76 : rad_out_t ! type for diagnostic outputs
77 :
78 : integer,public, allocatable :: cosp_cnt(:) ! counter for cosp
79 : integer,public :: cosp_cnt_init = 0 !initial value for cosp counter
80 :
81 : real(r8), public, protected :: nextsw_cday ! future radiation calday for surface models
82 :
83 : type rad_out_t
84 : real(r8) :: solin(pcols) ! Solar incident flux
85 :
86 : real(r8) :: qrsc(pcols,pver)
87 :
88 : real(r8) :: fsnsc(pcols) ! Clear sky surface abs solar flux
89 : real(r8) :: fsntc(pcols) ! Clear sky total column abs solar flux
90 : real(r8) :: fsdsc(pcols) ! Clear sky surface downwelling solar flux
91 :
92 : real(r8) :: fsntoa(pcols) ! Net solar flux at TOA
93 : real(r8) :: fsntoac(pcols) ! Clear sky net solar flux at TOA
94 : real(r8) :: fsutoa(pcols) ! upwelling solar flux at TOA
95 :
96 : real(r8) :: fsnirt(pcols) ! Near-IR flux absorbed at toa
97 : real(r8) :: fsnrtc(pcols) ! Clear sky near-IR flux absorbed at toa
98 : real(r8) :: fsnirtsq(pcols) ! Near-IR flux absorbed at toa >= 0.7 microns
99 :
100 : real(r8) :: fsn200(pcols) ! Net SW flux interpolated to 200 mb
101 : real(r8) :: fsn200c(pcols) ! Net clear-sky SW flux interpolated to 200 mb
102 : real(r8) :: fsnr(pcols) ! Net SW flux interpolated to tropopause
103 :
104 : real(r8) :: flux_sw_up(pcols,pverp) ! upward shortwave flux on interfaces
105 : real(r8) :: flux_sw_clr_up(pcols,pverp) ! upward shortwave clearsky flux
106 : real(r8) :: flux_sw_dn(pcols,pverp) ! downward flux
107 : real(r8) :: flux_sw_clr_dn(pcols,pverp) ! downward clearsky flux
108 :
109 : real(r8) :: flux_lw_up(pcols,pverp) ! upward longwave flux on interfaces
110 : real(r8) :: flux_lw_clr_up(pcols,pverp) ! upward longwave clearsky flux
111 : real(r8) :: flux_lw_dn(pcols,pverp) ! downward flux
112 : real(r8) :: flux_lw_clr_dn(pcols,pverp) ! downward clearsky flux
113 :
114 : real(r8) :: qrlc(pcols,pver)
115 :
116 : real(r8) :: flntc(pcols) ! Clear sky lw flux at model top
117 : real(r8) :: flut(pcols) ! Upward flux at top of model
118 : real(r8) :: flutc(pcols) ! Upward Clear Sky flux at top of model
119 : real(r8) :: lwcf(pcols) ! longwave cloud forcing
120 :
121 : real(r8) :: fln200(pcols) ! net longwave flux interpolated to 200 mb
122 : real(r8) :: fln200c(pcols) ! net clearsky longwave flux interpolated to 200 mb
123 : real(r8) :: flnr(pcols) ! net longwave flux interpolated to tropopause
124 :
125 : real(r8) :: flnsc(pcols) ! Clear sky lw flux at srf (up-down)
126 : real(r8) :: fldsc(pcols) ! Clear sky lw flux at srf (down)
127 :
128 : real(r8) :: tot_cld_vistau(pcols,pver) ! gbx water+ice cloud optical depth (only during day, night = fillvalue)
129 : real(r8) :: tot_icld_vistau(pcols,pver) ! in-cld water+ice cloud optical depth (only during day, night = fillvalue)
130 : real(r8) :: liq_icld_vistau(pcols,pver) ! in-cld liq cloud optical depth (only during day, night = fillvalue)
131 : real(r8) :: ice_icld_vistau(pcols,pver) ! in-cld ice cloud optical depth (only during day, night = fillvalue)
132 : real(r8) :: snow_icld_vistau(pcols,pver) ! snow in-cloud visible sw optical depth for output on history files
133 : real(r8) :: grau_icld_vistau(pcols,pver) ! Graupel in-cloud visible sw optical depth for output on history files
134 : end type rad_out_t
135 :
136 : ! Control variables set via namelist
137 : character(len=cl) :: coefs_lw_file ! filepath for lw coefficients
138 : character(len=cl) :: coefs_sw_file ! filepath for sw coefficients
139 :
140 : integer :: iradsw = -1 ! freq. of shortwave radiation calc in time steps (positive)
141 : ! or hours (negative).
142 : integer :: iradlw = -1 ! frequency of longwave rad. calc. in time steps (positive)
143 : ! or hours (negative).
144 :
145 : integer :: irad_always = 0 ! Specifies length of time in timesteps (positive)
146 : ! or hours (negative) SW/LW radiation will be
147 : ! run continuously from the start of an
148 : ! initial or restart run
149 : logical :: use_rad_dt_cosz = .false. ! if true, use radiation dt for all cosz calculations
150 : logical :: spectralflux = .false. ! calculate fluxes (up and down) per band.
151 : logical :: graupel_in_rad = .false. ! graupel in radiation code
152 : logical :: use_rad_uniform_angle = .false. ! if true, use the namelist rad_uniform_angle for the coszrs calculation
153 :
154 : ! active_calls is set by a rad_constituents method after parsing namelist input
155 : ! for the rad_climate and rad_diag_N entries.
156 : logical :: active_calls(0:N_DIAG)
157 :
158 : ! Physics buffer indices
159 : integer :: qrs_idx = 0
160 : integer :: qrl_idx = 0
161 : integer :: su_idx = 0
162 : integer :: sd_idx = 0
163 : integer :: lu_idx = 0
164 : integer :: ld_idx = 0
165 : integer :: fsds_idx = 0
166 : integer :: fsns_idx = 0
167 : integer :: fsnt_idx = 0
168 : integer :: flns_idx = 0
169 : integer :: flnt_idx = 0
170 : integer :: cld_idx = 0
171 : integer :: cldfsnow_idx = 0
172 : integer :: cldfgrau_idx = 0
173 :
174 : character(len=4) :: diag(0:N_DIAG) =(/' ','_d1 ','_d2 ','_d3 ','_d4 ','_d5 ',&
175 : '_d6 ','_d7 ','_d8 ','_d9 ','_d10'/)
176 :
177 : ! averaging time interval for zenith angle
178 : real(r8) :: dt_avg = 0._r8
179 : real(r8) :: rad_uniform_angle = -99._r8
180 :
181 : ! Number of layers in radiation calculations.
182 : integer :: nlay
183 :
184 : ! Number of CAM layers in radiation calculations. Is either equal to nlay, or is
185 : ! 1 less than nlay if "extra layer" is used in the radiation calculations.
186 : integer :: nlaycam
187 :
188 : ! Indices for copying data between CAM/WACCM and RRTMGP arrays. Since RRTMGP is
189 : ! vertical order agnostic we can send data using the top to bottom order used
190 : ! in CAM/WACCM. But the number of layers that RRTMGP does computations for
191 : ! may not match the number of layers in CAM/WACCM for two reasons:
192 : ! 1. If the CAM model top is below 1 Pa, then RRTMGP does calculations for an
193 : ! extra layer that is added between 1 Pa and the model top.
194 : ! 2. If the WACCM model top is above 1 Pa, then RRMTGP only does calculations
195 : ! for those model layers that are below 1 Pa.
196 : integer :: ktopcam ! Index in CAM arrays of top level (layer or interface) at which
197 : ! RRTMGP is active.
198 : integer :: ktoprad ! Index in RRTMGP arrays of the layer or interface corresponding
199 : ! to CAM's top layer or interface.
200 : ! Note: for CAM's top to bottom indexing, the index of a given layer
201 : ! (midpoint) and the upper interface of that layer, are the same.
202 :
203 : ! Gas optics objects contain the data read from the coefficients files.
204 : type(ty_gas_optics_rrtmgp) :: kdist_sw
205 : type(ty_gas_optics_rrtmgp) :: kdist_lw
206 :
207 : ! lower case version of gaslist for RRTMGP
208 : character(len=gasnamelength) :: gaslist_lc(nradgas)
209 :
210 : type(var_desc_t) :: cospcnt_desc ! cosp
211 : type(var_desc_t) :: nextsw_cday_desc
212 :
213 : !=========================================================================================
214 : contains
215 : !=========================================================================================
216 :
217 1536 : subroutine radiation_readnl(nlfile)
218 :
219 : ! Read radiation_nl namelist group.
220 :
221 : use namelist_utils, only: find_group_name
222 : use spmd_utils, only: mpicom, mstrid=>masterprocid, mpi_integer, mpi_logical, &
223 : mpi_character, mpi_real8
224 :
225 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
226 :
227 : ! Local variables
228 : integer :: unitn, ierr
229 : integer :: dtime ! timestep size
230 : character(len=32) :: errmsg
231 : character(len=*), parameter :: sub = 'radiation_readnl'
232 :
233 : character(len=cl) :: rrtmgp_coefs_lw_file, rrtmgp_coefs_sw_file
234 :
235 : namelist /radiation_nl/ &
236 : rrtmgp_coefs_lw_file, rrtmgp_coefs_sw_file, iradsw, iradlw, &
237 : irad_always, use_rad_dt_cosz, spectralflux, use_rad_uniform_angle, &
238 : rad_uniform_angle, graupel_in_rad
239 : !-----------------------------------------------------------------------------
240 :
241 1536 : if (masterproc) then
242 2 : open( newunit=unitn, file=trim(nlfile), status='old' )
243 2 : call find_group_name(unitn, 'radiation_nl', status=ierr)
244 2 : if (ierr == 0) then
245 2 : read(unitn, radiation_nl, iostat=ierr)
246 2 : if (ierr /= 0) then
247 0 : write(errmsg,'(a,i5)') 'iostat =', ierr
248 0 : call endrun(sub//': ERROR reading namelist: '//trim(errmsg))
249 : end if
250 : end if
251 2 : close(unitn)
252 : end if
253 :
254 : ! Broadcast namelist variables
255 1536 : call mpi_bcast(rrtmgp_coefs_lw_file, cl, mpi_character, mstrid, mpicom, ierr)
256 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: rrtmgp_coefs_lw_file")
257 1536 : call mpi_bcast(rrtmgp_coefs_sw_file, cl, mpi_character, mstrid, mpicom, ierr)
258 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: rrtmgp_coefs_sw_file")
259 1536 : call mpi_bcast(iradsw, 1, mpi_integer, mstrid, mpicom, ierr)
260 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: iradsw")
261 1536 : call mpi_bcast(iradlw, 1, mpi_integer, mstrid, mpicom, ierr)
262 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: iradlw")
263 1536 : call mpi_bcast(irad_always, 1, mpi_integer, mstrid, mpicom, ierr)
264 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: irad_always")
265 1536 : call mpi_bcast(use_rad_dt_cosz, 1, mpi_logical, mstrid, mpicom, ierr)
266 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: use_rad_dt_cosz")
267 1536 : call mpi_bcast(spectralflux, 1, mpi_logical, mstrid, mpicom, ierr)
268 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: spectralflux")
269 1536 : call mpi_bcast(use_rad_uniform_angle, 1, mpi_logical, mstrid, mpicom, ierr)
270 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: use_rad_uniform_angle")
271 1536 : call mpi_bcast(rad_uniform_angle, 1, mpi_real8, mstrid, mpicom, ierr)
272 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: rad_uniform_angle")
273 1536 : call mpi_bcast(graupel_in_rad, 1, mpi_logical, mstrid, mpicom, ierr)
274 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: graupel_in_rad")
275 :
276 1536 : if (use_rad_uniform_angle .and. rad_uniform_angle == -99._r8) then
277 : call endrun(sub//': ERROR - use_rad_uniform_angle is set to .true,' &
278 0 : //' but rad_uniform_angle is not set ')
279 : end if
280 :
281 : ! Set module data
282 1536 : coefs_lw_file = rrtmgp_coefs_lw_file
283 1536 : coefs_sw_file = rrtmgp_coefs_sw_file
284 :
285 : ! Convert iradsw, iradlw and irad_always from hours to timesteps if necessary
286 1536 : dtime = get_step_size()
287 1536 : if (iradsw < 0) iradsw = nint((-iradsw *3600._r8)/dtime)
288 1536 : if (iradlw < 0) iradlw = nint((-iradlw *3600._r8)/dtime)
289 1536 : if (irad_always < 0) irad_always = nint((-irad_always*3600._r8)/dtime)
290 :
291 : !-----------------------------------------------------------------------
292 : ! Print runtime options to log.
293 : !-----------------------------------------------------------------------
294 :
295 1536 : if (masterproc) then
296 2 : write(iulog,*) 'RRTMGP radiation scheme parameters:'
297 2 : write(iulog,10) trim(coefs_lw_file), trim(coefs_sw_file), nlwbands, nswbands, &
298 4 : iradsw, iradlw, irad_always, use_rad_dt_cosz, spectralflux, graupel_in_rad
299 : end if
300 :
301 : 10 format(' LW coefficents file: ', a/, &
302 : ' SW coefficents file: ', a/, &
303 : ' Number of LW bands: ',i5/, &
304 : ' Number of SW bands: ',i5/, &
305 : ' Frequency (timesteps) of Shortwave Radiation calc: ',i5/, &
306 : ' Frequency (timesteps) of Longwave Radiation calc: ',i5/, &
307 : ' SW/LW calc done every timestep for first N steps. N=',i5/, &
308 : ' Use average zenith angle: ',l5/, &
309 : ' Output spectrally resolved fluxes: ',l5/, &
310 : ' Graupel in Radiation Code: ',l5/)
311 :
312 1536 : end subroutine radiation_readnl
313 :
314 : !================================================================================================
315 :
316 1536 : subroutine radiation_register
317 :
318 : ! Register radiation fields in the physics buffer
319 :
320 1536 : call pbuf_add_field('QRS' , 'global',dtype_r8,(/pcols,pver/), qrs_idx) ! shortwave radiative heating rate
321 1536 : call pbuf_add_field('QRL' , 'global',dtype_r8,(/pcols,pver/), qrl_idx) ! longwave radiative heating rate
322 :
323 1536 : call pbuf_add_field('FSDS' , 'global',dtype_r8,(/pcols/), fsds_idx) ! Surface solar downward flux
324 1536 : call pbuf_add_field('FSNS' , 'global',dtype_r8,(/pcols/), fsns_idx) ! Surface net shortwave flux
325 1536 : call pbuf_add_field('FSNT' , 'global',dtype_r8,(/pcols/), fsnt_idx) ! Top-of-model net shortwave flux
326 :
327 1536 : call pbuf_add_field('FLNS' , 'global',dtype_r8,(/pcols/), flns_idx) ! Surface net longwave flux
328 1536 : call pbuf_add_field('FLNT' , 'global',dtype_r8,(/pcols/), flnt_idx) ! Top-of-model net longwave flux
329 :
330 : ! If the namelist has been configured for preserving the spectral fluxes, then create
331 : ! physics buffer variables to store the results. This data is accessed by CARMA.
332 1536 : if (spectralflux) then
333 0 : call pbuf_add_field('SU' , 'global',dtype_r8,(/pcols,pverp,nswbands/), su_idx) ! shortwave upward flux (per band)
334 0 : call pbuf_add_field('SD' , 'global',dtype_r8,(/pcols,pverp,nswbands/), sd_idx) ! shortwave downward flux (per band)
335 0 : call pbuf_add_field('LU' , 'global',dtype_r8,(/pcols,pverp,nlwbands/), lu_idx) ! longwave upward flux (per band)
336 0 : call pbuf_add_field('LD' , 'global',dtype_r8,(/pcols,pverp,nlwbands/), ld_idx) ! longwave downward flux (per band)
337 : end if
338 :
339 : ! Register fields for offline radiation driver.
340 1536 : call rad_data_register()
341 :
342 1536 : end subroutine radiation_register
343 :
344 : !================================================================================================
345 :
346 226008 : function radiation_do(op, timestep)
347 :
348 : ! Return true if the specified operation is done this timestep.
349 :
350 : character(len=*), intent(in) :: op ! name of operation
351 : integer, intent(in), optional:: timestep
352 : logical :: radiation_do ! return value
353 :
354 : ! Local variables
355 : integer :: nstep ! current timestep number
356 : !-----------------------------------------------------------------------
357 :
358 226008 : if (present(timestep)) then
359 226008 : nstep = timestep
360 : else
361 0 : nstep = get_nstep()
362 : end if
363 :
364 164088 : select case (op)
365 : case ('sw') ! do a shortwave heating calc this timestep?
366 : radiation_do = nstep == 0 .or. iradsw == 1 &
367 : .or. (mod(nstep-1,iradsw) == 0 .and. nstep /= 1) &
368 164088 : .or. nstep <= irad_always
369 : case ('lw') ! do a longwave heating calc this timestep?
370 : radiation_do = nstep == 0 .or. iradlw == 1 &
371 : .or. (mod(nstep-1,iradlw) == 0 .and. nstep /= 1) &
372 61920 : .or. nstep <= irad_always
373 : case default
374 226008 : call endrun('radiation_do: unknown operation:'//op)
375 : end select
376 :
377 226008 : end function radiation_do
378 :
379 : !================================================================================================
380 :
381 61920 : real(r8) function radiation_nextsw_cday()
382 :
383 : ! If a SW radiation calculation will be done on the next time-step, then return
384 : ! the calendar day of that time-step. Otherwise return -1.0
385 :
386 : ! Local variables
387 : integer :: nstep ! timestep counter
388 : logical :: dosw ! true => do shosrtwave calc
389 : integer :: offset ! offset for calendar day calculation
390 : integer :: dtime ! integer timestep size
391 : real(r8):: caldayp1 ! calendar day of next time-step
392 :
393 : !-----------------------------------------------------------------------
394 :
395 61920 : radiation_nextsw_cday = -1._r8
396 61920 : dosw = .false.
397 61920 : nstep = get_nstep()
398 61920 : dtime = get_step_size()
399 61920 : offset = 0
400 : do while (.not. dosw)
401 102168 : nstep = nstep + 1
402 102168 : offset = offset + dtime
403 102168 : if (radiation_do('sw', nstep)) then
404 61920 : radiation_nextsw_cday = get_curr_calday(offset=offset)
405 : dosw = .true.
406 : end if
407 : end do
408 61920 : if(radiation_nextsw_cday == -1._r8) then
409 0 : call endrun('error in radiation_nextsw_cday')
410 : end if
411 :
412 : ! determine if next radiation time-step not equal to next time-step
413 61920 : if (get_nstep() >= 1) then
414 55728 : caldayp1 = get_curr_calday(offset=int(dtime))
415 55728 : if (caldayp1 /= radiation_nextsw_cday) radiation_nextsw_cday = -1._r8
416 : end if
417 :
418 61920 : end function radiation_nextsw_cday
419 :
420 : !================================================================================================
421 :
422 1536 : subroutine radiation_init(pbuf2d)
423 :
424 : ! Initialize the radiation and cloud optics.
425 : ! Add fields to the history buffer.
426 :
427 : ! arguments
428 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
429 :
430 : ! local variables
431 : character(len=128) :: errmsg
432 :
433 : ! names of gases that are available in the model
434 : ! -- needed for the kdist initialization routines
435 1536 : type(ty_gas_concs) :: available_gases
436 :
437 : integer :: i, icall
438 : integer :: nstep ! current timestep number
439 : logical :: history_amwg ! output the variables used by the AMWG diag package
440 : logical :: history_vdiag ! output the variables used by the AMWG variability diag package
441 : logical :: history_budget ! output tendencies and state variables for CAM4
442 : ! temperature, water vapor, cloud ice and cloud
443 : ! liquid budgets.
444 : integer :: history_budget_histfile_num ! history file number for budget fields
445 : integer :: ierr, istat
446 :
447 : integer :: dtime
448 :
449 : character(len=*), parameter :: sub = 'radiation_init'
450 : !-----------------------------------------------------------------------
451 :
452 : ! Number of layers in radiation calculation is capped by the number of
453 : ! pressure interfaces below 1 Pa. When the entire model atmosphere is
454 : ! below 1 Pa then an extra layer is added to the top of the model for
455 : ! the purpose of the radiation calculation.
456 :
457 145920 : nlay = count( pref_edge(:) > 1._r8 ) ! pascals (0.01 mbar)
458 :
459 1536 : if (nlay == pverp) then
460 : ! Top model interface is below 1 Pa. RRTMGP is active in all model layers plus
461 : ! 1 extra layer between model top and 1 Pa.
462 0 : ktopcam = 1
463 0 : ktoprad = 2
464 0 : nlaycam = pver
465 1536 : else if (nlay == (pverp-1)) then
466 : ! Special case nlay == (pverp-1) -- topmost interface outside bounds (CAM MT config), treat as if it is ok.
467 1536 : ktopcam = 1
468 1536 : ktoprad = 2
469 1536 : nlaycam = pver
470 1536 : nlay = nlay+1 ! reassign the value so later code understands to treat this case like nlay==pverp
471 1536 : write(iulog,*) 'RADIATION_INIT: Special case of 1 model interface at p < 1Pa. Top layer will be INCLUDED in radiation calculation.'
472 1536 : write(iulog,*) 'RADIATION_INIT: nlay = ',nlay, ' same as pverp: ',nlay==pverp
473 : else
474 : ! nlay < pverp. nlay layers are used in radiation calcs, and they are
475 : ! all CAM layers.
476 0 : ktopcam = pver - nlay + 1
477 0 : ktoprad = 1
478 0 : nlaycam = nlay
479 : end if
480 :
481 : ! Create lowercase version of the gaslist for RRTMGP. The ty_gas_concs objects
482 : ! work with CAM's uppercase names, but other objects that get input from the gas
483 : ! concs objects don't work.
484 13824 : do i = 1, nradgas
485 13824 : gaslist_lc(i) = to_lower(gaslist(i))
486 : end do
487 :
488 1536 : errmsg = available_gases%init(gaslist_lc)
489 1536 : call stop_on_err(errmsg, sub, 'available_gases%init')
490 :
491 : ! Read RRTMGP coefficients files and initialize kdist objects.
492 1536 : call coefs_init(coefs_sw_file, available_gases, kdist_sw)
493 1536 : call coefs_init(coefs_lw_file, available_gases, kdist_lw)
494 :
495 : ! Set the sw/lw band boundaries in radconstants. Also sets
496 : ! indicies of specific bands for diagnostic output and COSP input.
497 1536 : call set_wavenumber_bands(kdist_sw, kdist_lw)
498 1536 : call rad_solar_var_init()
499 :
500 : ! The spectral band boundaries need to be set before this init is called.
501 1536 : call rrtmgp_inputs_init(ktopcam, ktoprad)
502 :
503 : ! initialize output fields for offline driver
504 1536 : call rad_data_init(pbuf2d)
505 :
506 1536 : call cloud_rad_props_init()
507 :
508 1536 : cld_idx = pbuf_get_index('CLD')
509 1536 : cldfsnow_idx = pbuf_get_index('CLDFSNOW', errcode=ierr)
510 1536 : cldfgrau_idx = pbuf_get_index('CLDFGRAU', errcode=ierr)
511 :
512 1536 : if (is_first_step()) then
513 768 : call pbuf_set_field(pbuf2d, qrl_idx, 0._r8)
514 : end if
515 :
516 : ! Set the radiation timestep for cosz calculations if requested using
517 : ! the adjusted iradsw value from radiation
518 1536 : if (use_rad_dt_cosz) then
519 0 : dtime = get_step_size()
520 0 : dt_avg = iradsw*dtime
521 : end if
522 :
523 : ! Surface components to get radiation computed today
524 1536 : if (.not. is_first_restart_step()) then
525 768 : nextsw_cday = get_curr_calday()
526 : end if
527 :
528 : call phys_getopts(history_amwg_out = history_amwg, &
529 : history_vdiag_out = history_vdiag, &
530 : history_budget_out = history_budget, &
531 1536 : history_budget_histfile_num_out = history_budget_histfile_num)
532 :
533 : ! "irad_always" is number of time steps to execute radiation continuously from
534 : ! start of initial OR restart run
535 1536 : nstep = get_nstep()
536 1536 : if (irad_always > 0) then
537 0 : irad_always = irad_always + nstep
538 : end if
539 :
540 1536 : if (docosp) call cospsimulator_intr_init()
541 :
542 4608 : allocate(cosp_cnt(begchunk:endchunk), stat=istat)
543 1536 : call handle_allocate_error(istat, sub, 'cosp_cnt')
544 1536 : if (is_first_restart_step()) then
545 3864 : cosp_cnt(begchunk:endchunk) = cosp_cnt_init
546 : else
547 3864 : cosp_cnt(begchunk:endchunk) = 0
548 : end if
549 :
550 : ! Add fields to history buffer
551 :
552 : call addfld('TOT_CLD_VISTAU', (/ 'lev' /), 'A', '1', &
553 : 'Total gbx cloud extinction visible sw optical depth', &
554 3072 : sampling_seq='rad_lwsw', flag_xyfill=.true.)
555 : call addfld('TOT_ICLD_VISTAU', (/ 'lev' /), 'A', '1', &
556 : 'Total in-cloud extinction visible sw optical depth', &
557 3072 : sampling_seq='rad_lwsw', flag_xyfill=.true.)
558 : call addfld('LIQ_ICLD_VISTAU', (/ 'lev' /), 'A', '1', &
559 : 'Liquid in-cloud extinction visible sw optical depth', &
560 3072 : sampling_seq='rad_lwsw', flag_xyfill=.true.)
561 : call addfld('ICE_ICLD_VISTAU', (/ 'lev' /), 'A', '1', &
562 : 'Ice in-cloud extinction visible sw optical depth', &
563 3072 : sampling_seq='rad_lwsw', flag_xyfill=.true.)
564 :
565 1536 : if (cldfsnow_idx > 0) then
566 : call addfld('SNOW_ICLD_VISTAU', (/ 'lev' /), 'A', '1', &
567 : 'Snow in-cloud extinction visible sw optical depth', &
568 3072 : sampling_seq='rad_lwsw', flag_xyfill=.true.)
569 : end if
570 1536 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
571 : call addfld('GRAU_ICLD_VISTAU', (/ 'lev' /), 'A', '1', &
572 : 'Graupel in-cloud extinction visible sw optical depth', &
573 0 : sampling_seq='rad_lwsw', flag_xyfill=.true.)
574 : endif
575 :
576 : ! get list of active radiation calls
577 1536 : call rad_cnst_get_call_list(active_calls)
578 :
579 : ! Add shortwave radiation fields to history master field list.
580 :
581 18432 : do icall = 0, N_DIAG
582 :
583 18432 : if (active_calls(icall)) then
584 :
585 : call addfld('SOLIN'//diag(icall), horiz_only, 'A', 'W/m2', &
586 1536 : 'Solar insolation', sampling_seq='rad_lwsw')
587 : call addfld('QRS'//diag(icall), (/ 'lev' /), 'A', 'K/s', &
588 3072 : 'Solar heating rate', sampling_seq='rad_lwsw')
589 : call addfld('QRSC'//diag(icall), (/ 'lev' /), 'A', 'K/s', &
590 3072 : 'Clearsky solar heating rate', sampling_seq='rad_lwsw')
591 : call addfld('FSNT'//diag(icall), horiz_only, 'A', 'W/m2', &
592 1536 : 'Net solar flux at top of model', sampling_seq='rad_lwsw')
593 : call addfld('FSNTC'//diag(icall), horiz_only, 'A', 'W/m2', &
594 1536 : 'Clearsky net solar flux at top of model', sampling_seq='rad_lwsw')
595 : call addfld('FSNTOA'//diag(icall), horiz_only, 'A', 'W/m2', &
596 1536 : 'Net solar flux at top of atmosphere', sampling_seq='rad_lwsw')
597 : call addfld('FSNTOAC'//diag(icall), horiz_only, 'A', 'W/m2', &
598 1536 : 'Clearsky net solar flux at top of atmosphere', sampling_seq='rad_lwsw')
599 : call addfld('SWCF'//diag(icall), horiz_only, 'A', 'W/m2', &
600 1536 : 'Shortwave cloud forcing', sampling_seq='rad_lwsw')
601 : call addfld('FSUTOA'//diag(icall), horiz_only, 'A', 'W/m2', &
602 1536 : 'Upwelling solar flux at top of atmosphere', sampling_seq='rad_lwsw')
603 : call addfld('FSNIRTOA'//diag(icall), horiz_only, 'A', 'W/m2', &
604 1536 : 'Net near-infrared flux (Nimbus-7 WFOV) at top of atmosphere', sampling_seq='rad_lwsw')
605 : call addfld('FSNRTOAC'//diag(icall), horiz_only, 'A', 'W/m2', &
606 1536 : 'Clearsky net near-infrared flux (Nimbus-7 WFOV) at top of atmosphere', sampling_seq='rad_lwsw')
607 : call addfld('FSNRTOAS'//diag(icall), horiz_only, 'A', 'W/m2', &
608 1536 : 'Net near-infrared flux (>= 0.7 microns) at top of atmosphere', sampling_seq='rad_lwsw')
609 : call addfld('FSN200'//diag(icall), horiz_only, 'A', 'W/m2', &
610 1536 : 'Net shortwave flux at 200 mb', sampling_seq='rad_lwsw')
611 : call addfld('FSN200C'//diag(icall), horiz_only, 'A', 'W/m2', &
612 1536 : 'Clearsky net shortwave flux at 200 mb', sampling_seq='rad_lwsw')
613 : call addfld('FSNR'//diag(icall), horiz_only, 'A', 'W/m2', &
614 1536 : 'Net solar flux at tropopause', sampling_seq='rad_lwsw')
615 : call addfld('SOLL'//diag(icall), horiz_only, 'A', 'W/m2', &
616 1536 : 'Solar downward near infrared direct to surface', sampling_seq='rad_lwsw')
617 : call addfld('SOLS'//diag(icall), horiz_only, 'A', 'W/m2', &
618 1536 : 'Solar downward visible direct to surface', sampling_seq='rad_lwsw')
619 : call addfld('SOLLD'//diag(icall), horiz_only, 'A', 'W/m2', &
620 1536 : 'Solar downward near infrared diffuse to surface', sampling_seq='rad_lwsw')
621 : call addfld('SOLSD'//diag(icall), horiz_only, 'A', 'W/m2', &
622 1536 : 'Solar downward visible diffuse to surface', sampling_seq='rad_lwsw')
623 : call addfld('FSNS'//diag(icall), horiz_only, 'A', 'W/m2', &
624 1536 : 'Net solar flux at surface', sampling_seq='rad_lwsw')
625 : call addfld('FSNSC'//diag(icall), horiz_only, 'A', 'W/m2', &
626 1536 : 'Clearsky net solar flux at surface', sampling_seq='rad_lwsw')
627 : call addfld('FSDS'//diag(icall), horiz_only, 'A', 'W/m2', &
628 1536 : 'Downwelling solar flux at surface', sampling_seq='rad_lwsw')
629 : call addfld('FSDSC'//diag(icall), horiz_only, 'A', 'W/m2', &
630 1536 : 'Clearsky downwelling solar flux at surface', sampling_seq='rad_lwsw')
631 :
632 : ! Fluxes on CAM grid
633 : call addfld('FUS'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', &
634 3072 : 'Shortwave upward flux', sampling_seq='rad_lwsw')
635 : call addfld('FDS'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', &
636 3072 : 'Shortwave downward flux', sampling_seq='rad_lwsw')
637 : call addfld('FUSC'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', &
638 3072 : 'Shortwave clear-sky upward flux', sampling_seq='rad_lwsw')
639 : call addfld('FDSC'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', &
640 3072 : 'Shortwave clear-sky downward flux', sampling_seq='rad_lwsw')
641 :
642 1536 : if (history_amwg) then
643 1536 : call add_default('SOLIN'//diag(icall), 1, ' ')
644 1536 : call add_default('QRS'//diag(icall), 1, ' ')
645 1536 : call add_default('FSNT'//diag(icall), 1, ' ')
646 1536 : call add_default('FSNTC'//diag(icall), 1, ' ')
647 1536 : call add_default('FSNTOA'//diag(icall), 1, ' ')
648 1536 : call add_default('FSNTOAC'//diag(icall), 1, ' ')
649 1536 : call add_default('SWCF'//diag(icall), 1, ' ')
650 1536 : call add_default('FSNS'//diag(icall), 1, ' ')
651 1536 : call add_default('FSNSC'//diag(icall), 1, ' ')
652 1536 : call add_default('FSUTOA'//diag(icall), 1, ' ')
653 1536 : call add_default('FSDSC'//diag(icall), 1, ' ')
654 1536 : call add_default('FSDS'//diag(icall), 1, ' ')
655 : endif
656 :
657 : end if
658 : end do
659 :
660 1536 : if (scm_crm_mode) then
661 0 : call add_default('FUS ', 1, ' ')
662 0 : call add_default('FUSC ', 1, ' ')
663 0 : call add_default('FDS ', 1, ' ')
664 0 : call add_default('FDSC ', 1, ' ')
665 : endif
666 :
667 : ! Add longwave radiation fields to history master field list.
668 :
669 18432 : do icall = 0, N_DIAG
670 :
671 18432 : if (active_calls(icall)) then
672 : call addfld('QRL'//diag(icall), (/ 'lev' /), 'A', 'K/s', &
673 3072 : 'Longwave heating rate', sampling_seq='rad_lwsw')
674 : call addfld('QRLC'//diag(icall), (/ 'lev' /), 'A', 'K/s', &
675 3072 : 'Clearsky longwave heating rate', sampling_seq='rad_lwsw')
676 : call addfld('FLNT'//diag(icall), horiz_only, 'A', 'W/m2', &
677 1536 : 'Net longwave flux at top of model', sampling_seq='rad_lwsw')
678 : call addfld('FLNTC'//diag(icall), horiz_only, 'A', 'W/m2', &
679 1536 : 'Clearsky net longwave flux at top of model', sampling_seq='rad_lwsw')
680 : call addfld('FLUT'//diag(icall), horiz_only, 'A', 'W/m2', &
681 1536 : 'Upwelling longwave flux at top of model', sampling_seq='rad_lwsw')
682 : call addfld('FLUTC'//diag(icall), horiz_only, 'A', 'W/m2', &
683 1536 : 'Clearsky upwelling longwave flux at top of model', sampling_seq='rad_lwsw')
684 : call addfld('LWCF'//diag(icall), horiz_only, 'A', 'W/m2', &
685 1536 : 'Longwave cloud forcing', sampling_seq='rad_lwsw')
686 : call addfld('FLN200'//diag(icall), horiz_only, 'A', 'W/m2', &
687 1536 : 'Net longwave flux at 200 mb', sampling_seq='rad_lwsw')
688 : call addfld('FLN200C'//diag(icall), horiz_only, 'A', 'W/m2', &
689 1536 : 'Clearsky net longwave flux at 200 mb', sampling_seq='rad_lwsw')
690 : call addfld('FLNR'//diag(icall), horiz_only, 'A', 'W/m2', &
691 1536 : 'Net longwave flux at tropopause', sampling_seq='rad_lwsw')
692 : call addfld('FLNS'//diag(icall), horiz_only, 'A', 'W/m2', &
693 1536 : 'Net longwave flux at surface', sampling_seq='rad_lwsw')
694 : call addfld('FLNSC'//diag(icall), horiz_only, 'A', 'W/m2', &
695 1536 : 'Clearsky net longwave flux at surface', sampling_seq='rad_lwsw')
696 : call addfld('FLDS'//diag(icall), horiz_only, 'A', 'W/m2', &
697 1536 : 'Downwelling longwave flux at surface', sampling_seq='rad_lwsw')
698 : call addfld('FLDSC'//diag(icall), horiz_only, 'A', 'W/m2', &
699 1536 : 'Clearsky Downwelling longwave flux at surface', sampling_seq='rad_lwsw')
700 :
701 : ! Fluxes on CAM grid
702 : call addfld('FUL'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', &
703 3072 : 'Longwave upward flux', sampling_seq='rad_lwsw')
704 : call addfld('FDL'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', &
705 3072 : 'Longwave downward flux', sampling_seq='rad_lwsw')
706 : call addfld('FULC'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', &
707 3072 : 'Longwave clear-sky upward flux', sampling_seq='rad_lwsw')
708 : call addfld('FDLC'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', &
709 3072 : 'Longwave clear-sky downward flux', sampling_seq='rad_lwsw')
710 :
711 1536 : if (history_amwg) then
712 1536 : call add_default('QRL'//diag(icall), 1, ' ')
713 1536 : call add_default('FLNT'//diag(icall), 1, ' ')
714 1536 : call add_default('FLNTC'//diag(icall), 1, ' ')
715 1536 : call add_default('FLUT'//diag(icall), 1, ' ')
716 1536 : call add_default('FLUTC'//diag(icall), 1, ' ')
717 1536 : call add_default('LWCF'//diag(icall), 1, ' ')
718 1536 : call add_default('FLNS'//diag(icall), 1, ' ')
719 1536 : call add_default('FLNSC'//diag(icall), 1, ' ')
720 1536 : call add_default('FLDS'//diag(icall), 1, ' ')
721 : end if
722 :
723 : end if
724 : end do
725 :
726 3072 : call addfld('EMIS', (/ 'lev' /), 'A', '1', 'Cloud longwave emissivity')
727 :
728 1536 : if (scm_crm_mode) then
729 0 : call add_default ('FUL ', 1, ' ')
730 0 : call add_default ('FULC ', 1, ' ')
731 0 : call add_default ('FDL ', 1, ' ')
732 0 : call add_default ('FDLC ', 1, ' ')
733 : endif
734 :
735 : ! Heating rate needed for d(theta)/dt computation
736 3072 : call addfld ('HR',(/ 'lev' /), 'A','K/s','Heating rate needed for d(theta)/dt computation')
737 :
738 1536 : if ( history_budget .and. history_budget_histfile_num > 1 ) then
739 0 : call add_default ('QRL ', history_budget_histfile_num, ' ')
740 0 : call add_default ('QRS ', history_budget_histfile_num, ' ')
741 : end if
742 :
743 1536 : if (history_vdiag) then
744 0 : call add_default('FLUT', 2, ' ')
745 0 : call add_default('FLUT', 3, ' ')
746 : end if
747 :
748 3072 : end subroutine radiation_init
749 :
750 : !===============================================================================
751 :
752 1536 : subroutine radiation_define_restart(file)
753 :
754 : ! define variables to be written to restart file
755 :
756 : ! arguments
757 : type(file_desc_t), intent(inout) :: file
758 :
759 : ! local variables
760 : integer :: ierr
761 : !----------------------------------------------------------------------------
762 :
763 1536 : call pio_seterrorhandling(file, PIO_BCAST_ERROR)
764 :
765 1536 : ierr = pio_def_var(file, 'nextsw_cday', pio_double, nextsw_cday_desc)
766 1536 : ierr = pio_put_att(file, nextsw_cday_desc, 'long_name', 'future radiation calday for surface models')
767 1536 : if (docosp) then
768 0 : ierr = pio_def_var(File, 'cosp_cnt_init', pio_int, cospcnt_desc)
769 : end if
770 :
771 1536 : end subroutine radiation_define_restart
772 :
773 : !===============================================================================
774 :
775 1536 : subroutine radiation_write_restart(file)
776 :
777 : ! write variables to restart file
778 :
779 : ! arguments
780 : type(file_desc_t), intent(inout) :: file
781 :
782 : ! local variables
783 : integer :: ierr
784 : !----------------------------------------------------------------------------
785 3072 : ierr = pio_put_var(File, nextsw_cday_desc, (/ nextsw_cday /))
786 1536 : if (docosp) then
787 0 : ierr = pio_put_var(File, cospcnt_desc, (/cosp_cnt(begchunk)/))
788 : end if
789 :
790 1536 : end subroutine radiation_write_restart
791 :
792 : !===============================================================================
793 :
794 768 : subroutine radiation_read_restart(file)
795 :
796 : ! read variables from restart file
797 :
798 : ! arguments
799 : type(file_desc_t), intent(inout) :: file
800 :
801 : ! local variables
802 : integer :: ierr
803 : type(var_desc_t) :: vardesc
804 : integer :: err_handling
805 :
806 : !----------------------------------------------------------------------------
807 768 : if (docosp) then
808 0 : call pio_seterrorhandling(File, PIO_BCAST_ERROR, err_handling)
809 0 : ierr = pio_inq_varid(File, 'cosp_cnt_init', vardesc)
810 0 : call pio_seterrorhandling(File, err_handling)
811 0 : if (ierr /= PIO_NOERR) then
812 0 : cosp_cnt_init = 0
813 : else
814 0 : ierr = pio_get_var(File, vardesc, cosp_cnt_init)
815 : end if
816 : end if
817 :
818 768 : ierr = pio_inq_varid(file, 'nextsw_cday', vardesc)
819 768 : ierr = pio_get_var(file, vardesc, nextsw_cday)
820 :
821 :
822 768 : end subroutine radiation_read_restart
823 :
824 : !===============================================================================
825 :
826 2600640 : subroutine radiation_tend( &
827 : state, ptend, pbuf, cam_out, cam_in, net_flx, rd_out)
828 :
829 : !-----------------------------------------------------------------------
830 : !
831 : ! CAM driver for radiation computation.
832 : !
833 : !-----------------------------------------------------------------------
834 :
835 : ! Location/Orbital Parameters for cosine zenith angle
836 : use phys_grid, only: get_rlat_all_p, get_rlon_all_p
837 : use cam_control_mod, only: eccen, mvelpp, lambm0, obliqr
838 : use shr_orb_mod, only: shr_orb_decl, shr_orb_cosz
839 :
840 : use rrtmgp_inputs, only: rrtmgp_set_state, rrtmgp_set_gases_lw, rrtmgp_set_cloud_lw, &
841 : rrtmgp_set_aer_lw, rrtmgp_set_gases_sw, rrtmgp_set_cloud_sw, &
842 : rrtmgp_set_aer_sw
843 :
844 : ! RRTMGP drivers for flux calculations.
845 : use mo_rte_lw, only: rte_lw
846 : use mo_rte_sw, only: rte_sw
847 :
848 : use radheat, only: radheat_tend
849 :
850 : use radiation_data, only: rad_data_write
851 :
852 : use interpolate_data, only: vertinterp
853 : use tropopause, only: tropopause_find_cam, TROP_ALG_HYBSTOB, TROP_ALG_CLIMATE
854 : use cospsimulator_intr, only: docosp, cospsimulator_intr_run, cosp_nradsteps
855 :
856 :
857 : ! Arguments
858 : type(physics_state), intent(in), target :: state
859 : type(physics_ptend), intent(out) :: ptend
860 : type(physics_buffer_desc), pointer :: pbuf(:)
861 : type(cam_out_t), intent(inout) :: cam_out
862 : type(cam_in_t), intent(in) :: cam_in
863 : real(r8), intent(out) :: net_flx(pcols)
864 :
865 : type(rad_out_t), target, optional, intent(out) :: rd_out
866 :
867 :
868 : ! Local variables
869 : type(rad_out_t), pointer :: rd ! allow rd_out to be optional by allocating a local object
870 : ! if the argument is not present
871 : logical :: write_output
872 :
873 : integer :: i, k, istat
874 : integer :: lchnk, ncol
875 : logical :: dosw, dolw
876 : integer :: icall ! loop index for climate/diagnostic radiation calls
877 :
878 : real(r8) :: calday ! current calendar day
879 : real(r8) :: delta ! Solar declination angle in radians
880 : real(r8) :: eccf ! Earth orbit eccentricity factor
881 : real(r8) :: clat(pcols) ! current latitudes(radians)
882 : real(r8) :: clon(pcols) ! current longitudes(radians)
883 : real(r8) :: coszrs(pcols) ! Cosine solar zenith angle
884 :
885 : ! Gathered indices of day and night columns
886 : ! chunk_column_index = IdxDay(daylight_column_index)
887 : integer :: Nday ! Number of daylight columns
888 : integer :: Nnite ! Number of night columns
889 : integer :: IdxDay(pcols) ! chunk indices of daylight columns
890 : integer :: IdxNite(pcols) ! chunk indices of night columns
891 :
892 : integer :: itim_old
893 :
894 61920 : real(r8), pointer :: cld(:,:) ! cloud fraction
895 61920 : real(r8), pointer :: cldfsnow(:,:) ! cloud fraction of just "snow clouds"
896 61920 : real(r8), pointer :: cldfgrau(:,:) ! cloud fraction of just "graupel clouds"
897 : real(r8) :: cldfprime(pcols,pver) ! combined cloud fraction
898 61920 : real(r8), pointer :: qrs(:,:) ! shortwave radiative heating rate
899 61920 : real(r8), pointer :: qrl(:,:) ! longwave radiative heating rate
900 61920 : real(r8), pointer :: fsds(:) ! Surface solar down flux
901 61920 : real(r8), pointer :: fsns(:) ! Surface solar absorbed flux
902 61920 : real(r8), pointer :: fsnt(:) ! Net column abs solar flux at model top
903 61920 : real(r8), pointer :: flns(:) ! Srf longwave cooling (up-down) flux
904 61920 : real(r8), pointer :: flnt(:) ! Net outgoing lw flux at model top
905 :
906 : real(r8), pointer, dimension(:,:,:) :: su => NULL() ! shortwave spectral flux up
907 : real(r8), pointer, dimension(:,:,:) :: sd => NULL() ! shortwave spectral flux down
908 : real(r8), pointer, dimension(:,:,:) :: lu => NULL() ! longwave spectral flux up
909 : real(r8), pointer, dimension(:,:,:) :: ld => NULL() ! longwave spectral flux down
910 :
911 : ! tropopause diagnostic
912 : integer :: troplev(pcols)
913 : real(r8) :: p_trop(pcols)
914 :
915 : ! state data passed to radiation calc
916 61920 : real(r8), allocatable :: t_sfc(:)
917 61920 : real(r8), allocatable :: emis_sfc(:,:)
918 61920 : real(r8), allocatable :: t_rad(:,:)
919 61920 : real(r8), allocatable :: pmid_rad(:,:)
920 61920 : real(r8), allocatable :: pint_rad(:,:)
921 61920 : real(r8), allocatable :: t_day(:,:)
922 61920 : real(r8), allocatable :: pmid_day(:,:)
923 61920 : real(r8), allocatable :: pint_day(:,:)
924 61920 : real(r8), allocatable :: coszrs_day(:)
925 61920 : real(r8), allocatable :: alb_dir(:,:)
926 61920 : real(r8), allocatable :: alb_dif(:,:)
927 :
928 : ! in-cloud optical depths for COSP
929 : real(r8) :: cld_tau_cloudsim(pcols,pver) ! liq + ice
930 : real(r8) :: snow_tau_cloudsim(pcols,pver) ! snow
931 : real(r8) :: grau_tau_cloudsim(pcols,pver) ! graupel
932 : real(r8) :: cld_lw_abs_cloudsim(pcols,pver) ! liq + ice
933 : real(r8) :: snow_lw_abs_cloudsim(pcols,pver)! snow
934 : real(r8) :: grau_lw_abs_cloudsim(pcols,pver)! graupel
935 :
936 : ! Set vertical indexing in RRTMGP to be the same as CAM (top to bottom).
937 : logical, parameter :: top_at_1 = .true.
938 :
939 : ! TOA solar flux on RRTMGP g-points
940 61920 : real(r8), allocatable :: toa_flux(:,:)
941 : ! Scale factors based on spectral distribution from input irradiance dataset
942 61920 : real(r8), allocatable :: sfac(:,:)
943 :
944 : ! Planck sources for LW.
945 61920 : type(ty_source_func_lw) :: sources_lw
946 :
947 : ! Gas volume mixing ratios. Use separate objects for LW and SW because SW only does
948 : ! calculations for daylight columns.
949 : ! These objects have a final method which deallocates the internal memory when they
950 : ! go out of scope (i.e., when radiation_tend returns), so no need for explicit deallocation.
951 61920 : type(ty_gas_concs) :: gas_concs_lw
952 61920 : type(ty_gas_concs) :: gas_concs_sw
953 :
954 : ! Atmosphere optics. This object is initialized with gas optics, then is incremented
955 : ! by the aerosol optics for the clear-sky radiative flux calculations, and then
956 : ! incremented again by the cloud optics for the all-sky radiative flux calculations.
957 61920 : type(ty_optical_props_1scl) :: atm_optics_lw
958 61920 : type(ty_optical_props_2str) :: atm_optics_sw
959 :
960 : ! Cloud optical properties objects (McICA sampling of cloud optical properties).
961 61920 : type(ty_optical_props_1scl) :: cloud_lw
962 61920 : type(ty_optical_props_2str) :: cloud_sw
963 :
964 : ! Aerosol optical properties objects.
965 61920 : type(ty_optical_props_1scl) :: aer_lw
966 61920 : type(ty_optical_props_2str) :: aer_sw
967 :
968 : ! Flux objects contain all fluxes computed by RRTMGP.
969 : ! SW allsky fluxes always include spectrally resolved fluxes needed for surface models.
970 : type(ty_fluxes_byband) :: fsw
971 : ! LW allsky fluxes only need spectrally resolved fluxes when spectralflux=.true.
972 : type(ty_fluxes_byband) :: flw
973 : ! Only broadband fluxes needed for clear sky (diagnostics).
974 : type(ty_fluxes_broadband) :: fswc, flwc
975 :
976 : ! Arrays for output diagnostics on CAM grid.
977 : real(r8) :: fns(pcols,pverp) ! net shortwave flux
978 : real(r8) :: fcns(pcols,pverp) ! net clear-sky shortwave flux
979 : real(r8) :: fnl(pcols,pverp) ! net longwave flux
980 : real(r8) :: fcnl(pcols,pverp) ! net clear-sky longwave flux
981 :
982 : ! for COSP
983 : real(r8) :: emis(pcols,pver) ! Cloud longwave emissivity
984 : real(r8) :: gb_snow_tau(pcols,pver) ! grid-box mean snow_tau
985 : real(r8) :: gb_snow_lw(pcols,pver) ! grid-box mean LW snow optical depth
986 :
987 : real(r8) :: ftem(pcols,pver) ! Temporary workspace for outfld variables
988 :
989 : character(len=128) :: errmsg
990 : character(len=*), parameter :: sub = 'radiation_tend'
991 : !--------------------------------------------------------------------------------------
992 :
993 61920 : lchnk = state%lchnk
994 61920 : ncol = state%ncol
995 :
996 61920 : if (present(rd_out)) then
997 0 : rd => rd_out
998 0 : write_output = .false.
999 : else
1000 61920 : allocate(rd, stat=istat)
1001 61920 : call handle_allocate_error(istat, sub, 'rd')
1002 61920 : write_output = .true.
1003 : end if
1004 :
1005 61920 : dosw = radiation_do('sw', get_nstep()) ! do shortwave radiation calc this timestep?
1006 61920 : dolw = radiation_do('lw', get_nstep()) ! do longwave radiation calc this timestep?
1007 :
1008 : ! Cosine solar zenith angle for current time step
1009 61920 : calday = get_curr_calday()
1010 61920 : call get_rlat_all_p(lchnk, ncol, clat)
1011 61920 : call get_rlon_all_p(lchnk, ncol, clon)
1012 :
1013 : call shr_orb_decl(calday, eccen, mvelpp, lambm0, obliqr, &
1014 61920 : delta, eccf)
1015 :
1016 61920 : if (use_rad_uniform_angle) then
1017 0 : do i = 1, ncol
1018 0 : coszrs(i) = shr_orb_cosz(calday, clat(i), clon(i), delta, dt_avg, &
1019 0 : uniform_angle=rad_uniform_angle)
1020 : end do
1021 : else
1022 1033920 : do i = 1, ncol
1023 : ! if dt_avg /= 0, it triggers using avg coszrs
1024 1033920 : coszrs(i) = shr_orb_cosz(calday, clat(i), clon(i), delta, dt_avg)
1025 : end do
1026 : end if
1027 :
1028 : ! Gather night/day column indices.
1029 61920 : Nday = 0
1030 61920 : Nnite = 0
1031 1033920 : do i = 1, ncol
1032 1033920 : if ( coszrs(i) > 0.0_r8 ) then
1033 486000 : Nday = Nday + 1
1034 486000 : IdxDay(Nday) = i
1035 : else
1036 486000 : Nnite = Nnite + 1
1037 486000 : IdxNite(Nnite) = i
1038 : end if
1039 : end do
1040 :
1041 : ! Associate pointers to physics buffer fields
1042 61920 : itim_old = pbuf_old_tim_idx()
1043 61920 : nullify(cldfsnow)
1044 61920 : if (cldfsnow_idx > 0) then
1045 247680 : call pbuf_get_field(pbuf, cldfsnow_idx, cldfsnow, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
1046 : end if
1047 61920 : nullify(cldfgrau)
1048 61920 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
1049 0 : call pbuf_get_field(pbuf, cldfgrau_idx, cldfgrau, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
1050 : endif
1051 :
1052 247680 : call pbuf_get_field(pbuf, cld_idx, cld, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
1053 :
1054 61920 : call pbuf_get_field(pbuf, qrs_idx, qrs)
1055 61920 : call pbuf_get_field(pbuf, qrl_idx, qrl)
1056 :
1057 61920 : call pbuf_get_field(pbuf, fsns_idx, fsns)
1058 61920 : call pbuf_get_field(pbuf, fsnt_idx, fsnt)
1059 61920 : call pbuf_get_field(pbuf, fsds_idx, fsds)
1060 61920 : call pbuf_get_field(pbuf, flns_idx, flns)
1061 61920 : call pbuf_get_field(pbuf, flnt_idx, flnt)
1062 :
1063 61920 : if (spectralflux) then
1064 0 : call pbuf_get_field(pbuf, su_idx, su)
1065 0 : call pbuf_get_field(pbuf, sd_idx, sd)
1066 0 : call pbuf_get_field(pbuf, lu_idx, lu)
1067 0 : call pbuf_get_field(pbuf, ld_idx, ld)
1068 : end if
1069 :
1070 : ! Allocate the flux arrays and init to zero.
1071 61920 : call initialize_rrtmgp_fluxes(nday, nlay+1, nswbands, fsw, do_direct=.true.)
1072 61920 : call initialize_rrtmgp_fluxes(nday, nlay+1, nswbands, fswc, do_direct=.true.)
1073 61920 : call initialize_rrtmgp_fluxes(ncol, nlay+1, nlwbands, flw)
1074 61920 : call initialize_rrtmgp_fluxes(ncol, nlay+1, nlwbands, flwc)
1075 :
1076 : ! For CRM, make cloud equal to input observations:
1077 61920 : if (scm_crm_mode .and. have_cld) then
1078 0 : do k = 1, pver
1079 0 : cld(:ncol,k)= cldobs(k)
1080 : end do
1081 : end if
1082 :
1083 : ! Find tropopause height if needed for diagnostic output
1084 61920 : if (hist_fld_active('FSNR') .or. hist_fld_active('FLNR')) then
1085 : !REMOVECAM - no longer need this when CAM is retired and pcols no longer exists
1086 0 : troplev(:) = 0
1087 0 : p_trop(:) = 0._r8
1088 : !REMOVECAM_END
1089 : call tropopause_find_cam(state, troplev, tropP=p_trop, primary=TROP_ALG_HYBSTOB, &
1090 0 : backup=TROP_ALG_CLIMATE)
1091 : end if
1092 :
1093 : ! Get time of next radiation calculation - albedos will need to be
1094 : ! calculated by each surface model at this time
1095 61920 : nextsw_cday = radiation_nextsw_cday()
1096 :
1097 61920 : if (dosw .or. dolw) then
1098 :
1099 : allocate( &
1100 : t_sfc(ncol), emis_sfc(nlwbands,ncol), toa_flux(nday,nswgpts), &
1101 : sfac(nday,nswgpts), &
1102 : t_rad(ncol,nlay), pmid_rad(ncol,nlay), pint_rad(ncol,nlay+1), &
1103 : t_day(nday,nlay), pmid_day(nday,nlay), pint_day(nday,nlay+1), &
1104 : coszrs_day(nday), alb_dir(nswbands,nday), alb_dif(nswbands,nday), &
1105 841288 : stat=istat)
1106 30960 : call handle_allocate_error(istat, sub, 't_sfc,..,alb_dif')
1107 :
1108 : ! Prepares state variables, daylit columns, albedos for RRTMGP
1109 : call rrtmgp_set_state( &
1110 : state, cam_in, ncol, nlay, nday, &
1111 : idxday, coszrs, kdist_sw, t_sfc, emis_sfc, &
1112 : t_rad, pmid_rad, pint_rad, t_day, pmid_day, &
1113 30960 : pint_day, coszrs_day, alb_dir, alb_dif)
1114 :
1115 : ! Output the mass per layer, and total column burdens for gas and aerosol
1116 : ! constituents in the climate list.
1117 30960 : call rad_cnst_out(0, state, pbuf)
1118 :
1119 : ! Modified cloud fraction accounts for radiatively active snow and/or graupel
1120 30960 : call modified_cloud_fraction(ncol, cld, cldfsnow, cldfgrau, cldfprime)
1121 :
1122 : !========================!
1123 : ! SHORTWAVE calculations !
1124 : !========================!
1125 :
1126 30960 : if (dosw) then
1127 :
1128 : ! Set cloud optical properties in cloud_sw object.
1129 : call rrtmgp_set_cloud_sw( &
1130 : state, pbuf, nlay, nday, idxday, &
1131 : nnite, idxnite, pmid_day, cld, cldfsnow, &
1132 : cldfgrau, cldfprime, graupel_in_rad, kdist_sw, cloud_sw, &
1133 : rd%tot_cld_vistau, rd%tot_icld_vistau, rd%liq_icld_vistau, &
1134 : rd%ice_icld_vistau, rd%snow_icld_vistau, rd%grau_icld_vistau, &
1135 30960 : cld_tau_cloudsim, snow_tau_cloudsim, grau_tau_cloudsim )
1136 :
1137 30960 : if (write_output) then
1138 30960 : call radiation_output_cld(lchnk, rd)
1139 : end if
1140 :
1141 : ! If no daylight columns, can't create empty RRTMGP objects
1142 30960 : if (nday > 0) then
1143 :
1144 : ! Initialize object for gas concentrations.
1145 16151 : errmsg = gas_concs_sw%init(gaslist_lc)
1146 16151 : call stop_on_err(errmsg, sub, 'gas_concs_sw%init')
1147 :
1148 : ! Initialize object for combined gas + aerosol + cloud optics.
1149 : ! Allocates arrays for properties represented on g-points.
1150 16151 : errmsg = atm_optics_sw%alloc_2str(nday, nlay, kdist_sw)
1151 16151 : call stop_on_err(errmsg, sub, 'atm_optics_sw%alloc_2str')
1152 :
1153 : ! Initialize object for SW aerosol optics. Allocates arrays
1154 : ! for properties represented by band.
1155 16151 : errmsg = aer_sw%alloc_2str(nday, nlay, kdist_sw%get_band_lims_wavenumber())
1156 16151 : call stop_on_err(errmsg, sub, 'aer_sw%alloc_2str')
1157 :
1158 : end if
1159 :
1160 : ! The climate (icall==0) calculation must occur last.
1161 371520 : do icall = N_DIAG, 0, -1
1162 371520 : if (active_calls(icall)) then
1163 :
1164 30960 : if (nday > 0) then
1165 :
1166 : ! Set gas volume mixing ratios for this call in gas_concs_sw.
1167 : call rrtmgp_set_gases_sw( &
1168 : icall, state, pbuf, nlay, nday, &
1169 16151 : idxday, gas_concs_sw)
1170 :
1171 : ! Compute the gas optics (stored in atm_optics_sw).
1172 : ! toa_flux is the reference solar source from RRTMGP data.
1173 : !$acc data copyin(kdist_sw,pmid_day,pint_day,t_day,gas_concs_sw) &
1174 : !$acc copy(atm_optics_sw) &
1175 : !$acc copyout(toa_flux)
1176 : errmsg = kdist_sw%gas_optics( &
1177 : pmid_day, pint_day, t_day, gas_concs_sw, atm_optics_sw, &
1178 16151 : toa_flux)
1179 : !$acc end data
1180 16151 : call stop_on_err(errmsg, sub, 'kdist_sw%gas_optics')
1181 :
1182 : ! Scale the solar source
1183 16151 : call get_variability(toa_flux, sfac)
1184 29057214 : toa_flux = toa_flux * sfac * eccf
1185 :
1186 : end if
1187 :
1188 : ! Set SW aerosol optical properties in the aer_sw object.
1189 : ! This call made even when no daylight columns because it does some
1190 : ! diagnostic aerosol output.
1191 : call rrtmgp_set_aer_sw( &
1192 30960 : icall, state, pbuf, nday, idxday, nnite, idxnite, aer_sw)
1193 :
1194 30960 : if (nday > 0) then
1195 :
1196 : ! Increment the gas optics (in atm_optics_sw) by the aerosol optics in aer_sw.
1197 : !$acc data copyin(coszrs_day, toa_flux, alb_dir, alb_dif, &
1198 : !$acc atm_optics_sw, atm_optics_sw%tau, &
1199 : !$acc atm_optics_sw%ssa, atm_optics_sw%g, &
1200 : !$acc aer_sw, aer_sw%tau, &
1201 : !$acc aer_sw%ssa, aer_sw%g, &
1202 : !$acc cloud_sw, cloud_sw%tau, &
1203 : !$acc cloud_sw%ssa, cloud_sw%g) &
1204 : !$acc copy(fswc, fswc%flux_net,fswc%flux_up,fswc%flux_dn, &
1205 : !$acc fsw, fsw%flux_net, fsw%flux_up, fsw%flux_dn)
1206 16151 : errmsg = aer_sw%increment(atm_optics_sw)
1207 16151 : call stop_on_err(errmsg, sub, 'aer_sw%increment')
1208 :
1209 : ! Compute clear-sky fluxes.
1210 : errmsg = rte_sw(&
1211 : atm_optics_sw, top_at_1, coszrs_day, toa_flux, &
1212 16151 : alb_dir, alb_dif, fswc)
1213 16151 : call stop_on_err(errmsg, sub, 'clear-sky rte_sw')
1214 :
1215 : ! Increment the aerosol+gas optics (in atm_optics_sw) by the cloud optics in cloud_sw.
1216 16151 : errmsg = cloud_sw%increment(atm_optics_sw)
1217 16151 : call stop_on_err(errmsg, sub, 'cloud_sw%increment')
1218 :
1219 : ! Compute all-sky fluxes.
1220 : errmsg = rte_sw(&
1221 : atm_optics_sw, top_at_1, coszrs_day, toa_flux, &
1222 16151 : alb_dir, alb_dif, fsw)
1223 16151 : call stop_on_err(errmsg, sub, 'all-sky rte_sw')
1224 : !$acc end data
1225 : end if
1226 :
1227 : ! Transform RRTMGP outputs to CAM outputs and compute heating rates.
1228 30960 : call set_sw_diags()
1229 :
1230 30960 : if (write_output) then
1231 30960 : call radiation_output_sw(lchnk, ncol, icall, rd, pbuf, cam_out)
1232 : end if
1233 :
1234 : end if ! (active_calls(icall))
1235 : end do ! loop over diagnostic calcs (icall)
1236 :
1237 : else
1238 : ! SW calc not done. pbuf carries Q*dp across timesteps.
1239 : ! Convert to Q before calling radheat_tend.
1240 0 : qrs(:ncol,:) = qrs(:ncol,:) / state%pdel(:ncol,:)
1241 : end if ! if (dosw)
1242 :
1243 : !=======================!
1244 : ! LONGWAVE calculations !
1245 : !=======================!
1246 :
1247 30960 : if (dolw) then
1248 :
1249 : ! Initialize object for Planck sources.
1250 30960 : errmsg = sources_lw%alloc(ncol, nlay, kdist_lw)
1251 30960 : call stop_on_err(errmsg, sub, 'sources_lw%alloc')
1252 :
1253 : ! Set cloud optical properties in cloud_lw object.
1254 : call rrtmgp_set_cloud_lw( &
1255 : state, pbuf, ncol, nlay, nlaycam, &
1256 : cld, cldfsnow, cldfgrau, cldfprime, graupel_in_rad, &
1257 30960 : kdist_lw, cloud_lw, cld_lw_abs_cloudsim, snow_lw_abs_cloudsim, grau_lw_abs_cloudsim)
1258 :
1259 : ! Initialize object for gas concentrations
1260 30960 : errmsg = gas_concs_lw%init(gaslist_lc)
1261 30960 : call stop_on_err(errmsg, sub, 'gas_concs_lw%init')
1262 :
1263 : ! Initialize object for combined gas + aerosol + cloud optics.
1264 30960 : errmsg = atm_optics_lw%alloc_1scl(ncol, nlay, kdist_lw)
1265 30960 : call stop_on_err(errmsg, sub, 'atm_optics_lw%alloc_1scl')
1266 :
1267 : ! Initialize object for LW aerosol optics.
1268 30960 : errmsg = aer_lw%alloc_1scl(ncol, nlay, kdist_lw%get_band_lims_wavenumber())
1269 30960 : call stop_on_err(errmsg, sub, 'aer_lw%alloc_1scl')
1270 :
1271 : ! The climate (icall==0) calculation must occur last.
1272 371520 : do icall = N_DIAG, 0, -1
1273 :
1274 371520 : if (active_calls(icall)) then
1275 :
1276 : ! Set gas volume mixing ratios for this call in gas_concs_lw.
1277 30960 : call rrtmgp_set_gases_lw(icall, state, pbuf, nlay, gas_concs_lw)
1278 :
1279 : ! Compute the gas optics and Planck sources.
1280 : !$acc data copyin(kdist_lw, pmid_rad, pint_rad, &
1281 : !$acc t_rad, t_sfc, gas_concs_lw) &
1282 : !$acc copy(atm_optics_lw, atm_optics_lw%tau, &
1283 : !$acc sources_lw, sources_lw%lay_source, &
1284 : !$acc sources_lw%sfc_source, sources_lw%lev_source_inc, &
1285 : !$acc sources_lw%lev_source_dec, sources_lw%sfc_source_jac)
1286 : errmsg = kdist_lw%gas_optics( &
1287 : pmid_rad, pint_rad, t_rad, t_sfc, gas_concs_lw, &
1288 30960 : atm_optics_lw, sources_lw)
1289 : !$acc end data
1290 30960 : call stop_on_err(errmsg, sub, 'kdist_lw%gas_optics')
1291 :
1292 : ! Set LW aerosol optical properties in the aer_lw object.
1293 30960 : call rrtmgp_set_aer_lw(icall, state, pbuf, aer_lw)
1294 :
1295 : ! Increment the gas optics by the aerosol optics.
1296 : !$acc data copyin(atm_optics_lw, atm_optics_lw%tau, &
1297 : !$acc aer_lw, aer_lw%tau, &
1298 : !$acc cloud_lw, cloud_lw%tau, &
1299 : !$acc sources_lw, sources_lw%lay_source, &
1300 : !$acc sources_lw%sfc_source, sources_lw%lev_source_inc, &
1301 : !$acc sources_lw%lev_source_dec, sources_lw%sfc_source_Jac, &
1302 : !$acc emis_sfc) &
1303 : !$acc copy(flwc, flwc%flux_net, flwc%flux_up, flwc%flux_dn, &
1304 : !$acc flw, flw%flux_net, flw%flux_up, flw%flux_dn)
1305 30960 : errmsg = aer_lw%increment(atm_optics_lw)
1306 30960 : call stop_on_err(errmsg, sub, 'aer_lw%increment')
1307 :
1308 : ! Compute clear-sky LW fluxes
1309 30960 : errmsg = rte_lw(atm_optics_lw, top_at_1, sources_lw, emis_sfc, flwc)
1310 30960 : call stop_on_err(errmsg, sub, 'clear-sky rte_lw')
1311 :
1312 : ! Increment the gas+aerosol optics by the cloud optics.
1313 30960 : errmsg = cloud_lw%increment(atm_optics_lw)
1314 30960 : call stop_on_err(errmsg, sub, 'cloud_lw%increment')
1315 :
1316 : ! Compute all-sky LW fluxes
1317 30960 : errmsg = rte_lw(atm_optics_lw, top_at_1, sources_lw, emis_sfc, flw)
1318 30960 : call stop_on_err(errmsg, sub, 'all-sky rte_lw')
1319 : !$acc end data
1320 :
1321 : ! Transform RRTMGP outputs to CAM outputs and compute heating rates.
1322 30960 : call set_lw_diags()
1323 :
1324 30960 : if (write_output) then
1325 30960 : call radiation_output_lw(lchnk, ncol, icall, rd, pbuf, cam_out)
1326 : end if
1327 :
1328 : end if ! (active_calls(icall))
1329 : end do ! loop over diagnostic calcs (icall)
1330 :
1331 : else
1332 : ! LW calc not done. pbuf carries Q*dp across timesteps.
1333 : ! Convert to Q before calling radheat_tend.
1334 0 : qrl(:ncol,:) = qrl(:ncol,:) / state%pdel(:ncol,:)
1335 : end if ! if (dolw)
1336 :
1337 0 : deallocate( &
1338 0 : t_sfc, emis_sfc, toa_flux, sfac, t_rad, pmid_rad, pint_rad, &
1339 30960 : t_day, pmid_day, pint_day, coszrs_day, alb_dir, alb_dif)
1340 :
1341 : !================!
1342 : ! COSP simulator !
1343 : !================!
1344 :
1345 30960 : if (docosp) then
1346 :
1347 0 : emis(:,:) = 0._r8
1348 0 : emis(:ncol,:) = 1._r8 - exp(-cld_lw_abs_cloudsim(:ncol,:))
1349 0 : call outfld('EMIS', emis, pcols, lchnk)
1350 :
1351 : ! compute grid-box mean SW and LW snow optical depth for use by COSP
1352 0 : gb_snow_tau(:,:) = 0._r8
1353 0 : gb_snow_lw(:,:) = 0._r8
1354 0 : if (cldfsnow_idx > 0) then
1355 0 : do i = 1, ncol
1356 0 : do k = 1, pver
1357 0 : if (cldfsnow(i,k) > 0._r8) then
1358 :
1359 : ! Add graupel to snow tau for cosp
1360 0 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
1361 0 : gb_snow_tau(i,k) = snow_tau_cloudsim(i,k)*cldfsnow(i,k) + &
1362 0 : grau_tau_cloudsim(i,k)*cldfgrau(i,k)
1363 : gb_snow_lw(i,k) = snow_lw_abs_cloudsim(i,k)*cldfsnow(i,k) + &
1364 0 : grau_lw_abs_cloudsim(i,k)*cldfgrau(i,k)
1365 : else
1366 0 : gb_snow_tau(i,k) = snow_tau_cloudsim(i,k)*cldfsnow(i,k)
1367 0 : gb_snow_lw(i,k) = snow_lw_abs_cloudsim(i,k)*cldfsnow(i,k)
1368 : end if
1369 : end if
1370 : end do
1371 : end do
1372 : end if
1373 :
1374 : ! advance counter for this timestep (chunk dimension required for thread safety)
1375 0 : cosp_cnt(lchnk) = cosp_cnt(lchnk) + 1
1376 :
1377 : ! if counter is the same as cosp_nradsteps, run cosp and reset counter
1378 0 : if (cosp_nradsteps .eq. cosp_cnt(lchnk)) then
1379 :
1380 : ! N.B.: For snow optical properties, the GRID-BOX MEAN shortwave and longwave
1381 : ! optical depths are passed.
1382 : call cospsimulator_intr_run( &
1383 : state, pbuf, cam_in, emis, coszrs, &
1384 : cld_swtau_in=cld_tau_cloudsim, snow_tau_in=gb_snow_tau, &
1385 0 : snow_emis_in=gb_snow_lw)
1386 0 : cosp_cnt(lchnk) = 0
1387 : end if
1388 : end if ! docosp
1389 :
1390 : else
1391 : ! Radiative flux calculations not done. The quantity Q*dp is carried by the
1392 : ! physics buffer across timesteps. It must be converted to Q (dry static energy
1393 : ! tendency) before being passed to radheat_tend.
1394 48108240 : qrs(:ncol,:) = qrs(:ncol,:) / state%pdel(:ncol,:)
1395 48108240 : qrl(:ncol,:) = qrl(:ncol,:) / state%pdel(:ncol,:)
1396 :
1397 : end if ! if (dosw .or. dolw) then
1398 :
1399 : ! Output for PORT: Parallel Offline Radiative Transport
1400 61920 : call rad_data_write(pbuf, state, cam_in, coszrs)
1401 :
1402 : ! Compute net radiative heating tendency. Note that the WACCM version
1403 : ! of radheat_tend merges upper atmosphere heating rates with those calculated
1404 : ! by RRTMGP.
1405 : call radheat_tend(state, pbuf, ptend, qrl, qrs, fsns, &
1406 61920 : fsnt, flns, flnt, cam_in%asdir, net_flx)
1407 :
1408 61920 : if (write_output) then
1409 : ! Compute heating rate for dtheta/dt
1410 5820480 : do k = 1, pver
1411 96216480 : do i = 1, ncol
1412 96154560 : ftem(i,k) = (qrs(i,k) + qrl(i,k))/cpair * (1.e5_r8/state%pmid(i,k))**cappa
1413 : end do
1414 : end do
1415 61920 : call outfld('HR', ftem, pcols, lchnk)
1416 : end if
1417 :
1418 : ! The radiative heating rates are carried in the physics buffer across timesteps
1419 : ! as Q*dp (for energy conservation).
1420 96216480 : qrs(:ncol,:) = qrs(:ncol,:) * state%pdel(:ncol,:)
1421 96216480 : qrl(:ncol,:) = qrl(:ncol,:) * state%pdel(:ncol,:)
1422 :
1423 61920 : if (.not. present(rd_out)) then
1424 61920 : deallocate(rd)
1425 : end if
1426 61920 : call free_optics_sw(atm_optics_sw)
1427 61920 : call free_optics_sw(cloud_sw)
1428 61920 : call free_optics_sw(aer_sw)
1429 61920 : call free_fluxes(fsw)
1430 61920 : call free_fluxes(fswc)
1431 :
1432 61920 : call sources_lw%finalize()
1433 61920 : call free_optics_lw(cloud_lw)
1434 61920 : call free_optics_lw(aer_lw)
1435 61920 : call free_fluxes(flw)
1436 123840 : call free_fluxes(flwc)
1437 :
1438 : !-------------------------------------------------------------------------------
1439 : contains
1440 : !-------------------------------------------------------------------------------
1441 :
1442 30960 : subroutine set_sw_diags()
1443 :
1444 : ! Transform RRTMGP output for CAM and compute heating rates.
1445 : ! SW fluxes from RRTMGP are on daylight columns only, so expand to
1446 : ! full chunks for output to CAM history.
1447 :
1448 : integer :: i
1449 : real(r8), dimension(size(fsw%bnd_flux_dn,1), &
1450 : size(fsw%bnd_flux_dn,2), &
1451 61920 : size(fsw%bnd_flux_dn,3)) :: flux_dn_diffuse
1452 : !-------------------------------------------------------------------------
1453 :
1454 : ! Initialize to provide 0.0 values for night columns.
1455 30960 : fns = 0._r8 ! net sw flux
1456 30960 : fcns = 0._r8 ! net sw clearsky flux
1457 526320 : fsds = 0._r8 ! downward sw flux at surface
1458 526320 : rd%fsdsc = 0._r8 ! downward sw clearsky flux at surface
1459 526320 : rd%fsutoa = 0._r8 ! upward sw flux at TOA
1460 526320 : rd%fsntoa = 0._r8 ! net sw at TOA
1461 526320 : rd%fsntoac = 0._r8 ! net sw clearsky flux at TOA
1462 526320 : rd%solin = 0._r8 ! solar irradiance at TOA
1463 49505040 : rd%flux_sw_up = 0._r8
1464 49505040 : rd%flux_sw_dn = 0._r8
1465 49505040 : rd%flux_sw_clr_up = 0._r8
1466 49505040 : rd%flux_sw_clr_dn = 0._r8
1467 :
1468 48978720 : qrs = 0._r8
1469 526320 : fsns = 0._r8
1470 526320 : fsnt = 0._r8
1471 48978720 : rd%qrsc = 0._r8
1472 526320 : rd%fsnsc = 0._r8
1473 526320 : rd%fsntc = 0._r8
1474 :
1475 273960 : do i = 1, nday
1476 23085000 : fns(idxday(i),ktopcam:) = fsw%flux_net(i, ktoprad:)
1477 23085000 : fcns(idxday(i),ktopcam:) = fswc%flux_net(i,ktoprad:)
1478 243000 : fsds(idxday(i)) = fsw%flux_dn(i, nlay+1)
1479 243000 : rd%fsdsc(idxday(i)) = fswc%flux_dn(i, nlay+1)
1480 243000 : rd%fsutoa(idxday(i)) = fsw%flux_up(i, 1)
1481 243000 : rd%fsntoa(idxday(i)) = fsw%flux_net(i, 1)
1482 243000 : rd%fsntoac(idxday(i)) = fswc%flux_net(i, 1)
1483 243000 : rd%solin(idxday(i)) = fswc%flux_dn(i, 1)
1484 23085000 : rd%flux_sw_up(idxday(i),ktopcam:) = fsw%flux_up(i,ktoprad:)
1485 23085000 : rd%flux_sw_dn(idxday(i),ktopcam:) = fsw%flux_dn(i,ktoprad:)
1486 23085000 : rd%flux_sw_clr_up(idxday(i),ktopcam:) = fswc%flux_up(i,ktoprad:)
1487 23115960 : rd%flux_sw_clr_dn(idxday(i),ktopcam:) = fswc%flux_dn(i,ktoprad:)
1488 : end do
1489 :
1490 : ! Compute heating rate as a dry static energy tendency.
1491 30960 : call heating_rate('SW', ncol, fns, qrs)
1492 30960 : call heating_rate('SW', ncol, fcns, rd%qrsc)
1493 :
1494 516960 : fsns(:ncol) = fns(:ncol,pverp) ! net sw flux at surface
1495 516960 : fsnt(:ncol) = fns(:ncol,ktopcam) ! net sw flux at top-of-model (w/o extra layer)
1496 547920 : rd%fsnsc(:ncol) = fcns(:ncol,pverp) ! net sw clearsky flux at surface
1497 547920 : rd%fsntc(:ncol) = fcns(:ncol,ktopcam) ! net sw clearsky flux at top
1498 :
1499 516960 : cam_out%netsw(:ncol) = fsns(:ncol)
1500 :
1501 : ! Output fluxes at 200 mb
1502 30960 : call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fns, rd%fsn200)
1503 30960 : call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fcns, rd%fsn200c)
1504 30960 : if (hist_fld_active('FSNR')) then
1505 0 : do i = 1,ncol
1506 0 : call vertinterp(1, 1, pverp, state%pint(i,:), p_trop(i), fns(i,:), rd%fsnr(i))
1507 : end do
1508 : end if
1509 :
1510 30960 : if (spectralflux) then
1511 0 : su = 0._r8
1512 0 : sd = 0._r8
1513 0 : do i = 1, nday
1514 0 : su(idxday(i),ktopcam:,:) = fsw%bnd_flux_up(i,ktoprad:,:)
1515 0 : sd(idxday(i),ktopcam:,:) = fsw%bnd_flux_dn(i,ktoprad:,:)
1516 : end do
1517 : end if
1518 :
1519 : ! Export surface fluxes
1520 : ! sols(pcols) Direct solar rad on surface (< 0.7)
1521 : ! soll(pcols) Direct solar rad on surface (>= 0.7)
1522 : ! RRTMGP: Near-IR bands (1-10), 820-16000 cm-1, 0.625-12.195 microns
1523 : ! Put half of band 10 in each of the UV/visible and near-IR values,
1524 : ! since this band straddles 0.7 microns:
1525 : ! UV/visible bands 10-13, 16000-50000 cm-1, 0.200-0.625 micron
1526 :
1527 : ! reset fluxes
1528 526320 : cam_out%sols = 0.0_r8
1529 526320 : cam_out%soll = 0.0_r8
1530 526320 : cam_out%solsd = 0.0_r8
1531 526320 : cam_out%solld = 0.0_r8
1532 :
1533 : ! Calculate diffuse flux from total and direct
1534 364831200 : flux_dn_diffuse = fsw%bnd_flux_dn - fsw%bnd_flux_dn_dir
1535 :
1536 273960 : do i = 1, nday
1537 729000 : cam_out%soll(idxday(i)) = sum(fsw%bnd_flux_dn_dir(i,nlay+1,1:9)) &
1538 3159000 : + 0.5_r8 * fsw%bnd_flux_dn_dir(i,nlay+1,10)
1539 :
1540 243000 : cam_out%sols(idxday(i)) = 0.5_r8 * fsw%bnd_flux_dn_dir(i,nlay+1,10) &
1541 1458000 : + sum(fsw%bnd_flux_dn_dir(i,nlay+1,11:14))
1542 :
1543 486000 : cam_out%solld(idxday(i)) = sum(flux_dn_diffuse(i,nlay+1,1:9)) &
1544 2916000 : + 0.5_r8 * flux_dn_diffuse(i,nlay+1,10)
1545 :
1546 243000 : cam_out%solsd(idxday(i)) = 0.5_r8 * flux_dn_diffuse(i, nlay+1, 10) &
1547 1731960 : + sum(flux_dn_diffuse(i,nlay+1,11:14))
1548 : end do
1549 :
1550 61920 : end subroutine set_sw_diags
1551 :
1552 : !-------------------------------------------------------------------------------
1553 :
1554 30960 : subroutine set_lw_diags()
1555 :
1556 : ! Set CAM LW diagnostics
1557 : !----------------------------------------------------------------------------
1558 :
1559 30960 : fnl = 0._r8
1560 30960 : fcnl = 0._r8
1561 :
1562 : ! RTE-RRTMGP convention for net is (down - up) **CAM assumes (up - down) !!
1563 48625200 : fnl(:ncol,ktopcam:) = -1._r8 * flw%flux_net( :, ktoprad:)
1564 48625200 : fcnl(:ncol,ktopcam:) = -1._r8 * flwc%flux_net( :, ktoprad:)
1565 :
1566 48625200 : rd%flux_lw_up(:ncol,ktopcam:) = flw%flux_up( :, ktoprad:)
1567 48625200 : rd%flux_lw_clr_up(:ncol,ktopcam:) = flwc%flux_up(:, ktoprad:)
1568 48625200 : rd%flux_lw_dn(:ncol,ktopcam:) = flw%flux_dn( :, ktoprad:)
1569 48625200 : rd%flux_lw_clr_dn(:ncol,ktopcam:) = flwc%flux_dn(:, ktoprad:)
1570 :
1571 30960 : call heating_rate('LW', ncol, fnl, qrl)
1572 30960 : call heating_rate('LW', ncol, fcnl, rd%qrlc)
1573 :
1574 516960 : flns(:ncol) = fnl(:ncol, pverp)
1575 516960 : flnt(:ncol) = fnl(:ncol, ktopcam)
1576 :
1577 547920 : rd%flnsc(:ncol) = fcnl(:ncol, pverp)
1578 547920 : rd%flntc(:ncol) = fcnl(:ncol, ktopcam) ! net lw flux at top-of-model
1579 :
1580 516960 : cam_out%flwds(:ncol) = flw%flux_dn(:, nlay+1)
1581 516960 : rd%fldsc(:ncol) = flwc%flux_dn(:, nlay+1)
1582 :
1583 516960 : rd%flut(:ncol) = flw%flux_up(:, ktoprad)
1584 516960 : rd%flutc(:ncol) = flwc%flux_up(:, ktoprad)
1585 :
1586 : ! Output fluxes at 200 mb
1587 30960 : call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fnl, rd%fln200)
1588 30960 : call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fcnl, rd%fln200c)
1589 30960 : if (hist_fld_active('FLNR')) then
1590 0 : do i = 1,ncol
1591 0 : call vertinterp(1, 1, pverp, state%pint(i,:), p_trop(i), fnl(i,:), rd%flnr(i))
1592 : end do
1593 : end if
1594 :
1595 30960 : if (spectralflux) then
1596 0 : lu = 0._r8
1597 0 : ld = 0._r8
1598 0 : lu(:ncol, ktopcam:, :) = flw%bnd_flux_up(:, ktoprad:, :)
1599 0 : ld(:ncol, ktopcam:, :) = flw%bnd_flux_dn(:, ktoprad:, :)
1600 : end if
1601 :
1602 30960 : end subroutine set_lw_diags
1603 :
1604 : !-------------------------------------------------------------------------------
1605 :
1606 123840 : subroutine heating_rate(type, ncol, flux_net, hrate)
1607 :
1608 : ! Compute heating rate as a dry static energy tendency
1609 :
1610 : ! arguments
1611 : character(2), intent(in) :: type ! either LW or SW
1612 : integer, intent(in) :: ncol
1613 : real(r8), intent(in) :: flux_net(pcols,pverp) ! W/m^2
1614 : real(r8), intent(out) :: hrate(pcols,pver) ! J/kg/s
1615 :
1616 : ! local vars
1617 : integer :: k
1618 :
1619 : ! Initialize for layers where RRTMGP is not providing fluxes.
1620 123840 : hrate = 0.0_r8
1621 :
1622 61920 : select case (type)
1623 : case ('LW')
1624 :
1625 5820480 : do k = ktopcam, pver
1626 : ! (flux divergence as bottom-MINUS-top) * g/dp
1627 11517120 : hrate(:ncol,k) = (flux_net(:ncol,k+1) - flux_net(:ncol,k)) * &
1628 107733600 : gravit * state%rpdel(:ncol,k)
1629 : end do
1630 :
1631 : case ('SW')
1632 :
1633 5882400 : do k = ktopcam, pver
1634 : ! top - bottom
1635 11517120 : hrate(:ncol,k) = (flux_net(:ncol,k) - flux_net(:ncol,k+1)) * &
1636 107733600 : gravit * state%rpdel(:ncol,k)
1637 : end do
1638 :
1639 : end select
1640 :
1641 123840 : end subroutine heating_rate
1642 :
1643 : !----------------------------------------------------------------------------
1644 : ! -- end contains statement of radiation_tend --
1645 : !----------------------------------------------------------------------------
1646 : end subroutine radiation_tend
1647 :
1648 : !===============================================================================
1649 :
1650 30960 : subroutine radiation_output_sw(lchnk, ncol, icall, rd, pbuf, cam_out)
1651 :
1652 : ! Dump shortwave radiation information to history buffer.
1653 :
1654 : integer , intent(in) :: lchnk
1655 : integer, intent(in) :: ncol
1656 : integer, intent(in) :: icall
1657 : type(rad_out_t), intent(in) :: rd
1658 : type(physics_buffer_desc), pointer :: pbuf(:)
1659 : type(cam_out_t), intent(in) :: cam_out
1660 :
1661 : ! local variables
1662 30960 : real(r8), pointer :: qrs(:,:)
1663 30960 : real(r8), pointer :: fsnt(:)
1664 30960 : real(r8), pointer :: fsns(:)
1665 30960 : real(r8), pointer :: fsds(:)
1666 :
1667 : real(r8) :: ftem(pcols)
1668 : !----------------------------------------------------------------------------
1669 :
1670 30960 : call pbuf_get_field(pbuf, qrs_idx, qrs)
1671 30960 : call pbuf_get_field(pbuf, fsnt_idx, fsnt)
1672 30960 : call pbuf_get_field(pbuf, fsns_idx, fsns)
1673 30960 : call pbuf_get_field(pbuf, fsds_idx, fsds)
1674 :
1675 30960 : call outfld('SOLIN'//diag(icall), rd%solin, pcols, lchnk)
1676 :
1677 : ! QRS is output as temperature tendency.
1678 48108240 : call outfld('QRS'//diag(icall), qrs(:ncol,:)/cpair, ncol, lchnk)
1679 48108240 : call outfld('QRSC'//diag(icall), rd%qrsc(:ncol,:)/cpair, ncol, lchnk)
1680 :
1681 30960 : call outfld('FSNT'//diag(icall), fsnt, pcols, lchnk)
1682 30960 : call outfld('FSNTC'//diag(icall), rd%fsntc, pcols, lchnk)
1683 30960 : call outfld('FSNTOA'//diag(icall), rd%fsntoa, pcols, lchnk)
1684 30960 : call outfld('FSNTOAC'//diag(icall), rd%fsntoac, pcols, lchnk)
1685 :
1686 516960 : ftem(:ncol) = rd%fsntoa(:ncol) - rd%fsntoac(:ncol)
1687 30960 : call outfld('SWCF'//diag(icall), ftem, pcols, lchnk)
1688 :
1689 30960 : call outfld('FSUTOA'//diag(icall), rd%fsutoa, pcols, lchnk)
1690 :
1691 30960 : call outfld('FSNIRTOA'//diag(icall), rd%fsnirt, pcols, lchnk)
1692 30960 : call outfld('FSNRTOAC'//diag(icall), rd%fsnrtc, pcols, lchnk)
1693 30960 : call outfld('FSNRTOAS'//diag(icall), rd%fsnirtsq, pcols, lchnk)
1694 :
1695 30960 : call outfld('FSN200'//diag(icall), rd%fsn200, pcols, lchnk)
1696 30960 : call outfld('FSN200C'//diag(icall), rd%fsn200c, pcols, lchnk)
1697 :
1698 30960 : call outfld('FSNR'//diag(icall), rd%fsnr, pcols, lchnk)
1699 :
1700 30960 : call outfld('SOLS'//diag(icall), cam_out%sols, pcols, lchnk)
1701 30960 : call outfld('SOLL'//diag(icall), cam_out%soll, pcols, lchnk)
1702 30960 : call outfld('SOLSD'//diag(icall), cam_out%solsd, pcols, lchnk)
1703 30960 : call outfld('SOLLD'//diag(icall), cam_out%solld, pcols, lchnk)
1704 :
1705 30960 : call outfld('FSNS'//diag(icall), fsns, pcols, lchnk)
1706 30960 : call outfld('FSNSC'//diag(icall), rd%fsnsc, pcols, lchnk)
1707 :
1708 30960 : call outfld('FSDS'//diag(icall), fsds, pcols, lchnk)
1709 30960 : call outfld('FSDSC'//diag(icall), rd%fsdsc, pcols, lchnk)
1710 :
1711 30960 : call outfld('FUS'//diag(icall), rd%flux_sw_up, pcols, lchnk)
1712 30960 : call outfld('FUSC'//diag(icall), rd%flux_sw_clr_up, pcols, lchnk)
1713 30960 : call outfld('FDS'//diag(icall), rd%flux_sw_dn, pcols, lchnk)
1714 30960 : call outfld('FDSC'//diag(icall), rd%flux_sw_clr_dn, pcols, lchnk)
1715 :
1716 30960 : end subroutine radiation_output_sw
1717 :
1718 : !===============================================================================
1719 :
1720 30960 : subroutine radiation_output_cld(lchnk, rd)
1721 :
1722 : ! Dump shortwave cloud optics information to history buffer.
1723 :
1724 : integer , intent(in) :: lchnk
1725 : type(rad_out_t), intent(in) :: rd
1726 : !----------------------------------------------------------------------------
1727 :
1728 30960 : call outfld('TOT_CLD_VISTAU', rd%tot_cld_vistau, pcols, lchnk)
1729 30960 : call outfld('TOT_ICLD_VISTAU', rd%tot_icld_vistau, pcols, lchnk)
1730 30960 : call outfld('LIQ_ICLD_VISTAU', rd%liq_icld_vistau, pcols, lchnk)
1731 30960 : call outfld('ICE_ICLD_VISTAU', rd%ice_icld_vistau, pcols, lchnk)
1732 30960 : if (cldfsnow_idx > 0) then
1733 30960 : call outfld('SNOW_ICLD_VISTAU', rd%snow_icld_vistau, pcols, lchnk)
1734 : endif
1735 30960 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
1736 0 : call outfld('GRAU_ICLD_VISTAU', rd%grau_icld_vistau , pcols, lchnk)
1737 : endif
1738 :
1739 30960 : end subroutine radiation_output_cld
1740 :
1741 : !===============================================================================
1742 :
1743 30960 : subroutine radiation_output_lw(lchnk, ncol, icall, rd, pbuf, cam_out)
1744 :
1745 : ! Dump longwave radiation information to history buffer
1746 :
1747 : integer, intent(in) :: lchnk
1748 : integer, intent(in) :: ncol
1749 : integer, intent(in) :: icall ! icall=0 for climate diagnostics
1750 : type(rad_out_t), intent(in) :: rd
1751 : type(physics_buffer_desc), pointer :: pbuf(:)
1752 : type(cam_out_t), intent(in) :: cam_out
1753 :
1754 : ! local variables
1755 30960 : real(r8), pointer :: qrl(:,:)
1756 30960 : real(r8), pointer :: flnt(:)
1757 30960 : real(r8), pointer :: flns(:)
1758 :
1759 : real(r8) :: ftem(pcols)
1760 : !----------------------------------------------------------------------------
1761 :
1762 30960 : call pbuf_get_field(pbuf, qrl_idx, qrl)
1763 30960 : call pbuf_get_field(pbuf, flnt_idx, flnt)
1764 30960 : call pbuf_get_field(pbuf, flns_idx, flns)
1765 :
1766 48108240 : call outfld('QRL'//diag(icall), qrl(:ncol,:)/cpair, ncol, lchnk)
1767 48108240 : call outfld('QRLC'//diag(icall), rd%qrlc(:ncol,:)/cpair, ncol, lchnk)
1768 :
1769 30960 : call outfld('FLNT'//diag(icall), flnt, pcols, lchnk)
1770 30960 : call outfld('FLNTC'//diag(icall), rd%flntc, pcols, lchnk)
1771 :
1772 30960 : call outfld('FLUT'//diag(icall), rd%flut, pcols, lchnk)
1773 30960 : call outfld('FLUTC'//diag(icall), rd%flutc, pcols, lchnk)
1774 :
1775 516960 : ftem(:ncol) = rd%flutc(:ncol) - rd%flut(:ncol)
1776 30960 : call outfld('LWCF'//diag(icall), ftem, pcols, lchnk)
1777 :
1778 30960 : call outfld('FLN200'//diag(icall), rd%fln200, pcols, lchnk)
1779 30960 : call outfld('FLN200C'//diag(icall), rd%fln200c, pcols, lchnk)
1780 :
1781 30960 : call outfld('FLNR'//diag(icall), rd%flnr, pcols, lchnk)
1782 :
1783 30960 : call outfld('FLNS'//diag(icall), flns, pcols, lchnk)
1784 30960 : call outfld('FLNSC'//diag(icall), rd%flnsc, pcols, lchnk)
1785 :
1786 30960 : call outfld('FLDS'//diag(icall), cam_out%flwds, pcols, lchnk)
1787 30960 : call outfld('FLDSC'//diag(icall), rd%fldsc, pcols, lchnk)
1788 :
1789 30960 : call outfld('FDL'//diag(icall), rd%flux_lw_dn, pcols, lchnk)
1790 30960 : call outfld('FDLC'//diag(icall), rd%flux_lw_clr_dn, pcols, lchnk)
1791 30960 : call outfld('FUL'//diag(icall), rd%flux_lw_up, pcols, lchnk)
1792 30960 : call outfld('FULC'//diag(icall), rd%flux_lw_clr_up, pcols, lchnk)
1793 :
1794 30960 : end subroutine radiation_output_lw
1795 :
1796 : !===============================================================================
1797 :
1798 3072 : subroutine coefs_init(coefs_file, available_gases, kdist)
1799 :
1800 : ! Read data from coefficients file. Initialize the kdist object.
1801 : ! available_gases object provides the gas names that CAM provides.
1802 :
1803 : ! arguments
1804 : character(len=*), intent(in) :: coefs_file
1805 : class(ty_gas_concs), intent(in) :: available_gases
1806 : class(ty_gas_optics_rrtmgp), intent(out) :: kdist
1807 :
1808 : ! local variables
1809 : type(file_desc_t) :: fh ! pio file handle
1810 : character(len=cl) :: locfn ! path to file on local storage
1811 :
1812 : ! File dimensions
1813 : integer :: &
1814 : absorber, &
1815 : atmos_layer, &
1816 : bnd, &
1817 : pressure, &
1818 : temperature, &
1819 : absorber_ext, &
1820 : pressure_interp, &
1821 : mixing_fraction, &
1822 : gpt, &
1823 : temperature_Planck
1824 :
1825 : integer :: i
1826 : integer :: did, vid
1827 : integer :: ierr, istat
1828 :
1829 3072 : character(32), dimension(:), allocatable :: gas_names
1830 3072 : integer, dimension(:,:,:), allocatable :: key_species
1831 3072 : integer, dimension(:,:), allocatable :: band2gpt
1832 3072 : real(r8), dimension(:,:), allocatable :: band_lims_wavenum
1833 3072 : real(r8), dimension(:), allocatable :: press_ref, temp_ref
1834 : real(r8) :: press_ref_trop, temp_ref_t, temp_ref_p
1835 3072 : real(r8), dimension(:,:,:), allocatable :: vmr_ref
1836 3072 : real(r8), dimension(:,:,:,:), allocatable :: kmajor
1837 3072 : real(r8), dimension(:,:,:), allocatable :: kminor_lower, kminor_upper
1838 3072 : real(r8), dimension(:,:), allocatable :: totplnk
1839 3072 : real(r8), dimension(:,:,:,:), allocatable :: planck_frac
1840 3072 : real(r8), dimension(:), allocatable :: solar_src_quiet, solar_src_facular, solar_src_sunspot
1841 : real(r8) :: tsi_default
1842 3072 : real(r8), dimension(:,:,:), allocatable :: rayl_lower, rayl_upper
1843 3072 : character(len=32), dimension(:), allocatable :: gas_minor, &
1844 3072 : identifier_minor, &
1845 3072 : minor_gases_lower, &
1846 3072 : minor_gases_upper, &
1847 3072 : scaling_gas_lower, &
1848 3072 : scaling_gas_upper
1849 3072 : integer, dimension(:,:), allocatable :: minor_limits_gpt_lower, &
1850 3072 : minor_limits_gpt_upper
1851 : ! Send these to RRTMGP as logicals,
1852 : ! but they have to be read from the netCDF as integers
1853 3072 : logical, dimension(:), allocatable :: minor_scales_with_density_lower, &
1854 3072 : minor_scales_with_density_upper
1855 3072 : logical, dimension(:), allocatable :: scale_by_complement_lower, &
1856 3072 : scale_by_complement_upper
1857 3072 : integer, dimension(:), allocatable :: int2log ! use this to convert integer-to-logical.
1858 3072 : integer, dimension(:), allocatable :: kminor_start_lower, kminor_start_upper
1859 3072 : real(r8), dimension(:,:), allocatable :: optimal_angle_fit
1860 : real(r8) :: mg_default, sb_default
1861 :
1862 : integer :: pairs, &
1863 : minorabsorbers, &
1864 : minor_absorber_intervals_lower, &
1865 : minor_absorber_intervals_upper, &
1866 : contributors_lower, &
1867 : contributors_upper, &
1868 : fit_coeffs
1869 :
1870 : character(len=128) :: error_msg
1871 : character(len=*), parameter :: sub = 'coefs_init'
1872 : !----------------------------------------------------------------------------
1873 :
1874 : ! Open file
1875 3072 : call getfil(coefs_file, locfn, 0)
1876 3072 : call cam_pio_openfile(fh, locfn, PIO_NOWRITE)
1877 :
1878 3072 : call pio_seterrorhandling(fh, PIO_BCAST_ERROR)
1879 :
1880 : ! Get dimensions
1881 :
1882 3072 : ierr = pio_inq_dimid(fh, 'absorber', did)
1883 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': absorber not found')
1884 3072 : ierr = pio_inq_dimlen(fh, did, absorber)
1885 :
1886 3072 : ierr = pio_inq_dimid(fh, 'atmos_layer', did)
1887 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': atmos_layer not found')
1888 3072 : ierr = pio_inq_dimlen(fh, did, atmos_layer)
1889 :
1890 3072 : ierr = pio_inq_dimid(fh, 'bnd', did)
1891 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': bnd not found')
1892 3072 : ierr = pio_inq_dimlen(fh, did, bnd)
1893 :
1894 3072 : ierr = pio_inq_dimid(fh, 'pressure', did)
1895 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': pressure not found')
1896 3072 : ierr = pio_inq_dimlen(fh, did, pressure)
1897 :
1898 3072 : ierr = pio_inq_dimid(fh, 'temperature', did)
1899 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': temperature not found')
1900 3072 : ierr = pio_inq_dimlen(fh, did, temperature)
1901 :
1902 3072 : ierr = pio_inq_dimid(fh, 'absorber_ext', did)
1903 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': absorber_ext not found')
1904 3072 : ierr = pio_inq_dimlen(fh, did, absorber_ext)
1905 :
1906 3072 : ierr = pio_inq_dimid(fh, 'pressure_interp', did)
1907 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': pressure_interp not found')
1908 3072 : ierr = pio_inq_dimlen(fh, did, pressure_interp)
1909 :
1910 3072 : ierr = pio_inq_dimid(fh, 'mixing_fraction', did)
1911 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': mixing_fraction not found')
1912 3072 : ierr = pio_inq_dimlen(fh, did, mixing_fraction)
1913 :
1914 3072 : ierr = pio_inq_dimid(fh, 'gpt', did)
1915 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': gpt not found')
1916 3072 : ierr = pio_inq_dimlen(fh, did, gpt)
1917 :
1918 3072 : temperature_Planck = 0
1919 3072 : ierr = pio_inq_dimid(fh, 'temperature_Planck', did)
1920 3072 : if (ierr == PIO_NOERR) then
1921 1536 : ierr = pio_inq_dimlen(fh, did, temperature_Planck)
1922 : end if
1923 3072 : ierr = pio_inq_dimid(fh, 'pair', did)
1924 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': pair not found')
1925 3072 : ierr = pio_inq_dimlen(fh, did, pairs)
1926 3072 : ierr = pio_inq_dimid(fh, 'minor_absorber', did)
1927 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': minor_absorber not found')
1928 3072 : ierr = pio_inq_dimlen(fh, did, minorabsorbers)
1929 3072 : ierr = pio_inq_dimid(fh, 'minor_absorber_intervals_lower', did)
1930 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': minor_absorber_intervals_lower not found')
1931 3072 : ierr = pio_inq_dimlen(fh, did, minor_absorber_intervals_lower)
1932 3072 : ierr = pio_inq_dimid(fh, 'minor_absorber_intervals_upper', did)
1933 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': minor_absorber_intervals_upper not found')
1934 3072 : ierr = pio_inq_dimlen(fh, did, minor_absorber_intervals_upper)
1935 3072 : ierr = pio_inq_dimid(fh, 'contributors_lower', did)
1936 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': contributors_lower not found')
1937 3072 : ierr = pio_inq_dimlen(fh, did, contributors_lower)
1938 3072 : ierr = pio_inq_dimid(fh, 'contributors_upper', did)
1939 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': contributors_upper not found')
1940 3072 : ierr = pio_inq_dimlen(fh, did, contributors_upper)
1941 :
1942 3072 : ierr = pio_inq_dimid(fh, 'fit_coeffs', did)
1943 3072 : if (ierr == PIO_NOERR) then
1944 1536 : ierr = pio_inq_dimlen(fh, did, fit_coeffs)
1945 : end if
1946 :
1947 : ! Get variables
1948 :
1949 : ! names of absorbing gases
1950 9216 : allocate(gas_names(absorber), stat=istat)
1951 3072 : call handle_allocate_error(istat, sub, 'gas_names')
1952 3072 : ierr = pio_inq_varid(fh, 'gas_names', vid)
1953 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': gas_names not found')
1954 3072 : ierr = pio_get_var(fh, vid, gas_names)
1955 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading gas_names')
1956 :
1957 : ! key species pair for each band
1958 12288 : allocate(key_species(2,atmos_layer,bnd), stat=istat)
1959 3072 : call handle_allocate_error(istat, sub, 'key_species')
1960 3072 : ierr = pio_inq_varid(fh, 'key_species', vid)
1961 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': key_species not found')
1962 3072 : ierr = pio_get_var(fh, vid, key_species)
1963 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading key_species')
1964 :
1965 : ! beginning and ending gpoint for each band
1966 9216 : allocate(band2gpt(2,bnd), stat=istat)
1967 3072 : call handle_allocate_error(istat, sub, 'band2gpt')
1968 3072 : ierr = pio_inq_varid(fh, 'bnd_limits_gpt', vid)
1969 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': bnd_limits_gpt not found')
1970 3072 : ierr = pio_get_var(fh, vid, band2gpt)
1971 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading bnd_limits_gpt')
1972 :
1973 : ! beginning and ending wavenumber for each band
1974 9216 : allocate(band_lims_wavenum(2,bnd), stat=istat)
1975 3072 : call handle_allocate_error(istat, sub, 'band_lims_wavenum')
1976 3072 : ierr = pio_inq_varid(fh, 'bnd_limits_wavenumber', vid)
1977 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': bnd_limits_wavenumber not found')
1978 3072 : ierr = pio_get_var(fh, vid, band_lims_wavenum)
1979 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading bnd_limits_wavenumber')
1980 :
1981 : ! pressures [hPa] for reference atmosphere; press_ref(# reference layers)
1982 9216 : allocate(press_ref(pressure), stat=istat)
1983 3072 : call handle_allocate_error(istat, sub, 'press_ref')
1984 3072 : ierr = pio_inq_varid(fh, 'press_ref', vid)
1985 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': press_ref not found')
1986 3072 : ierr = pio_get_var(fh, vid, press_ref)
1987 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading press_ref')
1988 :
1989 : ! reference pressure separating the lower and upper atmosphere
1990 3072 : ierr = pio_inq_varid(fh, 'press_ref_trop', vid)
1991 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': press_ref_trop not found')
1992 3072 : ierr = pio_get_var(fh, vid, press_ref_trop)
1993 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading press_ref_trop')
1994 :
1995 : ! temperatures [K] for reference atmosphere; temp_ref(# reference layers)
1996 9216 : allocate(temp_ref(temperature), stat=istat)
1997 3072 : call handle_allocate_error(istat, sub, 'temp_ref')
1998 3072 : ierr = pio_inq_varid(fh, 'temp_ref', vid)
1999 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': temp_ref not found')
2000 3072 : ierr = pio_get_var(fh, vid, temp_ref)
2001 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading temp_ref')
2002 :
2003 : ! standard spectroscopic reference temperature [K]
2004 3072 : ierr = pio_inq_varid(fh, 'absorption_coefficient_ref_T', vid)
2005 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': absorption_coefficient_ref_T not found')
2006 3072 : ierr = pio_get_var(fh, vid, temp_ref_t)
2007 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading absorption_coefficient_ref_T')
2008 :
2009 : ! standard spectroscopic reference pressure [hPa]
2010 3072 : ierr = pio_inq_varid(fh, 'absorption_coefficient_ref_P', vid)
2011 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': absorption_coefficient_ref_P not found')
2012 3072 : ierr = pio_get_var(fh, vid, temp_ref_p)
2013 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading absorption_coefficient_ref_P')
2014 :
2015 : ! volume mixing ratios for reference atmosphere
2016 15360 : allocate(vmr_ref(atmos_layer, absorber_ext, temperature), stat=istat)
2017 3072 : call handle_allocate_error(istat, sub, 'vmr_ref')
2018 3072 : ierr = pio_inq_varid(fh, 'vmr_ref', vid)
2019 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': vmr_ref not found')
2020 3072 : ierr = pio_get_var(fh, vid, vmr_ref)
2021 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading vmr_ref')
2022 :
2023 : ! absorption coefficients due to major absorbing gases
2024 18432 : allocate(kmajor(gpt,mixing_fraction,pressure_interp,temperature), stat=istat)
2025 3072 : call handle_allocate_error(istat, sub, 'kmajor')
2026 3072 : ierr = pio_inq_varid(fh, 'kmajor', vid)
2027 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': kmajor not found')
2028 3072 : ierr = pio_get_var(fh, vid, kmajor)
2029 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading kmajor')
2030 :
2031 : ! absorption coefficients due to minor absorbing gases in lower part of atmosphere
2032 15360 : allocate(kminor_lower(contributors_lower, mixing_fraction, temperature), stat=istat)
2033 3072 : call handle_allocate_error(istat, sub, 'kminor_lower')
2034 3072 : ierr = pio_inq_varid(fh, 'kminor_lower', vid)
2035 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': kminor_lower not found')
2036 3072 : ierr = pio_get_var(fh, vid, kminor_lower)
2037 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading kminor_lower')
2038 :
2039 : ! absorption coefficients due to minor absorbing gases in upper part of atmosphere
2040 15360 : allocate(kminor_upper(contributors_upper, mixing_fraction, temperature), stat=istat)
2041 3072 : call handle_allocate_error(istat, sub, 'kminor_upper')
2042 3072 : ierr = pio_inq_varid(fh, 'kminor_upper', vid)
2043 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': kminor_upper not found')
2044 3072 : ierr = pio_get_var(fh, vid, kminor_upper)
2045 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading kminor_upper')
2046 :
2047 : ! integrated Planck function by band
2048 3072 : ierr = pio_inq_varid(fh, 'totplnk', vid)
2049 3072 : if (ierr == PIO_NOERR) then
2050 6144 : allocate(totplnk(temperature_Planck,bnd), stat=istat)
2051 1536 : call handle_allocate_error(istat, sub, 'totplnk')
2052 1536 : ierr = pio_get_var(fh, vid, totplnk)
2053 1536 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading totplnk')
2054 : end if
2055 :
2056 : ! Planck fractions
2057 3072 : ierr = pio_inq_varid(fh, 'plank_fraction', vid)
2058 3072 : if (ierr == PIO_NOERR) then
2059 9216 : allocate(planck_frac(gpt,mixing_fraction,pressure_interp,temperature), stat=istat)
2060 1536 : call handle_allocate_error(istat, sub, 'planck_frac')
2061 1536 : ierr = pio_get_var(fh, vid, planck_frac)
2062 1536 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading plank_fraction')
2063 : end if
2064 :
2065 3072 : ierr = pio_inq_varid(fh, 'optimal_angle_fit', vid)
2066 3072 : if (ierr == PIO_NOERR) then
2067 6144 : allocate(optimal_angle_fit(fit_coeffs, bnd), stat=istat)
2068 1536 : call handle_allocate_error(istat, sub, 'optiman_angle_fit')
2069 1536 : ierr = pio_get_var(fh, vid, optimal_angle_fit)
2070 1536 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading optimal_angle_fit')
2071 : end if
2072 :
2073 3072 : ierr = pio_inq_varid(fh, 'solar_source_quiet', vid)
2074 3072 : if (ierr == PIO_NOERR) then
2075 4608 : allocate(solar_src_quiet(gpt), stat=istat)
2076 1536 : call handle_allocate_error(istat, sub, 'solar_src_quiet')
2077 1536 : ierr = pio_get_var(fh, vid, solar_src_quiet)
2078 1536 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading solar_source_quiet')
2079 : end if
2080 :
2081 3072 : ierr = pio_inq_varid(fh, 'solar_source_facular', vid)
2082 3072 : if (ierr == PIO_NOERR) then
2083 4608 : allocate(solar_src_facular(gpt), stat=istat)
2084 1536 : call handle_allocate_error(istat, sub, 'solar_src_facular')
2085 1536 : ierr = pio_get_var(fh, vid, solar_src_facular)
2086 1536 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading solar_source_facular')
2087 : end if
2088 :
2089 3072 : ierr = pio_inq_varid(fh, 'solar_source_sunspot', vid)
2090 3072 : if (ierr == PIO_NOERR) then
2091 4608 : allocate(solar_src_sunspot(gpt), stat=istat)
2092 1536 : call handle_allocate_error(istat, sub, 'solar_src_sunspot')
2093 1536 : ierr = pio_get_var(fh, vid, solar_src_sunspot)
2094 1536 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading solar_source_sunspot')
2095 : end if
2096 :
2097 3072 : ierr = pio_inq_varid(fh, 'tsi_default', vid)
2098 3072 : if (ierr == PIO_NOERR) then
2099 1536 : ierr = pio_get_var(fh, vid, tsi_default)
2100 1536 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading tsi_default')
2101 : end if
2102 :
2103 3072 : ierr = pio_inq_varid(fh, 'mg_default', vid)
2104 3072 : if (ierr == PIO_NOERR) then
2105 1536 : ierr = pio_get_var(fh, vid, mg_default)
2106 1536 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading mg_default')
2107 : end if
2108 :
2109 3072 : ierr = pio_inq_varid(fh, 'sb_default', vid)
2110 3072 : if (ierr == PIO_NOERR) then
2111 1536 : ierr = pio_get_var(fh, vid, sb_default)
2112 1536 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading sb_default')
2113 : end if
2114 :
2115 : ! rayleigh scattering contribution in lower part of atmosphere
2116 3072 : ierr = pio_inq_varid(fh, 'rayl_lower', vid)
2117 3072 : if (ierr == PIO_NOERR) then
2118 7680 : allocate(rayl_lower(gpt,mixing_fraction,temperature), stat=istat)
2119 1536 : call handle_allocate_error(istat, sub, 'rayl_lower')
2120 1536 : ierr = pio_get_var(fh, vid, rayl_lower)
2121 1536 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading rayl_lower')
2122 : end if
2123 :
2124 : ! rayleigh scattering contribution in upper part of atmosphere
2125 3072 : ierr = pio_inq_varid(fh, 'rayl_upper', vid)
2126 3072 : if (ierr == PIO_NOERR) then
2127 7680 : allocate(rayl_upper(gpt,mixing_fraction,temperature), stat=istat)
2128 1536 : call handle_allocate_error(istat, sub, 'rayl_upper')
2129 1536 : ierr = pio_get_var(fh, vid, rayl_upper)
2130 1536 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading rayl_upper')
2131 : end if
2132 :
2133 9216 : allocate(gas_minor(minorabsorbers), stat=istat)
2134 3072 : call handle_allocate_error(istat, sub, 'gas_minor')
2135 3072 : ierr = pio_inq_varid(fh, 'gas_minor', vid)
2136 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': gas_minor not found')
2137 3072 : ierr = pio_get_var(fh, vid, gas_minor)
2138 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading gas_minor')
2139 :
2140 9216 : allocate(identifier_minor(minorabsorbers), stat=istat)
2141 3072 : call handle_allocate_error(istat, sub, 'identifier_minor')
2142 3072 : ierr = pio_inq_varid(fh, 'identifier_minor', vid)
2143 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': identifier_minor not found')
2144 3072 : ierr = pio_get_var(fh, vid, identifier_minor)
2145 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading identifier_minor')
2146 :
2147 9216 : allocate(minor_gases_lower(minor_absorber_intervals_lower), stat=istat)
2148 3072 : call handle_allocate_error(istat, sub, 'minor_gases_lower')
2149 3072 : ierr = pio_inq_varid(fh, 'minor_gases_lower', vid)
2150 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': minor_gases_lower not found')
2151 3072 : ierr = pio_get_var(fh, vid, minor_gases_lower)
2152 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading minor_gases_lower')
2153 :
2154 9216 : allocate(minor_gases_upper(minor_absorber_intervals_upper), stat=istat)
2155 3072 : call handle_allocate_error(istat, sub, 'minor_gases_upper')
2156 3072 : ierr = pio_inq_varid(fh, 'minor_gases_upper', vid)
2157 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': minor_gases_upper not found')
2158 3072 : ierr = pio_get_var(fh, vid, minor_gases_upper)
2159 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading minor_gases_upper')
2160 :
2161 12288 : allocate(minor_limits_gpt_lower(pairs,minor_absorber_intervals_lower), stat=istat)
2162 3072 : call handle_allocate_error(istat, sub, 'minor_limits_gpt_lower')
2163 3072 : ierr = pio_inq_varid(fh, 'minor_limits_gpt_lower', vid)
2164 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': minor_limits_gpt_lower not found')
2165 3072 : ierr = pio_get_var(fh, vid, minor_limits_gpt_lower)
2166 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading minor_limits_gpt_lower')
2167 :
2168 12288 : allocate(minor_limits_gpt_upper(pairs,minor_absorber_intervals_upper), stat=istat)
2169 3072 : call handle_allocate_error(istat, sub, 'minor_limits_gpt_upper')
2170 3072 : ierr = pio_inq_varid(fh, 'minor_limits_gpt_upper', vid)
2171 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': minor_limits_gpt_upper not found')
2172 3072 : ierr = pio_get_var(fh, vid, minor_limits_gpt_upper)
2173 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading minor_limits_gpt_upper')
2174 :
2175 : ! Read as integer and convert to logical
2176 9216 : allocate(int2log(minor_absorber_intervals_lower), stat=istat)
2177 3072 : call handle_allocate_error(istat, sub, 'int2log for lower')
2178 :
2179 9216 : allocate(minor_scales_with_density_lower(minor_absorber_intervals_lower), stat=istat)
2180 3072 : call handle_allocate_error(istat, sub, 'minor_scales_with_density_lower')
2181 3072 : ierr = pio_inq_varid(fh, 'minor_scales_with_density_lower', vid)
2182 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': minor_scales_with_density_lower not found')
2183 3072 : ierr = pio_get_var(fh, vid, int2log)
2184 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading minor_scales_with_density_lower')
2185 147456 : do i = 1,minor_absorber_intervals_lower
2186 147456 : if (int2log(i) .eq. 0) then
2187 53760 : minor_scales_with_density_lower(i) = .false.
2188 : else
2189 90624 : minor_scales_with_density_lower(i) = .true.
2190 : end if
2191 : end do
2192 :
2193 9216 : allocate(scale_by_complement_lower(minor_absorber_intervals_lower), stat=istat)
2194 3072 : call handle_allocate_error(istat, sub, 'scale_by_complement_lower')
2195 3072 : ierr = pio_inq_varid(fh, 'scale_by_complement_lower', vid)
2196 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': scale_by_complement_lower not found')
2197 3072 : ierr = pio_get_var(fh, vid, int2log)
2198 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading scale_by_complement_lower')
2199 147456 : do i = 1,minor_absorber_intervals_lower
2200 147456 : if (int2log(i) .eq. 0) then
2201 104448 : scale_by_complement_lower(i) = .false.
2202 : else
2203 39936 : scale_by_complement_lower(i) = .true.
2204 : end if
2205 : end do
2206 :
2207 3072 : deallocate(int2log)
2208 :
2209 : ! Read as integer and convert to logical
2210 9216 : allocate(int2log(minor_absorber_intervals_upper), stat=istat)
2211 3072 : call handle_allocate_error(istat, sub, 'int2log for upper')
2212 :
2213 9216 : allocate(minor_scales_with_density_upper(minor_absorber_intervals_upper), stat=istat)
2214 3072 : call handle_allocate_error(istat, sub, 'minor_scales_with_density_upper')
2215 3072 : ierr = pio_inq_varid(fh, 'minor_scales_with_density_upper', vid)
2216 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': minor_scales_with_density_upper not found')
2217 3072 : ierr = pio_get_var(fh, vid, int2log)
2218 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading minor_scales_with_density_upper')
2219 92160 : do i = 1,minor_absorber_intervals_upper
2220 92160 : if (int2log(i) .eq. 0) then
2221 46080 : minor_scales_with_density_upper(i) = .false.
2222 : else
2223 43008 : minor_scales_with_density_upper(i) = .true.
2224 : end if
2225 : end do
2226 :
2227 9216 : allocate(scale_by_complement_upper(minor_absorber_intervals_upper), stat=istat)
2228 3072 : call handle_allocate_error(istat, sub, 'scale_by_complement_upper')
2229 3072 : ierr = pio_inq_varid(fh, 'scale_by_complement_upper', vid)
2230 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': scale_by_complement_upper not found')
2231 3072 : ierr = pio_get_var(fh, vid, int2log)
2232 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading scale_by_complement_upper')
2233 92160 : do i = 1,minor_absorber_intervals_upper
2234 92160 : if (int2log(i) .eq. 0) then
2235 70656 : scale_by_complement_upper(i) = .false.
2236 : else
2237 18432 : scale_by_complement_upper(i) = .true.
2238 : end if
2239 : end do
2240 :
2241 3072 : deallocate(int2log)
2242 :
2243 9216 : allocate(scaling_gas_lower(minor_absorber_intervals_lower), stat=istat)
2244 3072 : call handle_allocate_error(istat, sub, 'scaling_gas_lower')
2245 3072 : ierr = pio_inq_varid(fh, 'scaling_gas_lower', vid)
2246 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': scaling_gas_lower not found')
2247 3072 : ierr = pio_get_var(fh, vid, scaling_gas_lower)
2248 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading scaling_gas_lower')
2249 :
2250 9216 : allocate(scaling_gas_upper(minor_absorber_intervals_upper), stat=istat)
2251 3072 : call handle_allocate_error(istat, sub, 'scaling_gas_upper')
2252 3072 : ierr = pio_inq_varid(fh, 'scaling_gas_upper', vid)
2253 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': scaling_gas_upper not found')
2254 3072 : ierr = pio_get_var(fh, vid, scaling_gas_upper)
2255 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading scaling_gas_upper')
2256 :
2257 9216 : allocate(kminor_start_lower(minor_absorber_intervals_lower), stat=istat)
2258 3072 : call handle_allocate_error(istat, sub, 'kminor_start_lower')
2259 3072 : ierr = pio_inq_varid(fh, 'kminor_start_lower', vid)
2260 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': kminor_start_lower not found')
2261 3072 : ierr = pio_get_var(fh, vid, kminor_start_lower)
2262 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading kminor_start_lower')
2263 :
2264 9216 : allocate(kminor_start_upper(minor_absorber_intervals_upper), stat=istat)
2265 3072 : call handle_allocate_error(istat, sub, 'kminor_start_upper')
2266 3072 : ierr = pio_inq_varid(fh, 'kminor_start_upper', vid)
2267 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': kminor_start_upper not found')
2268 3072 : ierr = pio_get_var(fh, vid, kminor_start_upper)
2269 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading kminor_start_upper')
2270 :
2271 : ! Close file
2272 3072 : call pio_closefile(fh)
2273 :
2274 : ! Initialize the gas optics object with data. The calls are slightly different depending
2275 : ! on whether the radiation sources are internal to the atmosphere (longwave) or external (shortwave)
2276 :
2277 3072 : if (allocated(totplnk) .and. allocated(planck_frac)) then
2278 : error_msg = kdist%load( &
2279 : available_gases, gas_names, key_species, &
2280 : band2gpt, band_lims_wavenum, &
2281 : press_ref, press_ref_trop, temp_ref, &
2282 : temp_ref_p, temp_ref_t, vmr_ref, &
2283 : kmajor, kminor_lower, kminor_upper, &
2284 : gas_minor, identifier_minor, &
2285 : minor_gases_lower, minor_gases_upper, &
2286 : minor_limits_gpt_lower, minor_limits_gpt_upper, &
2287 : minor_scales_with_density_lower, &
2288 : minor_scales_with_density_upper, &
2289 : scaling_gas_lower, scaling_gas_upper, &
2290 : scale_by_complement_lower, scale_by_complement_upper, &
2291 : kminor_start_lower, kminor_start_upper, &
2292 : totplnk, planck_frac, rayl_lower, rayl_upper, &
2293 1536 : optimal_angle_fit)
2294 1536 : else if (allocated(solar_src_quiet)) then
2295 : error_msg = kdist%load( &
2296 : available_gases, gas_names, key_species, &
2297 : band2gpt, band_lims_wavenum, &
2298 : press_ref, press_ref_trop, temp_ref, &
2299 : temp_ref_p, temp_ref_t, vmr_ref, &
2300 : kmajor, kminor_lower, kminor_upper, &
2301 : gas_minor, identifier_minor, &
2302 : minor_gases_lower, minor_gases_upper, &
2303 : minor_limits_gpt_lower, minor_limits_gpt_upper, &
2304 : minor_scales_with_density_lower, &
2305 : minor_scales_with_density_upper, &
2306 : scaling_gas_lower, scaling_gas_upper, &
2307 : scale_by_complement_lower, &
2308 : scale_by_complement_upper, &
2309 : kminor_start_lower, &
2310 : kminor_start_upper, &
2311 : solar_src_quiet, solar_src_facular, solar_src_sunspot, &
2312 : tsi_default, mg_default, sb_default, &
2313 1536 : rayl_lower, rayl_upper)
2314 : else
2315 0 : error_msg = 'must supply either totplnk and planck_frac, or solar_src_[*]'
2316 : end if
2317 :
2318 3072 : call stop_on_err(error_msg, sub, 'kdist%load')
2319 :
2320 0 : deallocate( &
2321 0 : gas_names, key_species, &
2322 0 : band2gpt, band_lims_wavenum, &
2323 0 : press_ref, temp_ref, vmr_ref, &
2324 0 : kmajor, kminor_lower, kminor_upper, &
2325 0 : gas_minor, identifier_minor, &
2326 0 : minor_gases_lower, minor_gases_upper, &
2327 0 : minor_limits_gpt_lower, &
2328 0 : minor_limits_gpt_upper, &
2329 0 : minor_scales_with_density_lower, &
2330 0 : minor_scales_with_density_upper, &
2331 0 : scale_by_complement_lower, &
2332 0 : scale_by_complement_upper, &
2333 0 : scaling_gas_lower, scaling_gas_upper, &
2334 3072 : kminor_start_lower, kminor_start_upper)
2335 :
2336 3072 : if (allocated(totplnk)) deallocate(totplnk)
2337 3072 : if (allocated(planck_frac)) deallocate(planck_frac)
2338 3072 : if (allocated(optimal_angle_fit)) deallocate(optimal_angle_fit)
2339 3072 : if (allocated(solar_src_quiet)) deallocate(solar_src_quiet)
2340 3072 : if (allocated(solar_src_facular)) deallocate(solar_src_facular)
2341 3072 : if (allocated(solar_src_sunspot)) deallocate(solar_src_sunspot)
2342 3072 : if (allocated(rayl_lower)) deallocate(rayl_lower)
2343 3072 : if (allocated(rayl_upper)) deallocate(rayl_upper)
2344 :
2345 6144 : end subroutine coefs_init
2346 :
2347 : !=========================================================================================
2348 :
2349 247680 : subroutine initialize_rrtmgp_fluxes(ncol, nlevels, nbands, fluxes, do_direct)
2350 :
2351 : ! Allocate flux arrays and set values to zero.
2352 :
2353 : ! Arguments
2354 : integer, intent(in) :: ncol, nlevels, nbands
2355 : class(ty_fluxes_broadband), intent(inout) :: fluxes
2356 : logical, optional, intent(in) :: do_direct
2357 :
2358 : ! Local variables
2359 : logical :: do_direct_local
2360 : integer :: istat
2361 : character(len=*), parameter :: sub = 'initialize_rrtmgp_fluxes'
2362 : !----------------------------------------------------------------------------
2363 :
2364 247680 : if (present(do_direct)) then
2365 : do_direct_local = .true.
2366 : else
2367 123840 : do_direct_local = .false.
2368 : end if
2369 :
2370 : ! Broadband fluxes
2371 931464 : allocate(fluxes%flux_up(ncol, nlevels), stat=istat)
2372 247680 : call handle_allocate_error(istat, sub, 'fluxes%flux_up')
2373 683784 : allocate(fluxes%flux_dn(ncol, nlevels), stat=istat)
2374 247680 : call handle_allocate_error(istat, sub, 'fluxes%flux_dn')
2375 683784 : allocate(fluxes%flux_net(ncol, nlevels), stat=istat)
2376 247680 : call handle_allocate_error(istat, sub, 'fluxes%flux_net')
2377 247680 : if (do_direct_local) then
2378 312264 : allocate(fluxes%flux_dn_dir(ncol, nlevels), stat=istat)
2379 123840 : call handle_allocate_error(istat, sub, 'fluxes%flux_dn_dir')
2380 : end if
2381 :
2382 : select type (fluxes)
2383 : type is (ty_fluxes_byband)
2384 : ! Fluxes by band always needed for SW. Only allocate for LW
2385 : ! when spectralflux is true.
2386 123840 : if (nbands == nswbands .or. spectralflux) then
2387 279972 : allocate(fluxes%bnd_flux_up(ncol, nlevels, nbands), stat=istat)
2388 61920 : call handle_allocate_error(istat, sub, 'fluxes%bnd_flux_up')
2389 218052 : allocate(fluxes%bnd_flux_dn(ncol, nlevels, nbands), stat=istat)
2390 61920 : call handle_allocate_error(istat, sub, 'fluxes%bnd_flux_dn')
2391 218052 : allocate(fluxes%bnd_flux_net(ncol, nlevels, nbands), stat=istat)
2392 61920 : call handle_allocate_error(istat, sub, 'fluxes%bnd_flux_net')
2393 61920 : if (do_direct_local) then
2394 218052 : allocate(fluxes%bnd_flux_dn_dir(ncol, nlevels, nbands), stat=istat)
2395 61920 : call handle_allocate_error(istat, sub, 'fluxes%bnd_flux_dn_dir')
2396 : end if
2397 : end if
2398 : end select
2399 :
2400 : ! Initialize
2401 247680 : call reset_fluxes(fluxes)
2402 :
2403 247680 : end subroutine initialize_rrtmgp_fluxes
2404 :
2405 : !=========================================================================================
2406 :
2407 247680 : subroutine reset_fluxes(fluxes)
2408 :
2409 : ! Reset flux arrays to zero.
2410 :
2411 : class(ty_fluxes_broadband), intent(inout) :: fluxes
2412 : !----------------------------------------------------------------------------
2413 :
2414 : ! Reset broadband fluxes
2415 300797280 : fluxes%flux_up(:,:) = 0._r8
2416 300797280 : fluxes%flux_dn(:,:) = 0._r8
2417 300797280 : fluxes%flux_net(:,:) = 0._r8
2418 104352480 : if (associated(fluxes%flux_dn_dir)) fluxes%flux_dn_dir(:,:) = 0._r8
2419 :
2420 : select type (fluxes)
2421 : type is (ty_fluxes_byband)
2422 : ! Reset band-by-band fluxes
2423 729662400 : if (associated(fluxes%bnd_flux_up)) fluxes%bnd_flux_up(:,:,:) = 0._r8
2424 729724320 : if (associated(fluxes%bnd_flux_dn)) fluxes%bnd_flux_dn(:,:,:) = 0._r8
2425 729724320 : if (associated(fluxes%bnd_flux_net)) fluxes%bnd_flux_net(:,:,:) = 0._r8
2426 729848160 : if (associated(fluxes%bnd_flux_dn_dir)) fluxes%bnd_flux_dn_dir(:,:,:) = 0._r8
2427 : end select
2428 :
2429 247680 : end subroutine reset_fluxes
2430 :
2431 : !=========================================================================================
2432 :
2433 185760 : subroutine free_optics_sw(optics)
2434 :
2435 : type(ty_optical_props_2str), intent(inout) :: optics
2436 :
2437 185760 : if (allocated(optics%tau)) deallocate(optics%tau)
2438 185760 : if (allocated(optics%ssa)) deallocate(optics%ssa)
2439 185760 : if (allocated(optics%g)) deallocate(optics%g)
2440 185760 : call optics%finalize()
2441 :
2442 185760 : end subroutine free_optics_sw
2443 :
2444 : !=========================================================================================
2445 :
2446 123840 : subroutine free_optics_lw(optics)
2447 :
2448 : type(ty_optical_props_1scl), intent(inout) :: optics
2449 :
2450 123840 : if (allocated(optics%tau)) deallocate(optics%tau)
2451 123840 : call optics%finalize()
2452 :
2453 123840 : end subroutine free_optics_lw
2454 :
2455 : !=========================================================================================
2456 :
2457 247680 : subroutine free_fluxes(fluxes)
2458 :
2459 : class(ty_fluxes_broadband), intent(inout) :: fluxes
2460 :
2461 247680 : if (associated(fluxes%flux_up)) deallocate(fluxes%flux_up)
2462 247680 : if (associated(fluxes%flux_dn)) deallocate(fluxes%flux_dn)
2463 247680 : if (associated(fluxes%flux_net)) deallocate(fluxes%flux_net)
2464 247680 : if (associated(fluxes%flux_dn_dir)) deallocate(fluxes%flux_dn_dir)
2465 :
2466 : select type (fluxes)
2467 : type is (ty_fluxes_byband)
2468 61920 : if (associated(fluxes%bnd_flux_up)) deallocate(fluxes%bnd_flux_up)
2469 123840 : if (associated(fluxes%bnd_flux_dn)) deallocate(fluxes%bnd_flux_dn)
2470 123840 : if (associated(fluxes%bnd_flux_net)) deallocate(fluxes%bnd_flux_net)
2471 247680 : if (associated(fluxes%bnd_flux_dn_dir)) deallocate(fluxes%bnd_flux_dn_dir)
2472 : end select
2473 :
2474 247680 : end subroutine free_fluxes
2475 :
2476 : !=========================================================================================
2477 :
2478 30960 : subroutine modified_cloud_fraction(ncol, cld, cldfsnow, cldfgrau, cldfprime)
2479 :
2480 : ! Compute modified cloud fraction, cldfprime.
2481 : ! 1. initialize as cld
2482 : ! 2. modify for snow if cldfsnow is available. use max(cld, cldfsnow)
2483 : ! 3. modify for graupel if cldfgrau is available and graupel_in_rad is true.
2484 : ! use max(cldfprime, cldfgrau)
2485 :
2486 : ! Arguments
2487 : integer, intent(in) :: ncol
2488 : real(r8), pointer :: cld(:,:) ! cloud fraction
2489 : real(r8), pointer :: cldfsnow(:,:) ! cloud fraction of just "snow clouds"
2490 : real(r8), pointer :: cldfgrau(:,:) ! cloud fraction of just "graupel clouds"
2491 : real(r8), intent(out) :: cldfprime(:,:) ! modified cloud fraction
2492 :
2493 : ! Local variables
2494 : integer :: i, k
2495 : !----------------------------------------------------------------------------
2496 :
2497 30960 : if (associated(cldfsnow)) then
2498 2910240 : do k = 1, pver
2499 48108240 : do i = 1, ncol
2500 48077280 : cldfprime(i,k) = max(cld(i,k), cldfsnow(i,k))
2501 : end do
2502 : end do
2503 : else
2504 0 : cldfprime(:ncol,:) = cld(:ncol,:)
2505 : end if
2506 :
2507 30960 : if (associated(cldfgrau) .and. graupel_in_rad) then
2508 0 : do k = 1, pver
2509 0 : do i = 1, ncol
2510 0 : cldfprime(i,k) = max(cldfprime(i,k), cldfgrau(i,k))
2511 : end do
2512 : end do
2513 : end if
2514 :
2515 30960 : end subroutine modified_cloud_fraction
2516 :
2517 : !=========================================================================================
2518 :
2519 412456 : subroutine stop_on_err(errmsg, sub, info)
2520 :
2521 : ! call endrun if RRTMGP function returns non-empty error message.
2522 :
2523 : character(len=*), intent(in) :: errmsg ! return message from RRTMGP function
2524 : character(len=*), intent(in) :: sub ! name of calling subroutine
2525 : character(len=*), intent(in) :: info ! name of called function
2526 :
2527 412456 : if (len_trim(errmsg) > 0) then
2528 0 : call endrun(trim(sub)//': ERROR: '//trim(info)//': '//trim(errmsg))
2529 : end if
2530 :
2531 412456 : end subroutine stop_on_err
2532 :
2533 : !=========================================================================================
2534 :
2535 743040 : end module radiation
2536 :
|
