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 : use solar_irrad_data, only: sol_tsi
20 :
21 : use time_manager, only: get_nstep, is_first_step, is_first_restart_step, &
22 : get_curr_calday, get_step_size
23 :
24 : use rad_constituents, only: N_DIAG, rad_cnst_get_call_list, rad_cnst_get_gas, rad_cnst_out
25 :
26 : use rrtmgp_inputs, only: rrtmgp_inputs_init
27 :
28 : use radconstants, only: nradgas, gasnamelength, gaslist, nswbands, nlwbands, &
29 : nswgpts, set_wavenumber_bands
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 5232240 : 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 5232240 : if (present(timestep)) then
359 5232240 : nstep = timestep
360 : else
361 0 : nstep = get_nstep()
362 : end if
363 :
364 3739968 : 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 3739968 : .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 1492272 : .or. nstep <= irad_always
373 : case default
374 5232240 : call endrun('radiation_do: unknown operation:'//op)
375 : end select
376 :
377 5232240 : end function radiation_do
378 :
379 : !================================================================================================
380 :
381 1492272 : 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 1492272 : radiation_nextsw_cday = -1._r8
396 1492272 : dosw = .false.
397 1492272 : nstep = get_nstep()
398 1492272 : dtime = get_step_size()
399 1492272 : offset = 0
400 : do while (.not. dosw)
401 2247696 : nstep = nstep + 1
402 2247696 : offset = offset + dtime
403 2247696 : if (radiation_do('sw', nstep)) then
404 1492272 : radiation_nextsw_cday = get_curr_calday(offset=offset)
405 : dosw = .true.
406 : end if
407 : end do
408 1492272 : 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 1492272 : if (get_nstep() >= 1) then
414 1486080 : caldayp1 = get_curr_calday(offset=int(dtime))
415 1486080 : if (caldayp1 /= radiation_nextsw_cday) radiation_nextsw_cday = -1._r8
416 : end if
417 :
418 1492272 : 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 :
499 : ! The spectral band boundaries need to be set before this init is called.
500 1536 : call rrtmgp_inputs_init(ktopcam, ktoprad)
501 :
502 : ! initialize output fields for offline driver
503 1536 : call rad_data_init(pbuf2d)
504 :
505 1536 : call cloud_rad_props_init()
506 :
507 1536 : cld_idx = pbuf_get_index('CLD')
508 1536 : cldfsnow_idx = pbuf_get_index('CLDFSNOW', errcode=ierr)
509 1536 : cldfgrau_idx = pbuf_get_index('CLDFGRAU', errcode=ierr)
510 :
511 1536 : if (is_first_step()) then
512 768 : call pbuf_set_field(pbuf2d, qrl_idx, 0._r8)
513 : end if
514 :
515 : ! Set the radiation timestep for cosz calculations if requested using
516 : ! the adjusted iradsw value from radiation
517 1536 : if (use_rad_dt_cosz) then
518 0 : dtime = get_step_size()
519 0 : dt_avg = iradsw*dtime
520 : end if
521 :
522 : ! Surface components to get radiation computed today
523 1536 : if (.not. is_first_restart_step()) then
524 768 : nextsw_cday = get_curr_calday()
525 : end if
526 :
527 : call phys_getopts(history_amwg_out = history_amwg, &
528 : history_vdiag_out = history_vdiag, &
529 : history_budget_out = history_budget, &
530 1536 : history_budget_histfile_num_out = history_budget_histfile_num)
531 :
532 : ! "irad_always" is number of time steps to execute radiation continuously from
533 : ! start of initial OR restart run
534 1536 : nstep = get_nstep()
535 1536 : if (irad_always > 0) then
536 0 : irad_always = irad_always + nstep
537 : end if
538 :
539 1536 : if (docosp) call cospsimulator_intr_init()
540 :
541 4608 : allocate(cosp_cnt(begchunk:endchunk), stat=istat)
542 1536 : call handle_allocate_error(istat, sub, 'cosp_cnt')
543 1536 : if (is_first_restart_step()) then
544 3864 : cosp_cnt(begchunk:endchunk) = cosp_cnt_init
545 : else
546 3864 : cosp_cnt(begchunk:endchunk) = 0
547 : end if
548 :
549 : ! Add fields to history buffer
550 :
551 : call addfld('TOT_CLD_VISTAU', (/ 'lev' /), 'A', '1', &
552 : 'Total gbx cloud extinction visible sw optical depth', &
553 3072 : sampling_seq='rad_lwsw', flag_xyfill=.true.)
554 : call addfld('TOT_ICLD_VISTAU', (/ 'lev' /), 'A', '1', &
555 : 'Total in-cloud extinction visible sw optical depth', &
556 3072 : sampling_seq='rad_lwsw', flag_xyfill=.true.)
557 : call addfld('LIQ_ICLD_VISTAU', (/ 'lev' /), 'A', '1', &
558 : 'Liquid in-cloud extinction visible sw optical depth', &
559 3072 : sampling_seq='rad_lwsw', flag_xyfill=.true.)
560 : call addfld('ICE_ICLD_VISTAU', (/ 'lev' /), 'A', '1', &
561 : 'Ice in-cloud extinction visible sw optical depth', &
562 3072 : sampling_seq='rad_lwsw', flag_xyfill=.true.)
563 :
564 1536 : if (cldfsnow_idx > 0) then
565 : call addfld('SNOW_ICLD_VISTAU', (/ 'lev' /), 'A', '1', &
566 : 'Snow in-cloud extinction visible sw optical depth', &
567 3072 : sampling_seq='rad_lwsw', flag_xyfill=.true.)
568 : end if
569 1536 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
570 : call addfld('GRAU_ICLD_VISTAU', (/ 'lev' /), 'A', '1', &
571 : 'Graupel in-cloud extinction visible sw optical depth', &
572 0 : sampling_seq='rad_lwsw', flag_xyfill=.true.)
573 : endif
574 :
575 : ! get list of active radiation calls
576 1536 : call rad_cnst_get_call_list(active_calls)
577 :
578 : ! Add shortwave radiation fields to history master field list.
579 :
580 18432 : do icall = 0, N_DIAG
581 :
582 18432 : if (active_calls(icall)) then
583 :
584 : call addfld('SOLIN'//diag(icall), horiz_only, 'A', 'W/m2', &
585 1536 : 'Solar insolation', sampling_seq='rad_lwsw')
586 : call addfld('QRS'//diag(icall), (/ 'lev' /), 'A', 'K/s', &
587 3072 : 'Solar heating rate', sampling_seq='rad_lwsw')
588 : call addfld('QRSC'//diag(icall), (/ 'lev' /), 'A', 'K/s', &
589 3072 : 'Clearsky solar heating rate', sampling_seq='rad_lwsw')
590 : call addfld('FSNT'//diag(icall), horiz_only, 'A', 'W/m2', &
591 1536 : 'Net solar flux at top of model', sampling_seq='rad_lwsw')
592 : call addfld('FSNTC'//diag(icall), horiz_only, 'A', 'W/m2', &
593 1536 : 'Clearsky net solar flux at top of model', sampling_seq='rad_lwsw')
594 : call addfld('FSNTOA'//diag(icall), horiz_only, 'A', 'W/m2', &
595 1536 : 'Net solar flux at top of atmosphere', sampling_seq='rad_lwsw')
596 : call addfld('FSNTOAC'//diag(icall), horiz_only, 'A', 'W/m2', &
597 1536 : 'Clearsky net solar flux at top of atmosphere', sampling_seq='rad_lwsw')
598 : call addfld('SWCF'//diag(icall), horiz_only, 'A', 'W/m2', &
599 1536 : 'Shortwave cloud forcing', sampling_seq='rad_lwsw')
600 : call addfld('FSUTOA'//diag(icall), horiz_only, 'A', 'W/m2', &
601 1536 : 'Upwelling solar flux at top of atmosphere', sampling_seq='rad_lwsw')
602 : call addfld('FSNIRTOA'//diag(icall), horiz_only, 'A', 'W/m2', &
603 1536 : 'Net near-infrared flux (Nimbus-7 WFOV) at top of atmosphere', sampling_seq='rad_lwsw')
604 : call addfld('FSNRTOAC'//diag(icall), horiz_only, 'A', 'W/m2', &
605 1536 : 'Clearsky net near-infrared flux (Nimbus-7 WFOV) at top of atmosphere', sampling_seq='rad_lwsw')
606 : call addfld('FSNRTOAS'//diag(icall), horiz_only, 'A', 'W/m2', &
607 1536 : 'Net near-infrared flux (>= 0.7 microns) at top of atmosphere', sampling_seq='rad_lwsw')
608 : call addfld('FSN200'//diag(icall), horiz_only, 'A', 'W/m2', &
609 1536 : 'Net shortwave flux at 200 mb', sampling_seq='rad_lwsw')
610 : call addfld('FSN200C'//diag(icall), horiz_only, 'A', 'W/m2', &
611 1536 : 'Clearsky net shortwave flux at 200 mb', sampling_seq='rad_lwsw')
612 : call addfld('FSNR'//diag(icall), horiz_only, 'A', 'W/m2', &
613 1536 : 'Net solar flux at tropopause', sampling_seq='rad_lwsw')
614 : call addfld('SOLL'//diag(icall), horiz_only, 'A', 'W/m2', &
615 1536 : 'Solar downward near infrared direct to surface', sampling_seq='rad_lwsw')
616 : call addfld('SOLS'//diag(icall), horiz_only, 'A', 'W/m2', &
617 1536 : 'Solar downward visible direct to surface', sampling_seq='rad_lwsw')
618 : call addfld('SOLLD'//diag(icall), horiz_only, 'A', 'W/m2', &
619 1536 : 'Solar downward near infrared diffuse to surface', sampling_seq='rad_lwsw')
620 : call addfld('SOLSD'//diag(icall), horiz_only, 'A', 'W/m2', &
621 1536 : 'Solar downward visible diffuse to surface', sampling_seq='rad_lwsw')
622 : call addfld('FSNS'//diag(icall), horiz_only, 'A', 'W/m2', &
623 1536 : 'Net solar flux at surface', sampling_seq='rad_lwsw')
624 : call addfld('FSNSC'//diag(icall), horiz_only, 'A', 'W/m2', &
625 1536 : 'Clearsky net solar flux at surface', sampling_seq='rad_lwsw')
626 : call addfld('FSDS'//diag(icall), horiz_only, 'A', 'W/m2', &
627 1536 : 'Downwelling solar flux at surface', sampling_seq='rad_lwsw')
628 : call addfld('FSDSC'//diag(icall), horiz_only, 'A', 'W/m2', &
629 1536 : 'Clearsky downwelling solar flux at surface', sampling_seq='rad_lwsw')
630 :
631 : ! Fluxes on CAM grid
632 : call addfld('FUS'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', &
633 3072 : 'Shortwave upward flux', sampling_seq='rad_lwsw')
634 : call addfld('FDS'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', &
635 3072 : 'Shortwave downward flux', sampling_seq='rad_lwsw')
636 : call addfld('FUSC'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', &
637 3072 : 'Shortwave clear-sky upward flux', sampling_seq='rad_lwsw')
638 : call addfld('FDSC'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', &
639 3072 : 'Shortwave clear-sky downward flux', sampling_seq='rad_lwsw')
640 :
641 1536 : if (history_amwg) then
642 1536 : call add_default('SOLIN'//diag(icall), 1, ' ')
643 1536 : call add_default('QRS'//diag(icall), 1, ' ')
644 1536 : call add_default('FSNT'//diag(icall), 1, ' ')
645 1536 : call add_default('FSNTC'//diag(icall), 1, ' ')
646 1536 : call add_default('FSNTOA'//diag(icall), 1, ' ')
647 1536 : call add_default('FSNTOAC'//diag(icall), 1, ' ')
648 1536 : call add_default('SWCF'//diag(icall), 1, ' ')
649 1536 : call add_default('FSNS'//diag(icall), 1, ' ')
650 1536 : call add_default('FSNSC'//diag(icall), 1, ' ')
651 1536 : call add_default('FSUTOA'//diag(icall), 1, ' ')
652 1536 : call add_default('FSDSC'//diag(icall), 1, ' ')
653 1536 : call add_default('FSDS'//diag(icall), 1, ' ')
654 : endif
655 :
656 : end if
657 : end do
658 :
659 1536 : if (scm_crm_mode) then
660 0 : call add_default('FUS ', 1, ' ')
661 0 : call add_default('FUSC ', 1, ' ')
662 0 : call add_default('FDS ', 1, ' ')
663 0 : call add_default('FDSC ', 1, ' ')
664 : endif
665 :
666 : ! Add longwave radiation fields to history master field list.
667 :
668 18432 : do icall = 0, N_DIAG
669 :
670 18432 : if (active_calls(icall)) then
671 : call addfld('QRL'//diag(icall), (/ 'lev' /), 'A', 'K/s', &
672 3072 : 'Longwave heating rate', sampling_seq='rad_lwsw')
673 : call addfld('QRLC'//diag(icall), (/ 'lev' /), 'A', 'K/s', &
674 3072 : 'Clearsky longwave heating rate', sampling_seq='rad_lwsw')
675 : call addfld('FLNT'//diag(icall), horiz_only, 'A', 'W/m2', &
676 1536 : 'Net longwave flux at top of model', sampling_seq='rad_lwsw')
677 : call addfld('FLNTC'//diag(icall), horiz_only, 'A', 'W/m2', &
678 1536 : 'Clearsky net longwave flux at top of model', sampling_seq='rad_lwsw')
679 : call addfld('FLUT'//diag(icall), horiz_only, 'A', 'W/m2', &
680 1536 : 'Upwelling longwave flux at top of model', sampling_seq='rad_lwsw')
681 : call addfld('FLUTC'//diag(icall), horiz_only, 'A', 'W/m2', &
682 1536 : 'Clearsky upwelling longwave flux at top of model', sampling_seq='rad_lwsw')
683 : call addfld('LWCF'//diag(icall), horiz_only, 'A', 'W/m2', &
684 1536 : 'Longwave cloud forcing', sampling_seq='rad_lwsw')
685 : call addfld('FLN200'//diag(icall), horiz_only, 'A', 'W/m2', &
686 1536 : 'Net longwave flux at 200 mb', sampling_seq='rad_lwsw')
687 : call addfld('FLN200C'//diag(icall), horiz_only, 'A', 'W/m2', &
688 1536 : 'Clearsky net longwave flux at 200 mb', sampling_seq='rad_lwsw')
689 : call addfld('FLNR'//diag(icall), horiz_only, 'A', 'W/m2', &
690 1536 : 'Net longwave flux at tropopause', sampling_seq='rad_lwsw')
691 : call addfld('FLNS'//diag(icall), horiz_only, 'A', 'W/m2', &
692 1536 : 'Net longwave flux at surface', sampling_seq='rad_lwsw')
693 : call addfld('FLNSC'//diag(icall), horiz_only, 'A', 'W/m2', &
694 1536 : 'Clearsky net longwave flux at surface', sampling_seq='rad_lwsw')
695 : call addfld('FLDS'//diag(icall), horiz_only, 'A', 'W/m2', &
696 1536 : 'Downwelling longwave flux at surface', sampling_seq='rad_lwsw')
697 : call addfld('FLDSC'//diag(icall), horiz_only, 'A', 'W/m2', &
698 1536 : 'Clearsky Downwelling longwave flux at surface', sampling_seq='rad_lwsw')
699 :
700 : ! Fluxes on CAM grid
701 : call addfld('FUL'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', &
702 3072 : 'Longwave upward flux', sampling_seq='rad_lwsw')
703 : call addfld('FDL'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', &
704 3072 : 'Longwave downward flux', sampling_seq='rad_lwsw')
705 : call addfld('FULC'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', &
706 3072 : 'Longwave clear-sky upward flux', sampling_seq='rad_lwsw')
707 : call addfld('FDLC'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', &
708 3072 : 'Longwave clear-sky downward flux', sampling_seq='rad_lwsw')
709 :
710 1536 : if (history_amwg) then
711 1536 : call add_default('QRL'//diag(icall), 1, ' ')
712 1536 : call add_default('FLNT'//diag(icall), 1, ' ')
713 1536 : call add_default('FLNTC'//diag(icall), 1, ' ')
714 1536 : call add_default('FLUT'//diag(icall), 1, ' ')
715 1536 : call add_default('FLUTC'//diag(icall), 1, ' ')
716 1536 : call add_default('LWCF'//diag(icall), 1, ' ')
717 1536 : call add_default('FLNS'//diag(icall), 1, ' ')
718 1536 : call add_default('FLNSC'//diag(icall), 1, ' ')
719 1536 : call add_default('FLDS'//diag(icall), 1, ' ')
720 : end if
721 :
722 : end if
723 : end do
724 :
725 3072 : call addfld('EMIS', (/ 'lev' /), 'A', '1', 'Cloud longwave emissivity')
726 :
727 1536 : if (scm_crm_mode) then
728 0 : call add_default ('FUL ', 1, ' ')
729 0 : call add_default ('FULC ', 1, ' ')
730 0 : call add_default ('FDL ', 1, ' ')
731 0 : call add_default ('FDLC ', 1, ' ')
732 : endif
733 :
734 : ! Heating rate needed for d(theta)/dt computation
735 3072 : call addfld ('HR',(/ 'lev' /), 'A','K/s','Heating rate needed for d(theta)/dt computation')
736 :
737 1536 : if ( history_budget .and. history_budget_histfile_num > 1 ) then
738 0 : call add_default ('QRL ', history_budget_histfile_num, ' ')
739 0 : call add_default ('QRS ', history_budget_histfile_num, ' ')
740 : end if
741 :
742 1536 : if (history_vdiag) then
743 0 : call add_default('FLUT', 2, ' ')
744 0 : call add_default('FLUT', 3, ' ')
745 : end if
746 :
747 3072 : end subroutine radiation_init
748 :
749 : !===============================================================================
750 :
751 1536 : subroutine radiation_define_restart(file)
752 :
753 : ! define variables to be written to restart file
754 :
755 : ! arguments
756 : type(file_desc_t), intent(inout) :: file
757 :
758 : ! local variables
759 : integer :: ierr
760 : !----------------------------------------------------------------------------
761 :
762 1536 : call pio_seterrorhandling(file, PIO_BCAST_ERROR)
763 :
764 1536 : ierr = pio_def_var(file, 'nextsw_cday', pio_double, nextsw_cday_desc)
765 1536 : ierr = pio_put_att(file, nextsw_cday_desc, 'long_name', 'future radiation calday for surface models')
766 1536 : if (docosp) then
767 0 : ierr = pio_def_var(File, 'cosp_cnt_init', pio_int, cospcnt_desc)
768 : end if
769 :
770 1536 : end subroutine radiation_define_restart
771 :
772 : !===============================================================================
773 :
774 1536 : subroutine radiation_write_restart(file)
775 :
776 : ! write variables to restart file
777 :
778 : ! arguments
779 : type(file_desc_t), intent(inout) :: file
780 :
781 : ! local variables
782 : integer :: ierr
783 : !----------------------------------------------------------------------------
784 3072 : ierr = pio_put_var(File, nextsw_cday_desc, (/ nextsw_cday /))
785 1536 : if (docosp) then
786 0 : ierr = pio_put_var(File, cospcnt_desc, (/cosp_cnt(begchunk)/))
787 : end if
788 :
789 1536 : end subroutine radiation_write_restart
790 :
791 : !===============================================================================
792 :
793 768 : subroutine radiation_read_restart(file)
794 :
795 : ! read variables from restart file
796 :
797 : ! arguments
798 : type(file_desc_t), intent(inout) :: file
799 :
800 : ! local variables
801 : integer :: ierr
802 : type(var_desc_t) :: vardesc
803 : integer :: err_handling
804 :
805 : !----------------------------------------------------------------------------
806 768 : if (docosp) then
807 0 : call pio_seterrorhandling(File, PIO_BCAST_ERROR, err_handling)
808 0 : ierr = pio_inq_varid(File, 'cosp_cnt_init', vardesc)
809 0 : call pio_seterrorhandling(File, err_handling)
810 0 : if (ierr /= PIO_NOERR) then
811 0 : cosp_cnt_init = 0
812 : else
813 0 : ierr = pio_get_var(File, vardesc, cosp_cnt_init)
814 : end if
815 : end if
816 :
817 768 : ierr = pio_inq_varid(file, 'nextsw_cday', vardesc)
818 768 : ierr = pio_get_var(file, vardesc, nextsw_cday)
819 :
820 :
821 768 : end subroutine radiation_read_restart
822 :
823 : !===============================================================================
824 :
825 62675424 : subroutine radiation_tend( &
826 : state, ptend, pbuf, cam_out, cam_in, net_flx, rd_out)
827 :
828 : !-----------------------------------------------------------------------
829 : !
830 : ! CAM driver for radiation computation.
831 : !
832 : !-----------------------------------------------------------------------
833 :
834 : ! Location/Orbital Parameters for cosine zenith angle
835 : use phys_grid, only: get_rlat_all_p, get_rlon_all_p
836 : use cam_control_mod, only: eccen, mvelpp, lambm0, obliqr
837 : use shr_orb_mod, only: shr_orb_decl, shr_orb_cosz
838 :
839 : use rrtmgp_inputs, only: rrtmgp_set_state, rrtmgp_set_gases_lw, rrtmgp_set_cloud_lw, &
840 : rrtmgp_set_aer_lw, rrtmgp_set_gases_sw, rrtmgp_set_cloud_sw, &
841 : rrtmgp_set_aer_sw
842 :
843 : ! RRTMGP drivers for flux calculations.
844 : use mo_rte_lw, only: rte_lw
845 : use mo_rte_sw, only: rte_sw
846 :
847 : use radheat, only: radheat_tend
848 :
849 : use radiation_data, only: rad_data_write
850 :
851 : use interpolate_data, only: vertinterp
852 : use tropopause, only: tropopause_find_cam, TROP_ALG_HYBSTOB, TROP_ALG_CLIMATE
853 : use cospsimulator_intr, only: docosp, cospsimulator_intr_run, cosp_nradsteps
854 :
855 :
856 : ! Arguments
857 : type(physics_state), intent(in), target :: state
858 : type(physics_ptend), intent(out) :: ptend
859 : type(physics_buffer_desc), pointer :: pbuf(:)
860 : type(cam_out_t), intent(inout) :: cam_out
861 : type(cam_in_t), intent(in) :: cam_in
862 : real(r8), intent(out) :: net_flx(pcols)
863 :
864 : type(rad_out_t), target, optional, intent(out) :: rd_out
865 :
866 :
867 : ! Local variables
868 : type(rad_out_t), pointer :: rd ! allow rd_out to be optional by allocating a local object
869 : ! if the argument is not present
870 : logical :: write_output
871 :
872 : integer :: i, k, istat
873 : integer :: lchnk, ncol
874 : logical :: dosw, dolw
875 : integer :: icall ! loop index for climate/diagnostic radiation calls
876 :
877 : real(r8) :: calday ! current calendar day
878 : real(r8) :: delta ! Solar declination angle in radians
879 : real(r8) :: eccf ! Earth orbit eccentricity factor
880 : real(r8) :: clat(pcols) ! current latitudes(radians)
881 : real(r8) :: clon(pcols) ! current longitudes(radians)
882 : real(r8) :: coszrs(pcols) ! Cosine solar zenith angle
883 :
884 : ! Gathered indices of day and night columns
885 : ! chunk_column_index = IdxDay(daylight_column_index)
886 : integer :: Nday ! Number of daylight columns
887 : integer :: Nnite ! Number of night columns
888 : integer :: IdxDay(pcols) ! chunk indices of daylight columns
889 : integer :: IdxNite(pcols) ! chunk indices of night columns
890 :
891 : integer :: itim_old
892 :
893 1492272 : real(r8), pointer :: cld(:,:) ! cloud fraction
894 1492272 : real(r8), pointer :: cldfsnow(:,:) ! cloud fraction of just "snow clouds"
895 1492272 : real(r8), pointer :: cldfgrau(:,:) ! cloud fraction of just "graupel clouds"
896 : real(r8) :: cldfprime(pcols,pver) ! combined cloud fraction
897 1492272 : real(r8), pointer :: qrs(:,:) ! shortwave radiative heating rate
898 1492272 : real(r8), pointer :: qrl(:,:) ! longwave radiative heating rate
899 1492272 : real(r8), pointer :: fsds(:) ! Surface solar down flux
900 1492272 : real(r8), pointer :: fsns(:) ! Surface solar absorbed flux
901 1492272 : real(r8), pointer :: fsnt(:) ! Net column abs solar flux at model top
902 1492272 : real(r8), pointer :: flns(:) ! Srf longwave cooling (up-down) flux
903 1492272 : real(r8), pointer :: flnt(:) ! Net outgoing lw flux at model top
904 :
905 : real(r8), pointer, dimension(:,:,:) :: su => NULL() ! shortwave spectral flux up
906 : real(r8), pointer, dimension(:,:,:) :: sd => NULL() ! shortwave spectral flux down
907 : real(r8), pointer, dimension(:,:,:) :: lu => NULL() ! longwave spectral flux up
908 : real(r8), pointer, dimension(:,:,:) :: ld => NULL() ! longwave spectral flux down
909 :
910 : ! tropopause diagnostic
911 : integer :: troplev(pcols)
912 : real(r8) :: p_trop(pcols)
913 :
914 : ! state data passed to radiation calc
915 1492272 : real(r8), allocatable :: t_sfc(:)
916 1492272 : real(r8), allocatable :: emis_sfc(:,:)
917 1492272 : real(r8), allocatable :: t_rad(:,:)
918 1492272 : real(r8), allocatable :: pmid_rad(:,:)
919 1492272 : real(r8), allocatable :: pint_rad(:,:)
920 1492272 : real(r8), allocatable :: t_day(:,:)
921 1492272 : real(r8), allocatable :: pmid_day(:,:)
922 1492272 : real(r8), allocatable :: pint_day(:,:)
923 1492272 : real(r8), allocatable :: coszrs_day(:)
924 1492272 : real(r8), allocatable :: alb_dir(:,:)
925 1492272 : real(r8), allocatable :: alb_dif(:,:)
926 :
927 : ! in-cloud optical depths for COSP
928 : real(r8) :: cld_tau_cloudsim(pcols,pver) ! liq + ice
929 : real(r8) :: snow_tau_cloudsim(pcols,pver) ! snow
930 : real(r8) :: grau_tau_cloudsim(pcols,pver) ! graupel
931 : real(r8) :: cld_lw_abs_cloudsim(pcols,pver) ! liq + ice
932 : real(r8) :: snow_lw_abs_cloudsim(pcols,pver)! snow
933 : real(r8) :: grau_lw_abs_cloudsim(pcols,pver)! graupel
934 :
935 : ! Set vertical indexing in RRTMGP to be the same as CAM (top to bottom).
936 : logical, parameter :: top_at_1 = .true.
937 :
938 : ! TOA solar flux on RRTMGP g-points
939 1492272 : real(r8), allocatable :: toa_flux(:,:)
940 : ! TSI from RRTMGP data (from sum over g-point representation)
941 : real(r8) :: tsi_ref
942 :
943 : ! Planck sources for LW.
944 1492272 : type(ty_source_func_lw) :: sources_lw
945 :
946 : ! Gas volume mixing ratios. Use separate objects for LW and SW because SW only does
947 : ! calculations for daylight columns.
948 : ! These objects have a final method which deallocates the internal memory when they
949 : ! go out of scope (i.e., when radiation_tend returns), so no need for explicit deallocation.
950 1492272 : type(ty_gas_concs) :: gas_concs_lw
951 1492272 : type(ty_gas_concs) :: gas_concs_sw
952 :
953 : ! Atmosphere optics. This object is initialized with gas optics, then is incremented
954 : ! by the aerosol optics for the clear-sky radiative flux calculations, and then
955 : ! incremented again by the cloud optics for the all-sky radiative flux calculations.
956 1492272 : type(ty_optical_props_1scl) :: atm_optics_lw
957 1492272 : type(ty_optical_props_2str) :: atm_optics_sw
958 :
959 : ! Cloud optical properties objects (McICA sampling of cloud optical properties).
960 1492272 : type(ty_optical_props_1scl) :: cloud_lw
961 1492272 : type(ty_optical_props_2str) :: cloud_sw
962 :
963 : ! Aerosol optical properties objects.
964 1492272 : type(ty_optical_props_1scl) :: aer_lw
965 1492272 : type(ty_optical_props_2str) :: aer_sw
966 :
967 : ! Flux objects contain all fluxes computed by RRTMGP.
968 : ! SW allsky fluxes always include spectrally resolved fluxes needed for surface models.
969 : type(ty_fluxes_byband) :: fsw
970 : ! LW allsky fluxes only need spectrally resolved fluxes when spectralflux=.true.
971 : type(ty_fluxes_byband) :: flw
972 : ! Only broadband fluxes needed for clear sky (diagnostics).
973 : type(ty_fluxes_broadband) :: fswc, flwc
974 :
975 : ! Arrays for output diagnostics on CAM grid.
976 : real(r8) :: fns(pcols,pverp) ! net shortwave flux
977 : real(r8) :: fcns(pcols,pverp) ! net clear-sky shortwave flux
978 : real(r8) :: fnl(pcols,pverp) ! net longwave flux
979 : real(r8) :: fcnl(pcols,pverp) ! net clear-sky longwave flux
980 :
981 : ! for COSP
982 : real(r8) :: emis(pcols,pver) ! Cloud longwave emissivity
983 : real(r8) :: gb_snow_tau(pcols,pver) ! grid-box mean snow_tau
984 : real(r8) :: gb_snow_lw(pcols,pver) ! grid-box mean LW snow optical depth
985 :
986 : real(r8) :: ftem(pcols,pver) ! Temporary workspace for outfld variables
987 :
988 : character(len=128) :: errmsg
989 : character(len=*), parameter :: sub = 'radiation_tend'
990 : !--------------------------------------------------------------------------------------
991 :
992 1492272 : lchnk = state%lchnk
993 1492272 : ncol = state%ncol
994 :
995 1492272 : if (present(rd_out)) then
996 0 : rd => rd_out
997 0 : write_output = .false.
998 : else
999 1492272 : allocate(rd, stat=istat)
1000 1492272 : call handle_allocate_error(istat, sub, 'rd')
1001 1492272 : write_output = .true.
1002 : end if
1003 :
1004 1492272 : dosw = radiation_do('sw', get_nstep()) ! do shortwave radiation calc this timestep?
1005 1492272 : dolw = radiation_do('lw', get_nstep()) ! do longwave radiation calc this timestep?
1006 :
1007 : ! Cosine solar zenith angle for current time step
1008 1492272 : calday = get_curr_calday()
1009 1492272 : call get_rlat_all_p(lchnk, ncol, clat)
1010 1492272 : call get_rlon_all_p(lchnk, ncol, clon)
1011 :
1012 : call shr_orb_decl(calday, eccen, mvelpp, lambm0, obliqr, &
1013 1492272 : delta, eccf)
1014 :
1015 1492272 : if (use_rad_uniform_angle) then
1016 0 : do i = 1, ncol
1017 0 : coszrs(i) = shr_orb_cosz(calday, clat(i), clon(i), delta, dt_avg, &
1018 0 : uniform_angle=rad_uniform_angle)
1019 : end do
1020 : else
1021 24917472 : do i = 1, ncol
1022 : ! if dt_avg /= 0, it triggers using avg coszrs
1023 24917472 : coszrs(i) = shr_orb_cosz(calday, clat(i), clon(i), delta, dt_avg)
1024 : end do
1025 : end if
1026 :
1027 : ! Gather night/day column indices.
1028 1492272 : Nday = 0
1029 1492272 : Nnite = 0
1030 24917472 : do i = 1, ncol
1031 24917472 : if ( coszrs(i) > 0.0_r8 ) then
1032 11712600 : Nday = Nday + 1
1033 11712600 : IdxDay(Nday) = i
1034 : else
1035 11712600 : Nnite = Nnite + 1
1036 11712600 : IdxNite(Nnite) = i
1037 : end if
1038 : end do
1039 :
1040 : ! Associate pointers to physics buffer fields
1041 1492272 : itim_old = pbuf_old_tim_idx()
1042 1492272 : nullify(cldfsnow)
1043 1492272 : if (cldfsnow_idx > 0) then
1044 5969088 : call pbuf_get_field(pbuf, cldfsnow_idx, cldfsnow, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
1045 : end if
1046 1492272 : nullify(cldfgrau)
1047 1492272 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
1048 0 : call pbuf_get_field(pbuf, cldfgrau_idx, cldfgrau, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
1049 : endif
1050 :
1051 5969088 : call pbuf_get_field(pbuf, cld_idx, cld, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
1052 :
1053 1492272 : call pbuf_get_field(pbuf, qrs_idx, qrs)
1054 1492272 : call pbuf_get_field(pbuf, qrl_idx, qrl)
1055 :
1056 1492272 : call pbuf_get_field(pbuf, fsns_idx, fsns)
1057 1492272 : call pbuf_get_field(pbuf, fsnt_idx, fsnt)
1058 1492272 : call pbuf_get_field(pbuf, fsds_idx, fsds)
1059 1492272 : call pbuf_get_field(pbuf, flns_idx, flns)
1060 1492272 : call pbuf_get_field(pbuf, flnt_idx, flnt)
1061 :
1062 1492272 : if (spectralflux) then
1063 0 : call pbuf_get_field(pbuf, su_idx, su)
1064 0 : call pbuf_get_field(pbuf, sd_idx, sd)
1065 0 : call pbuf_get_field(pbuf, lu_idx, lu)
1066 0 : call pbuf_get_field(pbuf, ld_idx, ld)
1067 : end if
1068 :
1069 : ! Allocate the flux arrays and init to zero.
1070 1492272 : call initialize_rrtmgp_fluxes(nday, nlay+1, nswbands, fsw, do_direct=.true.)
1071 1492272 : call initialize_rrtmgp_fluxes(nday, nlay+1, nswbands, fswc, do_direct=.true.)
1072 1492272 : call initialize_rrtmgp_fluxes(ncol, nlay+1, nlwbands, flw)
1073 1492272 : call initialize_rrtmgp_fluxes(ncol, nlay+1, nlwbands, flwc)
1074 :
1075 : ! For CRM, make cloud equal to input observations:
1076 1492272 : if (scm_crm_mode .and. have_cld) then
1077 0 : do k = 1, pver
1078 0 : cld(:ncol,k)= cldobs(k)
1079 : end do
1080 : end if
1081 :
1082 : ! Find tropopause height if needed for diagnostic output
1083 1492272 : if (hist_fld_active('FSNR') .or. hist_fld_active('FLNR')) then
1084 : !REMOVECAM - no longer need this when CAM is retired and pcols no longer exists
1085 0 : troplev(:) = 0
1086 0 : p_trop(:) = 0._r8
1087 : !REMOVECAM_END
1088 : call tropopause_find_cam(state, troplev, tropP=p_trop, primary=TROP_ALG_HYBSTOB, &
1089 0 : backup=TROP_ALG_CLIMATE)
1090 : end if
1091 :
1092 : ! Get time of next radiation calculation - albedos will need to be
1093 : ! calculated by each surface model at this time
1094 1492272 : nextsw_cday = radiation_nextsw_cday()
1095 :
1096 1492272 : if (dosw .or. dolw) then
1097 :
1098 : allocate( &
1099 : t_sfc(ncol), emis_sfc(nlwbands,ncol), toa_flux(nday,nswgpts), &
1100 : t_rad(ncol,nlay), pmid_rad(ncol,nlay), pint_rad(ncol,nlay+1), &
1101 : t_day(nday,nlay), pmid_day(nday,nlay), pint_day(nday,nlay+1), &
1102 : coszrs_day(nday), alb_dir(nswbands,nday), alb_dif(nswbands,nday), &
1103 19139833 : stat=istat)
1104 746136 : call handle_allocate_error(istat, sub, 't_sfc,..,alb_dif')
1105 :
1106 : ! Prepares state variables, daylit columns, albedos for RRTMGP
1107 : call rrtmgp_set_state( &
1108 : state, cam_in, ncol, nlay, nday, &
1109 : idxday, coszrs, kdist_sw, t_sfc, emis_sfc, &
1110 : t_rad, pmid_rad, pint_rad, t_day, pmid_day, &
1111 746136 : pint_day, coszrs_day, alb_dir, alb_dif)
1112 :
1113 : ! Output the mass per layer, and total column burdens for gas and aerosol
1114 : ! constituents in the climate list.
1115 746136 : call rad_cnst_out(0, state, pbuf)
1116 :
1117 : ! Modified cloud fraction accounts for radiatively active snow and/or graupel
1118 746136 : call modified_cloud_fraction(ncol, cld, cldfsnow, cldfgrau, cldfprime)
1119 :
1120 : !========================!
1121 : ! SHORTWAVE calculations !
1122 : !========================!
1123 :
1124 746136 : if (dosw) then
1125 :
1126 : ! Set cloud optical properties in cloud_sw object.
1127 : call rrtmgp_set_cloud_sw( &
1128 : state, pbuf, nlay, nday, idxday, &
1129 : nnite, idxnite, pmid_day, cld, cldfsnow, &
1130 : cldfgrau, cldfprime, graupel_in_rad, kdist_sw, cloud_sw, &
1131 : rd%tot_cld_vistau, rd%tot_icld_vistau, rd%liq_icld_vistau, &
1132 : rd%ice_icld_vistau, rd%snow_icld_vistau, rd%grau_icld_vistau, &
1133 746136 : cld_tau_cloudsim, snow_tau_cloudsim, grau_tau_cloudsim )
1134 :
1135 746136 : if (write_output) then
1136 746136 : call radiation_output_cld(lchnk, rd)
1137 : end if
1138 :
1139 : ! If no daylight columns, can't create empty RRTMGP objects
1140 746136 : if (nday > 0) then
1141 :
1142 : ! Initialize object for gas concentrations.
1143 389263 : errmsg = gas_concs_sw%init(gaslist_lc)
1144 389263 : call stop_on_err(errmsg, sub, 'gas_concs_sw%init')
1145 :
1146 : ! Initialize object for combined gas + aerosol + cloud optics.
1147 : ! Allocates arrays for properties represented on g-points.
1148 389263 : errmsg = atm_optics_sw%alloc_2str(nday, nlay, kdist_sw)
1149 389263 : call stop_on_err(errmsg, sub, 'atm_optics_sw%alloc_2str')
1150 :
1151 : ! Initialize object for SW aerosol optics. Allocates arrays
1152 : ! for properties represented by band.
1153 389263 : errmsg = aer_sw%alloc_2str(nday, nlay, kdist_sw%get_band_lims_wavenumber())
1154 389263 : call stop_on_err(errmsg, sub, 'aer_sw%alloc_2str')
1155 :
1156 : end if
1157 :
1158 : ! The climate (icall==0) calculation must occur last.
1159 8953632 : do icall = N_DIAG, 0, -1
1160 8953632 : if (active_calls(icall)) then
1161 :
1162 746136 : if (nday > 0) then
1163 :
1164 : ! Set gas volume mixing ratios for this call in gas_concs_sw.
1165 : call rrtmgp_set_gases_sw( &
1166 : icall, state, pbuf, nlay, nday, &
1167 389263 : idxday, gas_concs_sw)
1168 :
1169 : ! Compute the gas optics (stored in atm_optics_sw).
1170 : ! toa_flux is the reference solar source from RRTMGP data.
1171 : errmsg = kdist_sw%gas_optics( &
1172 : pmid_day, pint_day, t_day, gas_concs_sw, atm_optics_sw, &
1173 389263 : toa_flux)
1174 389263 : call stop_on_err(errmsg, sub, 'kdist_sw%gas_optics')
1175 :
1176 : ! Scale the solar source
1177 43986719 : tsi_ref = sum(toa_flux(1,:))
1178 699892319 : toa_flux = toa_flux * sol_tsi * eccf / tsi_ref
1179 :
1180 : end if
1181 :
1182 : ! Set SW aerosol optical properties in the aer_sw object.
1183 : ! This call made even when no daylight columns because it does some
1184 : ! diagnostic aerosol output.
1185 : call rrtmgp_set_aer_sw( &
1186 746136 : icall, state, pbuf, nday, idxday, nnite, idxnite, aer_sw)
1187 :
1188 746136 : if (nday > 0) then
1189 :
1190 : ! Increment the gas optics (in atm_optics_sw) by the aerosol optics in aer_sw.
1191 389263 : errmsg = aer_sw%increment(atm_optics_sw)
1192 389263 : call stop_on_err(errmsg, sub, 'aer_sw%increment')
1193 :
1194 : ! Compute clear-sky fluxes.
1195 : errmsg = rte_sw(&
1196 : atm_optics_sw, top_at_1, coszrs_day, toa_flux, &
1197 389263 : alb_dir, alb_dif, fswc)
1198 389263 : call stop_on_err(errmsg, sub, 'clear-sky rte_sw')
1199 :
1200 : ! Increment the aerosol+gas optics (in atm_optics_sw) by the cloud optics in cloud_sw.
1201 389263 : errmsg = cloud_sw%increment(atm_optics_sw)
1202 389263 : call stop_on_err(errmsg, sub, 'cloud_sw%increment')
1203 :
1204 : ! Compute all-sky fluxes.
1205 : errmsg = rte_sw(&
1206 : atm_optics_sw, top_at_1, coszrs_day, toa_flux, &
1207 389263 : alb_dir, alb_dif, fsw)
1208 389263 : call stop_on_err(errmsg, sub, 'all-sky rte_sw')
1209 :
1210 : end if
1211 :
1212 : ! Transform RRTMGP outputs to CAM outputs and compute heating rates.
1213 746136 : call set_sw_diags()
1214 :
1215 746136 : if (write_output) then
1216 746136 : call radiation_output_sw(lchnk, ncol, icall, rd, pbuf, cam_out)
1217 : end if
1218 :
1219 : end if ! (active_calls(icall))
1220 : end do ! loop over diagnostic calcs (icall)
1221 :
1222 : else
1223 : ! SW calc not done. pbuf carries Q*dp across timesteps.
1224 : ! Convert to Q before calling radheat_tend.
1225 0 : qrs(:ncol,:) = qrs(:ncol,:) / state%pdel(:ncol,:)
1226 : end if ! if (dosw)
1227 :
1228 : !=======================!
1229 : ! LONGWAVE calculations !
1230 : !=======================!
1231 :
1232 746136 : if (dolw) then
1233 :
1234 : ! Initialize object for Planck sources.
1235 746136 : errmsg = sources_lw%alloc(ncol, nlay, kdist_lw)
1236 746136 : call stop_on_err(errmsg, sub, 'sources_lw%alloc')
1237 :
1238 : ! Set cloud optical properties in cloud_lw object.
1239 : call rrtmgp_set_cloud_lw( &
1240 : state, pbuf, ncol, nlay, nlaycam, &
1241 : cld, cldfsnow, cldfgrau, cldfprime, graupel_in_rad, &
1242 746136 : kdist_lw, cloud_lw, cld_lw_abs_cloudsim, snow_lw_abs_cloudsim, grau_lw_abs_cloudsim)
1243 :
1244 : ! Initialize object for gas concentrations
1245 746136 : errmsg = gas_concs_lw%init(gaslist_lc)
1246 746136 : call stop_on_err(errmsg, sub, 'gas_concs_lw%init')
1247 :
1248 : ! Initialize object for combined gas + aerosol + cloud optics.
1249 746136 : errmsg = atm_optics_lw%alloc_1scl(ncol, nlay, kdist_lw)
1250 746136 : call stop_on_err(errmsg, sub, 'atm_optics_lw%alloc_1scl')
1251 :
1252 : ! Initialize object for LW aerosol optics.
1253 746136 : errmsg = aer_lw%alloc_1scl(ncol, nlay, kdist_lw%get_band_lims_wavenumber())
1254 746136 : call stop_on_err(errmsg, sub, 'aer_lw%alloc_1scl')
1255 :
1256 : ! The climate (icall==0) calculation must occur last.
1257 8953632 : do icall = N_DIAG, 0, -1
1258 :
1259 8953632 : if (active_calls(icall)) then
1260 :
1261 : ! Set gas volume mixing ratios for this call in gas_concs_lw.
1262 746136 : call rrtmgp_set_gases_lw(icall, state, pbuf, nlay, gas_concs_lw)
1263 :
1264 : ! Compute the gas optics and Planck sources.
1265 : errmsg = kdist_lw%gas_optics( &
1266 : pmid_rad, pint_rad, t_rad, t_sfc, gas_concs_lw, &
1267 746136 : atm_optics_lw, sources_lw)
1268 746136 : call stop_on_err(errmsg, sub, 'kdist_lw%gas_optics')
1269 :
1270 : ! Set LW aerosol optical properties in the aer_lw object.
1271 746136 : call rrtmgp_set_aer_lw(icall, state, pbuf, aer_lw)
1272 :
1273 : ! Increment the gas optics by the aerosol optics.
1274 746136 : errmsg = aer_lw%increment(atm_optics_lw)
1275 746136 : call stop_on_err(errmsg, sub, 'aer_lw%increment')
1276 :
1277 : ! Compute clear-sky LW fluxes
1278 746136 : errmsg = rte_lw(atm_optics_lw, top_at_1, sources_lw, emis_sfc, flwc)
1279 746136 : call stop_on_err(errmsg, sub, 'clear-sky rte_lw')
1280 :
1281 : ! Increment the gas+aerosol optics by the cloud optics.
1282 746136 : errmsg = cloud_lw%increment(atm_optics_lw)
1283 746136 : call stop_on_err(errmsg, sub, 'cloud_lw%increment')
1284 :
1285 : ! Compute all-sky LW fluxes
1286 746136 : errmsg = rte_lw(atm_optics_lw, top_at_1, sources_lw, emis_sfc, flw)
1287 746136 : call stop_on_err(errmsg, sub, 'all-sky rte_lw')
1288 :
1289 : ! Transform RRTMGP outputs to CAM outputs and compute heating rates.
1290 746136 : call set_lw_diags()
1291 :
1292 746136 : if (write_output) then
1293 746136 : call radiation_output_lw(lchnk, ncol, icall, rd, pbuf, cam_out)
1294 : end if
1295 :
1296 : end if ! (active_calls(icall))
1297 : end do ! loop over diagnostic calcs (icall)
1298 :
1299 : else
1300 : ! LW calc not done. pbuf carries Q*dp across timesteps.
1301 : ! Convert to Q before calling radheat_tend.
1302 0 : qrl(:ncol,:) = qrl(:ncol,:) / state%pdel(:ncol,:)
1303 : end if ! if (dolw)
1304 :
1305 0 : deallocate( &
1306 0 : t_sfc, emis_sfc, toa_flux, t_rad, pmid_rad, pint_rad, &
1307 746136 : t_day, pmid_day, pint_day, coszrs_day, alb_dir, alb_dif)
1308 :
1309 : !================!
1310 : ! COSP simulator !
1311 : !================!
1312 :
1313 746136 : if (docosp) then
1314 :
1315 0 : emis(:,:) = 0._r8
1316 0 : emis(:ncol,:) = 1._r8 - exp(-cld_lw_abs_cloudsim(:ncol,:))
1317 0 : call outfld('EMIS', emis, pcols, lchnk)
1318 :
1319 : ! compute grid-box mean SW and LW snow optical depth for use by COSP
1320 0 : gb_snow_tau(:,:) = 0._r8
1321 0 : gb_snow_lw(:,:) = 0._r8
1322 0 : if (cldfsnow_idx > 0) then
1323 0 : do i = 1, ncol
1324 0 : do k = 1, pver
1325 0 : if (cldfsnow(i,k) > 0._r8) then
1326 :
1327 : ! Add graupel to snow tau for cosp
1328 0 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
1329 0 : gb_snow_tau(i,k) = snow_tau_cloudsim(i,k)*cldfsnow(i,k) + &
1330 0 : grau_tau_cloudsim(i,k)*cldfgrau(i,k)
1331 : gb_snow_lw(i,k) = snow_lw_abs_cloudsim(i,k)*cldfsnow(i,k) + &
1332 0 : grau_lw_abs_cloudsim(i,k)*cldfgrau(i,k)
1333 : else
1334 0 : gb_snow_tau(i,k) = snow_tau_cloudsim(i,k)*cldfsnow(i,k)
1335 0 : gb_snow_lw(i,k) = snow_lw_abs_cloudsim(i,k)*cldfsnow(i,k)
1336 : end if
1337 : end if
1338 : end do
1339 : end do
1340 : end if
1341 :
1342 : ! advance counter for this timestep (chunk dimension required for thread safety)
1343 0 : cosp_cnt(lchnk) = cosp_cnt(lchnk) + 1
1344 :
1345 : ! if counter is the same as cosp_nradsteps, run cosp and reset counter
1346 0 : if (cosp_nradsteps .eq. cosp_cnt(lchnk)) then
1347 :
1348 : ! N.B.: For snow optical properties, the GRID-BOX MEAN shortwave and longwave
1349 : ! optical depths are passed.
1350 : call cospsimulator_intr_run( &
1351 : state, pbuf, cam_in, emis, coszrs, &
1352 : cld_swtau_in=cld_tau_cloudsim, snow_tau_in=gb_snow_tau, &
1353 0 : snow_emis_in=gb_snow_lw)
1354 0 : cosp_cnt(lchnk) = 0
1355 : end if
1356 : end if ! docosp
1357 :
1358 : else
1359 : ! Radiative flux calculations not done. The quantity Q*dp is carried by the
1360 : ! physics buffer across timesteps. It must be converted to Q (dry static energy
1361 : ! tendency) before being passed to radheat_tend.
1362 1159408584 : qrs(:ncol,:) = qrs(:ncol,:) / state%pdel(:ncol,:)
1363 1159408584 : qrl(:ncol,:) = qrl(:ncol,:) / state%pdel(:ncol,:)
1364 :
1365 : end if ! if (dosw .or. dolw) then
1366 :
1367 : ! Output for PORT: Parallel Offline Radiative Transport
1368 1492272 : call rad_data_write(pbuf, state, cam_in, coszrs)
1369 :
1370 : ! Compute net radiative heating tendency. Note that the WACCM version
1371 : ! of radheat_tend merges upper atmosphere heating rates with those calculated
1372 : ! by RRTMGP.
1373 : call radheat_tend(state, pbuf, ptend, qrl, qrs, fsns, &
1374 1492272 : fsnt, flns, flnt, cam_in%asdir, net_flx)
1375 :
1376 1492272 : if (write_output) then
1377 : ! Compute heating rate for dtheta/dt
1378 140273568 : do k = 1, pver
1379 2318817168 : do i = 1, ncol
1380 2317324896 : ftem(i,k) = (qrs(i,k) + qrl(i,k))/cpair * (1.e5_r8/state%pmid(i,k))**cappa
1381 : end do
1382 : end do
1383 1492272 : call outfld('HR', ftem, pcols, lchnk)
1384 : end if
1385 :
1386 : ! The radiative heating rates are carried in the physics buffer across timesteps
1387 : ! as Q*dp (for energy conservation).
1388 2318817168 : qrs(:ncol,:) = qrs(:ncol,:) * state%pdel(:ncol,:)
1389 2318817168 : qrl(:ncol,:) = qrl(:ncol,:) * state%pdel(:ncol,:)
1390 :
1391 1492272 : if (.not. present(rd_out)) then
1392 1492272 : deallocate(rd)
1393 : end if
1394 1492272 : call free_optics_sw(atm_optics_sw)
1395 1492272 : call free_optics_sw(cloud_sw)
1396 1492272 : call free_optics_sw(aer_sw)
1397 1492272 : call free_fluxes(fsw)
1398 1492272 : call free_fluxes(fswc)
1399 :
1400 1492272 : call sources_lw%finalize()
1401 1492272 : call free_optics_lw(cloud_lw)
1402 1492272 : call free_optics_lw(aer_lw)
1403 1492272 : call free_fluxes(flw)
1404 2984544 : call free_fluxes(flwc)
1405 :
1406 : !-------------------------------------------------------------------------------
1407 : contains
1408 : !-------------------------------------------------------------------------------
1409 :
1410 746136 : subroutine set_sw_diags()
1411 :
1412 : ! Transform RRTMGP output for CAM and compute heating rates.
1413 : ! SW fluxes from RRTMGP are on daylight columns only, so expand to
1414 : ! full chunks for output to CAM history.
1415 :
1416 : integer :: i
1417 : real(r8), dimension(size(fsw%bnd_flux_dn,1), &
1418 : size(fsw%bnd_flux_dn,2), &
1419 1492272 : size(fsw%bnd_flux_dn,3)) :: flux_dn_diffuse
1420 : !-------------------------------------------------------------------------
1421 :
1422 : ! Initialize to provide 0.0 values for night columns.
1423 746136 : fns = 0._r8 ! net sw flux
1424 746136 : fcns = 0._r8 ! net sw clearsky flux
1425 12684312 : fsds = 0._r8 ! downward sw flux at surface
1426 12684312 : rd%fsdsc = 0._r8 ! downward sw clearsky flux at surface
1427 12684312 : rd%fsutoa = 0._r8 ! upward sw flux at TOA
1428 12684312 : rd%fsntoa = 0._r8 ! net sw at TOA
1429 12684312 : rd%fsntoac = 0._r8 ! net sw clearsky flux at TOA
1430 12684312 : rd%solin = 0._r8 ! solar irradiance at TOA
1431 1193071464 : rd%flux_sw_up = 0._r8
1432 1193071464 : rd%flux_sw_dn = 0._r8
1433 1193071464 : rd%flux_sw_clr_up = 0._r8
1434 1193071464 : rd%flux_sw_clr_dn = 0._r8
1435 :
1436 1180387152 : qrs = 0._r8
1437 12684312 : fsns = 0._r8
1438 12684312 : fsnt = 0._r8
1439 1180387152 : rd%qrsc = 0._r8
1440 12684312 : rd%fsnsc = 0._r8
1441 12684312 : rd%fsntc = 0._r8
1442 :
1443 6602436 : do i = 1, nday
1444 556348500 : fns(idxday(i),ktopcam:) = fsw%flux_net(i, ktoprad:)
1445 556348500 : fcns(idxday(i),ktopcam:) = fswc%flux_net(i,ktoprad:)
1446 5856300 : fsds(idxday(i)) = fsw%flux_dn(i, nlay+1)
1447 5856300 : rd%fsdsc(idxday(i)) = fswc%flux_dn(i, nlay+1)
1448 5856300 : rd%fsutoa(idxday(i)) = fsw%flux_up(i, 1)
1449 5856300 : rd%fsntoa(idxday(i)) = fsw%flux_net(i, 1)
1450 5856300 : rd%fsntoac(idxday(i)) = fswc%flux_net(i, 1)
1451 5856300 : rd%solin(idxday(i)) = fswc%flux_dn(i, 1)
1452 556348500 : rd%flux_sw_up(idxday(i),ktopcam:) = fsw%flux_up(i,ktoprad:)
1453 556348500 : rd%flux_sw_dn(idxday(i),ktopcam:) = fsw%flux_dn(i,ktoprad:)
1454 556348500 : rd%flux_sw_clr_up(idxday(i),ktopcam:) = fswc%flux_up(i,ktoprad:)
1455 557094636 : rd%flux_sw_clr_dn(idxday(i),ktopcam:) = fswc%flux_dn(i,ktoprad:)
1456 : end do
1457 :
1458 : ! Compute heating rate as a dry static energy tendency.
1459 746136 : call heating_rate('SW', ncol, fns, qrs)
1460 746136 : call heating_rate('SW', ncol, fcns, rd%qrsc)
1461 :
1462 12458736 : fsns(:ncol) = fns(:ncol,pverp) ! net sw flux at surface
1463 12458736 : fsnt(:ncol) = fns(:ncol,ktopcam) ! net sw flux at top-of-model (w/o extra layer)
1464 13204872 : rd%fsnsc(:ncol) = fcns(:ncol,pverp) ! net sw clearsky flux at surface
1465 13204872 : rd%fsntc(:ncol) = fcns(:ncol,ktopcam) ! net sw clearsky flux at top
1466 :
1467 12458736 : cam_out%netsw(:ncol) = fsns(:ncol)
1468 :
1469 : ! Output fluxes at 200 mb
1470 746136 : call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fns, rd%fsn200)
1471 746136 : call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fcns, rd%fsn200c)
1472 746136 : if (hist_fld_active('FSNR')) then
1473 0 : do i = 1,ncol
1474 0 : call vertinterp(1, 1, pverp, state%pint(i,:), p_trop(i), fns(i,:), rd%fsnr(i))
1475 : end do
1476 : end if
1477 :
1478 746136 : if (spectralflux) then
1479 0 : su = 0._r8
1480 0 : sd = 0._r8
1481 0 : do i = 1, nday
1482 0 : su(idxday(i),ktopcam:,:) = fsw%bnd_flux_up(i,ktoprad:,:)
1483 0 : sd(idxday(i),ktopcam:,:) = fsw%bnd_flux_dn(i,ktoprad:,:)
1484 : end do
1485 : end if
1486 :
1487 : ! Export surface fluxes
1488 : ! sols(pcols) Direct solar rad on surface (< 0.7)
1489 : ! soll(pcols) Direct solar rad on surface (>= 0.7)
1490 : ! RRTMGP: Near-IR bands (1-10), 820-16000 cm-1, 0.625-12.195 microns
1491 : ! Put half of band 10 in each of the UV/visible and near-IR values,
1492 : ! since this band straddles 0.7 microns:
1493 : ! UV/visible bands 10-13, 16000-50000 cm-1, 0.200-0.625 micron
1494 :
1495 : ! reset fluxes
1496 12684312 : cam_out%sols = 0.0_r8
1497 12684312 : cam_out%soll = 0.0_r8
1498 12684312 : cam_out%solsd = 0.0_r8
1499 12684312 : cam_out%solld = 0.0_r8
1500 :
1501 : ! Calculate diffuse flux from total and direct
1502 8792431920 : flux_dn_diffuse = fsw%bnd_flux_dn - fsw%bnd_flux_dn_dir
1503 :
1504 6602436 : do i = 1, nday
1505 17568900 : cam_out%soll(idxday(i)) = sum(fsw%bnd_flux_dn_dir(i,nlay+1,1:9)) &
1506 76131900 : + 0.5_r8 * fsw%bnd_flux_dn_dir(i,nlay+1,10)
1507 :
1508 5856300 : cam_out%sols(idxday(i)) = 0.5_r8 * fsw%bnd_flux_dn_dir(i,nlay+1,10) &
1509 35137800 : + sum(fsw%bnd_flux_dn_dir(i,nlay+1,11:14))
1510 :
1511 11712600 : cam_out%solld(idxday(i)) = sum(flux_dn_diffuse(i,nlay+1,1:9)) &
1512 70275600 : + 0.5_r8 * flux_dn_diffuse(i,nlay+1,10)
1513 :
1514 5856300 : cam_out%solsd(idxday(i)) = 0.5_r8 * flux_dn_diffuse(i, nlay+1, 10) &
1515 41740236 : + sum(flux_dn_diffuse(i,nlay+1,11:14))
1516 : end do
1517 :
1518 1492272 : end subroutine set_sw_diags
1519 :
1520 : !-------------------------------------------------------------------------------
1521 :
1522 746136 : subroutine set_lw_diags()
1523 :
1524 : ! Set CAM LW diagnostics
1525 : !----------------------------------------------------------------------------
1526 :
1527 746136 : fnl = 0._r8
1528 746136 : fcnl = 0._r8
1529 :
1530 : ! RTE-RRTMGP convention for net is (down - up) **CAM assumes (up - down) !!
1531 1171867320 : fnl(:ncol,ktopcam:) = -1._r8 * flw%flux_net( :, ktoprad:)
1532 1171867320 : fcnl(:ncol,ktopcam:) = -1._r8 * flwc%flux_net( :, ktoprad:)
1533 :
1534 1171867320 : rd%flux_lw_up(:ncol,ktopcam:) = flw%flux_up( :, ktoprad:)
1535 1171867320 : rd%flux_lw_clr_up(:ncol,ktopcam:) = flwc%flux_up(:, ktoprad:)
1536 1171867320 : rd%flux_lw_dn(:ncol,ktopcam:) = flw%flux_dn( :, ktoprad:)
1537 1171867320 : rd%flux_lw_clr_dn(:ncol,ktopcam:) = flwc%flux_dn(:, ktoprad:)
1538 :
1539 746136 : call heating_rate('LW', ncol, fnl, qrl)
1540 746136 : call heating_rate('LW', ncol, fcnl, rd%qrlc)
1541 :
1542 12458736 : flns(:ncol) = fnl(:ncol, pverp)
1543 12458736 : flnt(:ncol) = fnl(:ncol, ktopcam)
1544 :
1545 13204872 : rd%flnsc(:ncol) = fcnl(:ncol, pverp)
1546 13204872 : rd%flntc(:ncol) = fcnl(:ncol, ktopcam) ! net lw flux at top-of-model
1547 :
1548 12458736 : cam_out%flwds(:ncol) = flw%flux_dn(:, nlay+1)
1549 12458736 : rd%fldsc(:ncol) = flwc%flux_dn(:, nlay+1)
1550 :
1551 12458736 : rd%flut(:ncol) = flw%flux_up(:, ktoprad)
1552 12458736 : rd%flutc(:ncol) = flwc%flux_up(:, ktoprad)
1553 :
1554 : ! Output fluxes at 200 mb
1555 746136 : call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fnl, rd%fln200)
1556 746136 : call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fcnl, rd%fln200c)
1557 746136 : if (hist_fld_active('FLNR')) then
1558 0 : do i = 1,ncol
1559 0 : call vertinterp(1, 1, pverp, state%pint(i,:), p_trop(i), fnl(i,:), rd%flnr(i))
1560 : end do
1561 : end if
1562 :
1563 746136 : if (spectralflux) then
1564 0 : lu = 0._r8
1565 0 : ld = 0._r8
1566 0 : lu(:ncol, ktopcam:, :) = flw%bnd_flux_up(:, ktoprad:, :)
1567 0 : ld(:ncol, ktopcam:, :) = flw%bnd_flux_dn(:, ktoprad:, :)
1568 : end if
1569 :
1570 746136 : end subroutine set_lw_diags
1571 :
1572 : !-------------------------------------------------------------------------------
1573 :
1574 2984544 : subroutine heating_rate(type, ncol, flux_net, hrate)
1575 :
1576 : ! Compute heating rate as a dry static energy tendency
1577 :
1578 : ! arguments
1579 : character(2), intent(in) :: type ! either LW or SW
1580 : integer, intent(in) :: ncol
1581 : real(r8), intent(in) :: flux_net(pcols,pverp) ! W/m^2
1582 : real(r8), intent(out) :: hrate(pcols,pver) ! J/kg/s
1583 :
1584 : ! local vars
1585 : integer :: k
1586 :
1587 : ! Initialize for layers where RRTMGP is not providing fluxes.
1588 2984544 : hrate = 0.0_r8
1589 :
1590 1492272 : select case (type)
1591 : case ('LW')
1592 :
1593 140273568 : do k = ktopcam, pver
1594 : ! (flux divergence as bottom-MINUS-top) * g/dp
1595 277562592 : hrate(:ncol,k) = (flux_net(:ncol,k+1) - flux_net(:ncol,k)) * &
1596 2596379760 : gravit * state%rpdel(:ncol,k)
1597 : end do
1598 :
1599 : case ('SW')
1600 :
1601 141765840 : do k = ktopcam, pver
1602 : ! top - bottom
1603 277562592 : hrate(:ncol,k) = (flux_net(:ncol,k) - flux_net(:ncol,k+1)) * &
1604 2596379760 : gravit * state%rpdel(:ncol,k)
1605 : end do
1606 :
1607 : end select
1608 :
1609 2984544 : end subroutine heating_rate
1610 :
1611 : !----------------------------------------------------------------------------
1612 : ! -- end contains statement of radiation_tend --
1613 : !----------------------------------------------------------------------------
1614 : end subroutine radiation_tend
1615 :
1616 : !===============================================================================
1617 :
1618 746136 : subroutine radiation_output_sw(lchnk, ncol, icall, rd, pbuf, cam_out)
1619 :
1620 : ! Dump shortwave radiation information to history buffer.
1621 :
1622 : integer , intent(in) :: lchnk
1623 : integer, intent(in) :: ncol
1624 : integer, intent(in) :: icall
1625 : type(rad_out_t), intent(in) :: rd
1626 : type(physics_buffer_desc), pointer :: pbuf(:)
1627 : type(cam_out_t), intent(in) :: cam_out
1628 :
1629 : ! local variables
1630 746136 : real(r8), pointer :: qrs(:,:)
1631 746136 : real(r8), pointer :: fsnt(:)
1632 746136 : real(r8), pointer :: fsns(:)
1633 746136 : real(r8), pointer :: fsds(:)
1634 :
1635 : real(r8) :: ftem(pcols)
1636 : !----------------------------------------------------------------------------
1637 :
1638 746136 : call pbuf_get_field(pbuf, qrs_idx, qrs)
1639 746136 : call pbuf_get_field(pbuf, fsnt_idx, fsnt)
1640 746136 : call pbuf_get_field(pbuf, fsns_idx, fsns)
1641 746136 : call pbuf_get_field(pbuf, fsds_idx, fsds)
1642 :
1643 746136 : call outfld('SOLIN'//diag(icall), rd%solin, pcols, lchnk)
1644 :
1645 : ! QRS is output as temperature tendency.
1646 1159408584 : call outfld('QRS'//diag(icall), qrs(:ncol,:)/cpair, ncol, lchnk)
1647 1159408584 : call outfld('QRSC'//diag(icall), rd%qrsc(:ncol,:)/cpair, ncol, lchnk)
1648 :
1649 746136 : call outfld('FSNT'//diag(icall), fsnt, pcols, lchnk)
1650 746136 : call outfld('FSNTC'//diag(icall), rd%fsntc, pcols, lchnk)
1651 746136 : call outfld('FSNTOA'//diag(icall), rd%fsntoa, pcols, lchnk)
1652 746136 : call outfld('FSNTOAC'//diag(icall), rd%fsntoac, pcols, lchnk)
1653 :
1654 12458736 : ftem(:ncol) = rd%fsntoa(:ncol) - rd%fsntoac(:ncol)
1655 746136 : call outfld('SWCF'//diag(icall), ftem, pcols, lchnk)
1656 :
1657 746136 : call outfld('FSUTOA'//diag(icall), rd%fsutoa, pcols, lchnk)
1658 :
1659 746136 : call outfld('FSNIRTOA'//diag(icall), rd%fsnirt, pcols, lchnk)
1660 746136 : call outfld('FSNRTOAC'//diag(icall), rd%fsnrtc, pcols, lchnk)
1661 746136 : call outfld('FSNRTOAS'//diag(icall), rd%fsnirtsq, pcols, lchnk)
1662 :
1663 746136 : call outfld('FSN200'//diag(icall), rd%fsn200, pcols, lchnk)
1664 746136 : call outfld('FSN200C'//diag(icall), rd%fsn200c, pcols, lchnk)
1665 :
1666 746136 : call outfld('FSNR'//diag(icall), rd%fsnr, pcols, lchnk)
1667 :
1668 746136 : call outfld('SOLS'//diag(icall), cam_out%sols, pcols, lchnk)
1669 746136 : call outfld('SOLL'//diag(icall), cam_out%soll, pcols, lchnk)
1670 746136 : call outfld('SOLSD'//diag(icall), cam_out%solsd, pcols, lchnk)
1671 746136 : call outfld('SOLLD'//diag(icall), cam_out%solld, pcols, lchnk)
1672 :
1673 746136 : call outfld('FSNS'//diag(icall), fsns, pcols, lchnk)
1674 746136 : call outfld('FSNSC'//diag(icall), rd%fsnsc, pcols, lchnk)
1675 :
1676 746136 : call outfld('FSDS'//diag(icall), fsds, pcols, lchnk)
1677 746136 : call outfld('FSDSC'//diag(icall), rd%fsdsc, pcols, lchnk)
1678 :
1679 746136 : call outfld('FUS'//diag(icall), rd%flux_sw_up, pcols, lchnk)
1680 746136 : call outfld('FUSC'//diag(icall), rd%flux_sw_clr_up, pcols, lchnk)
1681 746136 : call outfld('FDS'//diag(icall), rd%flux_sw_dn, pcols, lchnk)
1682 746136 : call outfld('FDSC'//diag(icall), rd%flux_sw_clr_dn, pcols, lchnk)
1683 :
1684 746136 : end subroutine radiation_output_sw
1685 :
1686 : !===============================================================================
1687 :
1688 746136 : subroutine radiation_output_cld(lchnk, rd)
1689 :
1690 : ! Dump shortwave cloud optics information to history buffer.
1691 :
1692 : integer , intent(in) :: lchnk
1693 : type(rad_out_t), intent(in) :: rd
1694 : !----------------------------------------------------------------------------
1695 :
1696 746136 : call outfld('TOT_CLD_VISTAU', rd%tot_cld_vistau, pcols, lchnk)
1697 746136 : call outfld('TOT_ICLD_VISTAU', rd%tot_icld_vistau, pcols, lchnk)
1698 746136 : call outfld('LIQ_ICLD_VISTAU', rd%liq_icld_vistau, pcols, lchnk)
1699 746136 : call outfld('ICE_ICLD_VISTAU', rd%ice_icld_vistau, pcols, lchnk)
1700 746136 : if (cldfsnow_idx > 0) then
1701 746136 : call outfld('SNOW_ICLD_VISTAU', rd%snow_icld_vistau, pcols, lchnk)
1702 : endif
1703 746136 : if (cldfgrau_idx > 0 .and. graupel_in_rad) then
1704 0 : call outfld('GRAU_ICLD_VISTAU', rd%grau_icld_vistau , pcols, lchnk)
1705 : endif
1706 :
1707 746136 : end subroutine radiation_output_cld
1708 :
1709 : !===============================================================================
1710 :
1711 746136 : subroutine radiation_output_lw(lchnk, ncol, icall, rd, pbuf, cam_out)
1712 :
1713 : ! Dump longwave radiation information to history buffer
1714 :
1715 : integer, intent(in) :: lchnk
1716 : integer, intent(in) :: ncol
1717 : integer, intent(in) :: icall ! icall=0 for climate diagnostics
1718 : type(rad_out_t), intent(in) :: rd
1719 : type(physics_buffer_desc), pointer :: pbuf(:)
1720 : type(cam_out_t), intent(in) :: cam_out
1721 :
1722 : ! local variables
1723 746136 : real(r8), pointer :: qrl(:,:)
1724 746136 : real(r8), pointer :: flnt(:)
1725 746136 : real(r8), pointer :: flns(:)
1726 :
1727 : real(r8) :: ftem(pcols)
1728 : !----------------------------------------------------------------------------
1729 :
1730 746136 : call pbuf_get_field(pbuf, qrl_idx, qrl)
1731 746136 : call pbuf_get_field(pbuf, flnt_idx, flnt)
1732 746136 : call pbuf_get_field(pbuf, flns_idx, flns)
1733 :
1734 1159408584 : call outfld('QRL'//diag(icall), qrl(:ncol,:)/cpair, ncol, lchnk)
1735 1159408584 : call outfld('QRLC'//diag(icall), rd%qrlc(:ncol,:)/cpair, ncol, lchnk)
1736 :
1737 746136 : call outfld('FLNT'//diag(icall), flnt, pcols, lchnk)
1738 746136 : call outfld('FLNTC'//diag(icall), rd%flntc, pcols, lchnk)
1739 :
1740 746136 : call outfld('FLUT'//diag(icall), rd%flut, pcols, lchnk)
1741 746136 : call outfld('FLUTC'//diag(icall), rd%flutc, pcols, lchnk)
1742 :
1743 12458736 : ftem(:ncol) = rd%flutc(:ncol) - rd%flut(:ncol)
1744 746136 : call outfld('LWCF'//diag(icall), ftem, pcols, lchnk)
1745 :
1746 746136 : call outfld('FLN200'//diag(icall), rd%fln200, pcols, lchnk)
1747 746136 : call outfld('FLN200C'//diag(icall), rd%fln200c, pcols, lchnk)
1748 :
1749 746136 : call outfld('FLNR'//diag(icall), rd%flnr, pcols, lchnk)
1750 :
1751 746136 : call outfld('FLNS'//diag(icall), flns, pcols, lchnk)
1752 746136 : call outfld('FLNSC'//diag(icall), rd%flnsc, pcols, lchnk)
1753 :
1754 746136 : call outfld('FLDS'//diag(icall), cam_out%flwds, pcols, lchnk)
1755 746136 : call outfld('FLDSC'//diag(icall), rd%fldsc, pcols, lchnk)
1756 :
1757 746136 : call outfld('FDL'//diag(icall), rd%flux_lw_dn, pcols, lchnk)
1758 746136 : call outfld('FDLC'//diag(icall), rd%flux_lw_clr_dn, pcols, lchnk)
1759 746136 : call outfld('FUL'//diag(icall), rd%flux_lw_up, pcols, lchnk)
1760 746136 : call outfld('FULC'//diag(icall), rd%flux_lw_clr_up, pcols, lchnk)
1761 :
1762 746136 : end subroutine radiation_output_lw
1763 :
1764 : !===============================================================================
1765 :
1766 3072 : subroutine coefs_init(coefs_file, available_gases, kdist)
1767 :
1768 : ! Read data from coefficients file. Initialize the kdist object.
1769 : ! available_gases object provides the gas names that CAM provides.
1770 :
1771 : ! arguments
1772 : character(len=*), intent(in) :: coefs_file
1773 : class(ty_gas_concs), intent(in) :: available_gases
1774 : class(ty_gas_optics_rrtmgp), intent(out) :: kdist
1775 :
1776 : ! local variables
1777 : type(file_desc_t) :: fh ! pio file handle
1778 : character(len=cl) :: locfn ! path to file on local storage
1779 :
1780 : ! File dimensions
1781 : integer :: &
1782 : absorber, &
1783 : atmos_layer, &
1784 : bnd, &
1785 : pressure, &
1786 : temperature, &
1787 : absorber_ext, &
1788 : pressure_interp, &
1789 : mixing_fraction, &
1790 : gpt, &
1791 : temperature_Planck
1792 :
1793 : integer :: i
1794 : integer :: did, vid
1795 : integer :: ierr, istat
1796 :
1797 3072 : character(32), dimension(:), allocatable :: gas_names
1798 3072 : integer, dimension(:,:,:), allocatable :: key_species
1799 3072 : integer, dimension(:,:), allocatable :: band2gpt
1800 3072 : real(r8), dimension(:,:), allocatable :: band_lims_wavenum
1801 3072 : real(r8), dimension(:), allocatable :: press_ref, temp_ref
1802 : real(r8) :: press_ref_trop, temp_ref_t, temp_ref_p
1803 3072 : real(r8), dimension(:,:,:), allocatable :: vmr_ref
1804 3072 : real(r8), dimension(:,:,:,:), allocatable :: kmajor
1805 3072 : real(r8), dimension(:,:,:), allocatable :: kminor_lower, kminor_upper
1806 3072 : real(r8), dimension(:,:), allocatable :: totplnk
1807 3072 : real(r8), dimension(:,:,:,:), allocatable :: planck_frac
1808 3072 : real(r8), dimension(:), allocatable :: solar_src_quiet, solar_src_facular, solar_src_sunspot
1809 : real(r8) :: tsi_default
1810 3072 : real(r8), dimension(:,:,:), allocatable :: rayl_lower, rayl_upper
1811 3072 : character(len=32), dimension(:), allocatable :: gas_minor, &
1812 3072 : identifier_minor, &
1813 3072 : minor_gases_lower, &
1814 3072 : minor_gases_upper, &
1815 3072 : scaling_gas_lower, &
1816 3072 : scaling_gas_upper
1817 3072 : integer, dimension(:,:), allocatable :: minor_limits_gpt_lower, &
1818 3072 : minor_limits_gpt_upper
1819 : ! Send these to RRTMGP as logicals,
1820 : ! but they have to be read from the netCDF as integers
1821 3072 : logical, dimension(:), allocatable :: minor_scales_with_density_lower, &
1822 3072 : minor_scales_with_density_upper
1823 3072 : logical, dimension(:), allocatable :: scale_by_complement_lower, &
1824 3072 : scale_by_complement_upper
1825 3072 : integer, dimension(:), allocatable :: int2log ! use this to convert integer-to-logical.
1826 3072 : integer, dimension(:), allocatable :: kminor_start_lower, kminor_start_upper
1827 3072 : real(r8), dimension(:,:), allocatable :: optimal_angle_fit
1828 : real(r8) :: mg_default, sb_default
1829 :
1830 : integer :: pairs, &
1831 : minorabsorbers, &
1832 : minor_absorber_intervals_lower, &
1833 : minor_absorber_intervals_upper, &
1834 : contributors_lower, &
1835 : contributors_upper, &
1836 : fit_coeffs
1837 :
1838 : character(len=128) :: error_msg
1839 : character(len=*), parameter :: sub = 'coefs_init'
1840 : !----------------------------------------------------------------------------
1841 :
1842 : ! Open file
1843 3072 : call getfil(coefs_file, locfn, 0)
1844 3072 : call cam_pio_openfile(fh, locfn, PIO_NOWRITE)
1845 :
1846 3072 : call pio_seterrorhandling(fh, PIO_BCAST_ERROR)
1847 :
1848 : ! Get dimensions
1849 :
1850 3072 : ierr = pio_inq_dimid(fh, 'absorber', did)
1851 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': absorber not found')
1852 3072 : ierr = pio_inq_dimlen(fh, did, absorber)
1853 :
1854 3072 : ierr = pio_inq_dimid(fh, 'atmos_layer', did)
1855 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': atmos_layer not found')
1856 3072 : ierr = pio_inq_dimlen(fh, did, atmos_layer)
1857 :
1858 3072 : ierr = pio_inq_dimid(fh, 'bnd', did)
1859 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': bnd not found')
1860 3072 : ierr = pio_inq_dimlen(fh, did, bnd)
1861 :
1862 3072 : ierr = pio_inq_dimid(fh, 'pressure', did)
1863 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': pressure not found')
1864 3072 : ierr = pio_inq_dimlen(fh, did, pressure)
1865 :
1866 3072 : ierr = pio_inq_dimid(fh, 'temperature', did)
1867 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': temperature not found')
1868 3072 : ierr = pio_inq_dimlen(fh, did, temperature)
1869 :
1870 3072 : ierr = pio_inq_dimid(fh, 'absorber_ext', did)
1871 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': absorber_ext not found')
1872 3072 : ierr = pio_inq_dimlen(fh, did, absorber_ext)
1873 :
1874 3072 : ierr = pio_inq_dimid(fh, 'pressure_interp', did)
1875 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': pressure_interp not found')
1876 3072 : ierr = pio_inq_dimlen(fh, did, pressure_interp)
1877 :
1878 3072 : ierr = pio_inq_dimid(fh, 'mixing_fraction', did)
1879 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': mixing_fraction not found')
1880 3072 : ierr = pio_inq_dimlen(fh, did, mixing_fraction)
1881 :
1882 3072 : ierr = pio_inq_dimid(fh, 'gpt', did)
1883 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': gpt not found')
1884 3072 : ierr = pio_inq_dimlen(fh, did, gpt)
1885 :
1886 3072 : temperature_Planck = 0
1887 3072 : ierr = pio_inq_dimid(fh, 'temperature_Planck', did)
1888 3072 : if (ierr == PIO_NOERR) then
1889 1536 : ierr = pio_inq_dimlen(fh, did, temperature_Planck)
1890 : end if
1891 3072 : ierr = pio_inq_dimid(fh, 'pair', did)
1892 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': pair not found')
1893 3072 : ierr = pio_inq_dimlen(fh, did, pairs)
1894 3072 : ierr = pio_inq_dimid(fh, 'minor_absorber', did)
1895 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': minor_absorber not found')
1896 3072 : ierr = pio_inq_dimlen(fh, did, minorabsorbers)
1897 3072 : ierr = pio_inq_dimid(fh, 'minor_absorber_intervals_lower', did)
1898 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': minor_absorber_intervals_lower not found')
1899 3072 : ierr = pio_inq_dimlen(fh, did, minor_absorber_intervals_lower)
1900 3072 : ierr = pio_inq_dimid(fh, 'minor_absorber_intervals_upper', did)
1901 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': minor_absorber_intervals_upper not found')
1902 3072 : ierr = pio_inq_dimlen(fh, did, minor_absorber_intervals_upper)
1903 3072 : ierr = pio_inq_dimid(fh, 'contributors_lower', did)
1904 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': contributors_lower not found')
1905 3072 : ierr = pio_inq_dimlen(fh, did, contributors_lower)
1906 3072 : ierr = pio_inq_dimid(fh, 'contributors_upper', did)
1907 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': contributors_upper not found')
1908 3072 : ierr = pio_inq_dimlen(fh, did, contributors_upper)
1909 :
1910 3072 : ierr = pio_inq_dimid(fh, 'fit_coeffs', did)
1911 3072 : if (ierr == PIO_NOERR) then
1912 1536 : ierr = pio_inq_dimlen(fh, did, fit_coeffs)
1913 : end if
1914 :
1915 : ! Get variables
1916 :
1917 : ! names of absorbing gases
1918 9216 : allocate(gas_names(absorber), stat=istat)
1919 3072 : call handle_allocate_error(istat, sub, 'gas_names')
1920 3072 : ierr = pio_inq_varid(fh, 'gas_names', vid)
1921 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': gas_names not found')
1922 3072 : ierr = pio_get_var(fh, vid, gas_names)
1923 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading gas_names')
1924 :
1925 : ! key species pair for each band
1926 12288 : allocate(key_species(2,atmos_layer,bnd), stat=istat)
1927 3072 : call handle_allocate_error(istat, sub, 'key_species')
1928 3072 : ierr = pio_inq_varid(fh, 'key_species', vid)
1929 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': key_species not found')
1930 3072 : ierr = pio_get_var(fh, vid, key_species)
1931 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading key_species')
1932 :
1933 : ! beginning and ending gpoint for each band
1934 9216 : allocate(band2gpt(2,bnd), stat=istat)
1935 3072 : call handle_allocate_error(istat, sub, 'band2gpt')
1936 3072 : ierr = pio_inq_varid(fh, 'bnd_limits_gpt', vid)
1937 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': bnd_limits_gpt not found')
1938 3072 : ierr = pio_get_var(fh, vid, band2gpt)
1939 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading bnd_limits_gpt')
1940 :
1941 : ! beginning and ending wavenumber for each band
1942 9216 : allocate(band_lims_wavenum(2,bnd), stat=istat)
1943 3072 : call handle_allocate_error(istat, sub, 'band_lims_wavenum')
1944 3072 : ierr = pio_inq_varid(fh, 'bnd_limits_wavenumber', vid)
1945 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': bnd_limits_wavenumber not found')
1946 3072 : ierr = pio_get_var(fh, vid, band_lims_wavenum)
1947 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading bnd_limits_wavenumber')
1948 :
1949 : ! pressures [hPa] for reference atmosphere; press_ref(# reference layers)
1950 9216 : allocate(press_ref(pressure), stat=istat)
1951 3072 : call handle_allocate_error(istat, sub, 'press_ref')
1952 3072 : ierr = pio_inq_varid(fh, 'press_ref', vid)
1953 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': press_ref not found')
1954 3072 : ierr = pio_get_var(fh, vid, press_ref)
1955 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading press_ref')
1956 :
1957 : ! reference pressure separating the lower and upper atmosphere
1958 3072 : ierr = pio_inq_varid(fh, 'press_ref_trop', vid)
1959 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': press_ref_trop not found')
1960 3072 : ierr = pio_get_var(fh, vid, press_ref_trop)
1961 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading press_ref_trop')
1962 :
1963 : ! temperatures [K] for reference atmosphere; temp_ref(# reference layers)
1964 9216 : allocate(temp_ref(temperature), stat=istat)
1965 3072 : call handle_allocate_error(istat, sub, 'temp_ref')
1966 3072 : ierr = pio_inq_varid(fh, 'temp_ref', vid)
1967 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': temp_ref not found')
1968 3072 : ierr = pio_get_var(fh, vid, temp_ref)
1969 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading temp_ref')
1970 :
1971 : ! standard spectroscopic reference temperature [K]
1972 3072 : ierr = pio_inq_varid(fh, 'absorption_coefficient_ref_T', vid)
1973 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': absorption_coefficient_ref_T not found')
1974 3072 : ierr = pio_get_var(fh, vid, temp_ref_t)
1975 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading absorption_coefficient_ref_T')
1976 :
1977 : ! standard spectroscopic reference pressure [hPa]
1978 3072 : ierr = pio_inq_varid(fh, 'absorption_coefficient_ref_P', vid)
1979 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': absorption_coefficient_ref_P not found')
1980 3072 : ierr = pio_get_var(fh, vid, temp_ref_p)
1981 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading absorption_coefficient_ref_P')
1982 :
1983 : ! volume mixing ratios for reference atmosphere
1984 15360 : allocate(vmr_ref(atmos_layer, absorber_ext, temperature), stat=istat)
1985 3072 : call handle_allocate_error(istat, sub, 'vmr_ref')
1986 3072 : ierr = pio_inq_varid(fh, 'vmr_ref', vid)
1987 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': vmr_ref not found')
1988 3072 : ierr = pio_get_var(fh, vid, vmr_ref)
1989 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading vmr_ref')
1990 :
1991 : ! absorption coefficients due to major absorbing gases
1992 18432 : allocate(kmajor(gpt,mixing_fraction,pressure_interp,temperature), stat=istat)
1993 3072 : call handle_allocate_error(istat, sub, 'kmajor')
1994 3072 : ierr = pio_inq_varid(fh, 'kmajor', vid)
1995 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': kmajor not found')
1996 3072 : ierr = pio_get_var(fh, vid, kmajor)
1997 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading kmajor')
1998 :
1999 : ! absorption coefficients due to minor absorbing gases in lower part of atmosphere
2000 15360 : allocate(kminor_lower(contributors_lower, mixing_fraction, temperature), stat=istat)
2001 3072 : call handle_allocate_error(istat, sub, 'kminor_lower')
2002 3072 : ierr = pio_inq_varid(fh, 'kminor_lower', vid)
2003 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': kminor_lower not found')
2004 3072 : ierr = pio_get_var(fh, vid, kminor_lower)
2005 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading kminor_lower')
2006 :
2007 : ! absorption coefficients due to minor absorbing gases in upper part of atmosphere
2008 15360 : allocate(kminor_upper(contributors_upper, mixing_fraction, temperature), stat=istat)
2009 3072 : call handle_allocate_error(istat, sub, 'kminor_upper')
2010 3072 : ierr = pio_inq_varid(fh, 'kminor_upper', vid)
2011 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': kminor_upper not found')
2012 3072 : ierr = pio_get_var(fh, vid, kminor_upper)
2013 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading kminor_upper')
2014 :
2015 : ! integrated Planck function by band
2016 3072 : ierr = pio_inq_varid(fh, 'totplnk', vid)
2017 3072 : if (ierr == PIO_NOERR) then
2018 6144 : allocate(totplnk(temperature_Planck,bnd), stat=istat)
2019 1536 : call handle_allocate_error(istat, sub, 'totplnk')
2020 1536 : ierr = pio_get_var(fh, vid, totplnk)
2021 1536 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading totplnk')
2022 : end if
2023 :
2024 : ! Planck fractions
2025 3072 : ierr = pio_inq_varid(fh, 'plank_fraction', vid)
2026 3072 : if (ierr == PIO_NOERR) then
2027 9216 : allocate(planck_frac(gpt,mixing_fraction,pressure_interp,temperature), stat=istat)
2028 1536 : call handle_allocate_error(istat, sub, 'planck_frac')
2029 1536 : ierr = pio_get_var(fh, vid, planck_frac)
2030 1536 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading plank_fraction')
2031 : end if
2032 :
2033 3072 : ierr = pio_inq_varid(fh, 'optimal_angle_fit', vid)
2034 3072 : if (ierr == PIO_NOERR) then
2035 6144 : allocate(optimal_angle_fit(fit_coeffs, bnd), stat=istat)
2036 1536 : call handle_allocate_error(istat, sub, 'optiman_angle_fit')
2037 1536 : ierr = pio_get_var(fh, vid, optimal_angle_fit)
2038 1536 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading optimal_angle_fit')
2039 : end if
2040 :
2041 3072 : ierr = pio_inq_varid(fh, 'solar_source_quiet', vid)
2042 3072 : if (ierr == PIO_NOERR) then
2043 4608 : allocate(solar_src_quiet(gpt), stat=istat)
2044 1536 : call handle_allocate_error(istat, sub, 'solar_src_quiet')
2045 1536 : ierr = pio_get_var(fh, vid, solar_src_quiet)
2046 1536 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading solar_source_quiet')
2047 : end if
2048 :
2049 3072 : ierr = pio_inq_varid(fh, 'solar_source_facular', vid)
2050 3072 : if (ierr == PIO_NOERR) then
2051 4608 : allocate(solar_src_facular(gpt), stat=istat)
2052 1536 : call handle_allocate_error(istat, sub, 'solar_src_facular')
2053 1536 : ierr = pio_get_var(fh, vid, solar_src_facular)
2054 1536 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading solar_source_facular')
2055 : end if
2056 :
2057 3072 : ierr = pio_inq_varid(fh, 'solar_source_sunspot', vid)
2058 3072 : if (ierr == PIO_NOERR) then
2059 4608 : allocate(solar_src_sunspot(gpt), stat=istat)
2060 1536 : call handle_allocate_error(istat, sub, 'solar_src_sunspot')
2061 1536 : ierr = pio_get_var(fh, vid, solar_src_sunspot)
2062 1536 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading solar_source_sunspot')
2063 : end if
2064 :
2065 3072 : ierr = pio_inq_varid(fh, 'tsi_default', vid)
2066 3072 : if (ierr == PIO_NOERR) then
2067 1536 : ierr = pio_get_var(fh, vid, tsi_default)
2068 1536 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading tsi_default')
2069 : end if
2070 :
2071 3072 : ierr = pio_inq_varid(fh, 'mg_default', vid)
2072 3072 : if (ierr == PIO_NOERR) then
2073 1536 : ierr = pio_get_var(fh, vid, mg_default)
2074 1536 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading mg_default')
2075 : end if
2076 :
2077 3072 : ierr = pio_inq_varid(fh, 'sb_default', vid)
2078 3072 : if (ierr == PIO_NOERR) then
2079 1536 : ierr = pio_get_var(fh, vid, sb_default)
2080 1536 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading sb_default')
2081 : end if
2082 :
2083 : ! rayleigh scattering contribution in lower part of atmosphere
2084 3072 : ierr = pio_inq_varid(fh, 'rayl_lower', vid)
2085 3072 : if (ierr == PIO_NOERR) then
2086 7680 : allocate(rayl_lower(gpt,mixing_fraction,temperature), stat=istat)
2087 1536 : call handle_allocate_error(istat, sub, 'rayl_lower')
2088 1536 : ierr = pio_get_var(fh, vid, rayl_lower)
2089 1536 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading rayl_lower')
2090 : end if
2091 :
2092 : ! rayleigh scattering contribution in upper part of atmosphere
2093 3072 : ierr = pio_inq_varid(fh, 'rayl_upper', vid)
2094 3072 : if (ierr == PIO_NOERR) then
2095 7680 : allocate(rayl_upper(gpt,mixing_fraction,temperature), stat=istat)
2096 1536 : call handle_allocate_error(istat, sub, 'rayl_upper')
2097 1536 : ierr = pio_get_var(fh, vid, rayl_upper)
2098 1536 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading rayl_upper')
2099 : end if
2100 :
2101 9216 : allocate(gas_minor(minorabsorbers), stat=istat)
2102 3072 : call handle_allocate_error(istat, sub, 'gas_minor')
2103 3072 : ierr = pio_inq_varid(fh, 'gas_minor', vid)
2104 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': gas_minor not found')
2105 3072 : ierr = pio_get_var(fh, vid, gas_minor)
2106 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading gas_minor')
2107 :
2108 9216 : allocate(identifier_minor(minorabsorbers), stat=istat)
2109 3072 : call handle_allocate_error(istat, sub, 'identifier_minor')
2110 3072 : ierr = pio_inq_varid(fh, 'identifier_minor', vid)
2111 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': identifier_minor not found')
2112 3072 : ierr = pio_get_var(fh, vid, identifier_minor)
2113 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading identifier_minor')
2114 :
2115 9216 : allocate(minor_gases_lower(minor_absorber_intervals_lower), stat=istat)
2116 3072 : call handle_allocate_error(istat, sub, 'minor_gases_lower')
2117 3072 : ierr = pio_inq_varid(fh, 'minor_gases_lower', vid)
2118 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': minor_gases_lower not found')
2119 3072 : ierr = pio_get_var(fh, vid, minor_gases_lower)
2120 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading minor_gases_lower')
2121 :
2122 9216 : allocate(minor_gases_upper(minor_absorber_intervals_upper), stat=istat)
2123 3072 : call handle_allocate_error(istat, sub, 'minor_gases_upper')
2124 3072 : ierr = pio_inq_varid(fh, 'minor_gases_upper', vid)
2125 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': minor_gases_upper not found')
2126 3072 : ierr = pio_get_var(fh, vid, minor_gases_upper)
2127 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading minor_gases_upper')
2128 :
2129 12288 : allocate(minor_limits_gpt_lower(pairs,minor_absorber_intervals_lower), stat=istat)
2130 3072 : call handle_allocate_error(istat, sub, 'minor_limits_gpt_lower')
2131 3072 : ierr = pio_inq_varid(fh, 'minor_limits_gpt_lower', vid)
2132 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': minor_limits_gpt_lower not found')
2133 3072 : ierr = pio_get_var(fh, vid, minor_limits_gpt_lower)
2134 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading minor_limits_gpt_lower')
2135 :
2136 12288 : allocate(minor_limits_gpt_upper(pairs,minor_absorber_intervals_upper), stat=istat)
2137 3072 : call handle_allocate_error(istat, sub, 'minor_limits_gpt_upper')
2138 3072 : ierr = pio_inq_varid(fh, 'minor_limits_gpt_upper', vid)
2139 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': minor_limits_gpt_upper not found')
2140 3072 : ierr = pio_get_var(fh, vid, minor_limits_gpt_upper)
2141 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading minor_limits_gpt_upper')
2142 :
2143 : ! Read as integer and convert to logical
2144 9216 : allocate(int2log(minor_absorber_intervals_lower), stat=istat)
2145 3072 : call handle_allocate_error(istat, sub, 'int2log for lower')
2146 :
2147 9216 : allocate(minor_scales_with_density_lower(minor_absorber_intervals_lower), stat=istat)
2148 3072 : call handle_allocate_error(istat, sub, 'minor_scales_with_density_lower')
2149 3072 : ierr = pio_inq_varid(fh, 'minor_scales_with_density_lower', vid)
2150 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': minor_scales_with_density_lower not found')
2151 3072 : ierr = pio_get_var(fh, vid, int2log)
2152 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading minor_scales_with_density_lower')
2153 147456 : do i = 1,minor_absorber_intervals_lower
2154 147456 : if (int2log(i) .eq. 0) then
2155 53760 : minor_scales_with_density_lower(i) = .false.
2156 : else
2157 90624 : minor_scales_with_density_lower(i) = .true.
2158 : end if
2159 : end do
2160 :
2161 9216 : allocate(scale_by_complement_lower(minor_absorber_intervals_lower), stat=istat)
2162 3072 : call handle_allocate_error(istat, sub, 'scale_by_complement_lower')
2163 3072 : ierr = pio_inq_varid(fh, 'scale_by_complement_lower', vid)
2164 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': scale_by_complement_lower not found')
2165 3072 : ierr = pio_get_var(fh, vid, int2log)
2166 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading scale_by_complement_lower')
2167 147456 : do i = 1,minor_absorber_intervals_lower
2168 147456 : if (int2log(i) .eq. 0) then
2169 104448 : scale_by_complement_lower(i) = .false.
2170 : else
2171 39936 : scale_by_complement_lower(i) = .true.
2172 : end if
2173 : end do
2174 :
2175 3072 : deallocate(int2log)
2176 :
2177 : ! Read as integer and convert to logical
2178 9216 : allocate(int2log(minor_absorber_intervals_upper), stat=istat)
2179 3072 : call handle_allocate_error(istat, sub, 'int2log for upper')
2180 :
2181 9216 : allocate(minor_scales_with_density_upper(minor_absorber_intervals_upper), stat=istat)
2182 3072 : call handle_allocate_error(istat, sub, 'minor_scales_with_density_upper')
2183 3072 : ierr = pio_inq_varid(fh, 'minor_scales_with_density_upper', vid)
2184 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': minor_scales_with_density_upper not found')
2185 3072 : ierr = pio_get_var(fh, vid, int2log)
2186 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading minor_scales_with_density_upper')
2187 92160 : do i = 1,minor_absorber_intervals_upper
2188 92160 : if (int2log(i) .eq. 0) then
2189 46080 : minor_scales_with_density_upper(i) = .false.
2190 : else
2191 43008 : minor_scales_with_density_upper(i) = .true.
2192 : end if
2193 : end do
2194 :
2195 9216 : allocate(scale_by_complement_upper(minor_absorber_intervals_upper), stat=istat)
2196 3072 : call handle_allocate_error(istat, sub, 'scale_by_complement_upper')
2197 3072 : ierr = pio_inq_varid(fh, 'scale_by_complement_upper', vid)
2198 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': scale_by_complement_upper not found')
2199 3072 : ierr = pio_get_var(fh, vid, int2log)
2200 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading scale_by_complement_upper')
2201 92160 : do i = 1,minor_absorber_intervals_upper
2202 92160 : if (int2log(i) .eq. 0) then
2203 70656 : scale_by_complement_upper(i) = .false.
2204 : else
2205 18432 : scale_by_complement_upper(i) = .true.
2206 : end if
2207 : end do
2208 :
2209 3072 : deallocate(int2log)
2210 :
2211 9216 : allocate(scaling_gas_lower(minor_absorber_intervals_lower), stat=istat)
2212 3072 : call handle_allocate_error(istat, sub, 'scaling_gas_lower')
2213 3072 : ierr = pio_inq_varid(fh, 'scaling_gas_lower', vid)
2214 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': scaling_gas_lower not found')
2215 3072 : ierr = pio_get_var(fh, vid, scaling_gas_lower)
2216 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading scaling_gas_lower')
2217 :
2218 9216 : allocate(scaling_gas_upper(minor_absorber_intervals_upper), stat=istat)
2219 3072 : call handle_allocate_error(istat, sub, 'scaling_gas_upper')
2220 3072 : ierr = pio_inq_varid(fh, 'scaling_gas_upper', vid)
2221 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': scaling_gas_upper not found')
2222 3072 : ierr = pio_get_var(fh, vid, scaling_gas_upper)
2223 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading scaling_gas_upper')
2224 :
2225 9216 : allocate(kminor_start_lower(minor_absorber_intervals_lower), stat=istat)
2226 3072 : call handle_allocate_error(istat, sub, 'kminor_start_lower')
2227 3072 : ierr = pio_inq_varid(fh, 'kminor_start_lower', vid)
2228 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': kminor_start_lower not found')
2229 3072 : ierr = pio_get_var(fh, vid, kminor_start_lower)
2230 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading kminor_start_lower')
2231 :
2232 9216 : allocate(kminor_start_upper(minor_absorber_intervals_upper), stat=istat)
2233 3072 : call handle_allocate_error(istat, sub, 'kminor_start_upper')
2234 3072 : ierr = pio_inq_varid(fh, 'kminor_start_upper', vid)
2235 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': kminor_start_upper not found')
2236 3072 : ierr = pio_get_var(fh, vid, kminor_start_upper)
2237 3072 : if (ierr /= PIO_NOERR) call endrun(sub//': error reading kminor_start_upper')
2238 :
2239 : ! Close file
2240 3072 : call pio_closefile(fh)
2241 :
2242 : ! Initialize the gas optics object with data. The calls are slightly different depending
2243 : ! on whether the radiation sources are internal to the atmosphere (longwave) or external (shortwave)
2244 :
2245 3072 : if (allocated(totplnk) .and. allocated(planck_frac)) then
2246 : error_msg = kdist%load( &
2247 : available_gases, gas_names, key_species, &
2248 : band2gpt, band_lims_wavenum, &
2249 : press_ref, press_ref_trop, temp_ref, &
2250 : temp_ref_p, temp_ref_t, vmr_ref, &
2251 : kmajor, kminor_lower, kminor_upper, &
2252 : gas_minor, identifier_minor, &
2253 : minor_gases_lower, minor_gases_upper, &
2254 : minor_limits_gpt_lower, minor_limits_gpt_upper, &
2255 : minor_scales_with_density_lower, &
2256 : minor_scales_with_density_upper, &
2257 : scaling_gas_lower, scaling_gas_upper, &
2258 : scale_by_complement_lower, scale_by_complement_upper, &
2259 : kminor_start_lower, kminor_start_upper, &
2260 : totplnk, planck_frac, rayl_lower, rayl_upper, &
2261 1536 : optimal_angle_fit)
2262 1536 : else if (allocated(solar_src_quiet)) then
2263 : error_msg = kdist%load( &
2264 : available_gases, gas_names, key_species, &
2265 : band2gpt, band_lims_wavenum, &
2266 : press_ref, press_ref_trop, temp_ref, &
2267 : temp_ref_p, temp_ref_t, vmr_ref, &
2268 : kmajor, kminor_lower, kminor_upper, &
2269 : gas_minor, identifier_minor, &
2270 : minor_gases_lower, minor_gases_upper, &
2271 : minor_limits_gpt_lower, minor_limits_gpt_upper, &
2272 : minor_scales_with_density_lower, &
2273 : minor_scales_with_density_upper, &
2274 : scaling_gas_lower, scaling_gas_upper, &
2275 : scale_by_complement_lower, &
2276 : scale_by_complement_upper, &
2277 : kminor_start_lower, &
2278 : kminor_start_upper, &
2279 : solar_src_quiet, solar_src_facular, solar_src_sunspot, &
2280 : tsi_default, mg_default, sb_default, &
2281 1536 : rayl_lower, rayl_upper)
2282 : else
2283 0 : error_msg = 'must supply either totplnk and planck_frac, or solar_src_[*]'
2284 : end if
2285 :
2286 3072 : call stop_on_err(error_msg, sub, 'kdist%load')
2287 :
2288 0 : deallocate( &
2289 0 : gas_names, key_species, &
2290 0 : band2gpt, band_lims_wavenum, &
2291 0 : press_ref, temp_ref, vmr_ref, &
2292 0 : kmajor, kminor_lower, kminor_upper, &
2293 0 : gas_minor, identifier_minor, &
2294 0 : minor_gases_lower, minor_gases_upper, &
2295 0 : minor_limits_gpt_lower, &
2296 0 : minor_limits_gpt_upper, &
2297 0 : minor_scales_with_density_lower, &
2298 0 : minor_scales_with_density_upper, &
2299 0 : scale_by_complement_lower, &
2300 0 : scale_by_complement_upper, &
2301 0 : scaling_gas_lower, scaling_gas_upper, &
2302 3072 : kminor_start_lower, kminor_start_upper)
2303 :
2304 3072 : if (allocated(totplnk)) deallocate(totplnk)
2305 3072 : if (allocated(planck_frac)) deallocate(planck_frac)
2306 3072 : if (allocated(optimal_angle_fit)) deallocate(optimal_angle_fit)
2307 3072 : if (allocated(solar_src_quiet)) deallocate(solar_src_quiet)
2308 3072 : if (allocated(solar_src_facular)) deallocate(solar_src_facular)
2309 3072 : if (allocated(solar_src_sunspot)) deallocate(solar_src_sunspot)
2310 3072 : if (allocated(rayl_lower)) deallocate(rayl_lower)
2311 3072 : if (allocated(rayl_upper)) deallocate(rayl_upper)
2312 :
2313 6144 : end subroutine coefs_init
2314 :
2315 : !=========================================================================================
2316 :
2317 5969088 : subroutine initialize_rrtmgp_fluxes(ncol, nlevels, nbands, fluxes, do_direct)
2318 :
2319 : ! Allocate flux arrays and set values to zero.
2320 :
2321 : ! Arguments
2322 : integer, intent(in) :: ncol, nlevels, nbands
2323 : class(ty_fluxes_broadband), intent(inout) :: fluxes
2324 : logical, optional, intent(in) :: do_direct
2325 :
2326 : ! Local variables
2327 : logical :: do_direct_local
2328 : integer :: istat
2329 : character(len=*), parameter :: sub = 'initialize_rrtmgp_fluxes'
2330 : !----------------------------------------------------------------------------
2331 :
2332 5969088 : if (present(do_direct)) then
2333 : do_direct_local = .true.
2334 : else
2335 2984544 : do_direct_local = .false.
2336 : end if
2337 :
2338 : ! Broadband fluxes
2339 22448964 : allocate(fluxes%flux_up(ncol, nlevels), stat=istat)
2340 5969088 : call handle_allocate_error(istat, sub, 'fluxes%flux_up')
2341 16479876 : allocate(fluxes%flux_dn(ncol, nlevels), stat=istat)
2342 5969088 : call handle_allocate_error(istat, sub, 'fluxes%flux_dn')
2343 16479876 : allocate(fluxes%flux_net(ncol, nlevels), stat=istat)
2344 5969088 : call handle_allocate_error(istat, sub, 'fluxes%flux_net')
2345 5969088 : if (do_direct_local) then
2346 7526244 : allocate(fluxes%flux_dn_dir(ncol, nlevels), stat=istat)
2347 2984544 : call handle_allocate_error(istat, sub, 'fluxes%flux_dn_dir')
2348 : end if
2349 :
2350 : select type (fluxes)
2351 : type is (ty_fluxes_byband)
2352 : ! Fluxes by band always needed for SW. Only allocate for LW
2353 : ! when spectralflux is true.
2354 2984544 : if (nbands == nswbands .or. spectralflux) then
2355 6747666 : allocate(fluxes%bnd_flux_up(ncol, nlevels, nbands), stat=istat)
2356 1492272 : call handle_allocate_error(istat, sub, 'fluxes%bnd_flux_up')
2357 5255394 : allocate(fluxes%bnd_flux_dn(ncol, nlevels, nbands), stat=istat)
2358 1492272 : call handle_allocate_error(istat, sub, 'fluxes%bnd_flux_dn')
2359 5255394 : allocate(fluxes%bnd_flux_net(ncol, nlevels, nbands), stat=istat)
2360 1492272 : call handle_allocate_error(istat, sub, 'fluxes%bnd_flux_net')
2361 1492272 : if (do_direct_local) then
2362 5255394 : allocate(fluxes%bnd_flux_dn_dir(ncol, nlevels, nbands), stat=istat)
2363 1492272 : call handle_allocate_error(istat, sub, 'fluxes%bnd_flux_dn_dir')
2364 : end if
2365 : end if
2366 : end select
2367 :
2368 : ! Initialize
2369 5969088 : call reset_fluxes(fluxes)
2370 :
2371 5969088 : end subroutine initialize_rrtmgp_fluxes
2372 :
2373 : !=========================================================================================
2374 :
2375 5969088 : subroutine reset_fluxes(fluxes)
2376 :
2377 : ! Reset flux arrays to zero.
2378 :
2379 : class(ty_fluxes_broadband), intent(inout) :: fluxes
2380 : !----------------------------------------------------------------------------
2381 :
2382 : ! Reset broadband fluxes
2383 7249214448 : fluxes%flux_up(:,:) = 0._r8
2384 7249214448 : fluxes%flux_dn(:,:) = 0._r8
2385 7249214448 : fluxes%flux_net(:,:) = 0._r8
2386 2514894768 : if (associated(fluxes%flux_dn_dir)) fluxes%flux_dn_dir(:,:) = 0._r8
2387 :
2388 : select type (fluxes)
2389 : type is (ty_fluxes_byband)
2390 : ! Reset band-by-band fluxes
2391 17584863840 : if (associated(fluxes%bnd_flux_up)) fluxes%bnd_flux_up(:,:,:) = 0._r8
2392 17586356112 : if (associated(fluxes%bnd_flux_dn)) fluxes%bnd_flux_dn(:,:,:) = 0._r8
2393 17586356112 : if (associated(fluxes%bnd_flux_net)) fluxes%bnd_flux_net(:,:,:) = 0._r8
2394 17589340656 : if (associated(fluxes%bnd_flux_dn_dir)) fluxes%bnd_flux_dn_dir(:,:,:) = 0._r8
2395 : end select
2396 :
2397 5969088 : end subroutine reset_fluxes
2398 :
2399 : !=========================================================================================
2400 :
2401 4476816 : subroutine free_optics_sw(optics)
2402 :
2403 : type(ty_optical_props_2str), intent(inout) :: optics
2404 :
2405 4476816 : if (allocated(optics%tau)) deallocate(optics%tau)
2406 4476816 : if (allocated(optics%ssa)) deallocate(optics%ssa)
2407 4476816 : if (allocated(optics%g)) deallocate(optics%g)
2408 4476816 : call optics%finalize()
2409 :
2410 4476816 : end subroutine free_optics_sw
2411 :
2412 : !=========================================================================================
2413 :
2414 2984544 : subroutine free_optics_lw(optics)
2415 :
2416 : type(ty_optical_props_1scl), intent(inout) :: optics
2417 :
2418 2984544 : if (allocated(optics%tau)) deallocate(optics%tau)
2419 2984544 : call optics%finalize()
2420 :
2421 2984544 : end subroutine free_optics_lw
2422 :
2423 : !=========================================================================================
2424 :
2425 5969088 : subroutine free_fluxes(fluxes)
2426 :
2427 : class(ty_fluxes_broadband), intent(inout) :: fluxes
2428 :
2429 5969088 : if (associated(fluxes%flux_up)) deallocate(fluxes%flux_up)
2430 5969088 : if (associated(fluxes%flux_dn)) deallocate(fluxes%flux_dn)
2431 5969088 : if (associated(fluxes%flux_net)) deallocate(fluxes%flux_net)
2432 5969088 : if (associated(fluxes%flux_dn_dir)) deallocate(fluxes%flux_dn_dir)
2433 :
2434 : select type (fluxes)
2435 : type is (ty_fluxes_byband)
2436 1492272 : if (associated(fluxes%bnd_flux_up)) deallocate(fluxes%bnd_flux_up)
2437 2984544 : if (associated(fluxes%bnd_flux_dn)) deallocate(fluxes%bnd_flux_dn)
2438 2984544 : if (associated(fluxes%bnd_flux_net)) deallocate(fluxes%bnd_flux_net)
2439 5969088 : if (associated(fluxes%bnd_flux_dn_dir)) deallocate(fluxes%bnd_flux_dn_dir)
2440 : end select
2441 :
2442 5969088 : end subroutine free_fluxes
2443 :
2444 : !=========================================================================================
2445 :
2446 746136 : subroutine modified_cloud_fraction(ncol, cld, cldfsnow, cldfgrau, cldfprime)
2447 :
2448 : ! Compute modified cloud fraction, cldfprime.
2449 : ! 1. initialize as cld
2450 : ! 2. modify for snow if cldfsnow is available. use max(cld, cldfsnow)
2451 : ! 3. modify for graupel if cldfgrau is available and graupel_in_rad is true.
2452 : ! use max(cldfprime, cldfgrau)
2453 :
2454 : ! Arguments
2455 : integer, intent(in) :: ncol
2456 : real(r8), pointer :: cld(:,:) ! cloud fraction
2457 : real(r8), pointer :: cldfsnow(:,:) ! cloud fraction of just "snow clouds"
2458 : real(r8), pointer :: cldfgrau(:,:) ! cloud fraction of just "graupel clouds"
2459 : real(r8), intent(out) :: cldfprime(:,:) ! modified cloud fraction
2460 :
2461 : ! Local variables
2462 : integer :: i, k
2463 : !----------------------------------------------------------------------------
2464 :
2465 746136 : if (associated(cldfsnow)) then
2466 70136784 : do k = 1, pver
2467 1159408584 : do i = 1, ncol
2468 1158662448 : cldfprime(i,k) = max(cld(i,k), cldfsnow(i,k))
2469 : end do
2470 : end do
2471 : else
2472 0 : cldfprime(:ncol,:) = cld(:ncol,:)
2473 : end if
2474 :
2475 746136 : if (associated(cldfgrau) .and. graupel_in_rad) then
2476 0 : do k = 1, pver
2477 0 : do i = 1, ncol
2478 0 : cldfprime(i,k) = max(cldfprime(i,k), cldfgrau(i,k))
2479 : end do
2480 : end do
2481 : end if
2482 :
2483 746136 : end subroutine modified_cloud_fraction
2484 :
2485 : !=========================================================================================
2486 :
2487 9833936 : subroutine stop_on_err(errmsg, sub, info)
2488 :
2489 : ! call endrun if RRTMGP function returns non-empty error message.
2490 :
2491 : character(len=*), intent(in) :: errmsg ! return message from RRTMGP function
2492 : character(len=*), intent(in) :: sub ! name of calling subroutine
2493 : character(len=*), intent(in) :: info ! name of called function
2494 :
2495 9833936 : if (len_trim(errmsg) > 0) then
2496 0 : call endrun(trim(sub)//': ERROR: '//trim(info)//': '//trim(errmsg))
2497 : end if
2498 :
2499 9833936 : end subroutine stop_on_err
2500 :
2501 : !=========================================================================================
2502 :
2503 17907264 : end module radiation
2504 :
|
