Line data Source code
1 : module rrtmgp_inputs
2 :
3 : !--------------------------------------------------------------------------------
4 : ! Transform data for inputs from CAM's data structures to those used by
5 : ! RRTMGP. Subset the number of model levels if CAM's top exceeds RRTMGP's
6 : ! valid domain. Add an extra layer if CAM's top is below 1 Pa.
7 : ! The vertical indexing increases from top to bottom of atmosphere in both
8 : ! CAM and RRTMGP arrays.
9 : !--------------------------------------------------------------------------------
10 :
11 : use shr_kind_mod, only: r8=>shr_kind_r8
12 : use ppgrid, only: pcols, pver, pverp
13 :
14 : use physconst, only: stebol, pi
15 :
16 : use physics_types, only: physics_state
17 : use physics_buffer, only: physics_buffer_desc
18 : use camsrfexch, only: cam_in_t
19 :
20 : use radconstants, only: nradgas, gaslist, nswbands, nlwbands, nswgpts, nlwgpts, &
21 : get_sw_spectral_boundaries, idx_sw_diag, idx_sw_cloudsim, &
22 : idx_lw_cloudsim
23 :
24 : use rad_constituents, only: rad_cnst_get_gas
25 :
26 : use cloud_rad_props, only: get_liquid_optics_sw, liquid_cloud_get_rad_props_lw, &
27 : get_ice_optics_sw, ice_cloud_get_rad_props_lw, &
28 : get_snow_optics_sw, snow_cloud_get_rad_props_lw, &
29 : get_grau_optics_sw, grau_cloud_get_rad_props_lw
30 :
31 : use mcica_subcol_gen, only: mcica_subcol_sw, mcica_subcol_lw
32 :
33 : use aer_rad_props, only: aer_rad_props_sw, aer_rad_props_lw
34 :
35 : use mo_gas_concentrations, only: ty_gas_concs
36 : use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp
37 : use mo_optical_props, only: ty_optical_props_2str, ty_optical_props_1scl
38 :
39 : use cam_history_support, only: fillvalue
40 : use cam_logfile, only: iulog
41 : use cam_abortutils, only: endrun
42 : use error_messages, only: alloc_err
43 :
44 : implicit none
45 : private
46 : save
47 :
48 : public :: &
49 : rrtmgp_inputs_init, &
50 : rrtmgp_set_state, &
51 : rrtmgp_set_gases_lw, &
52 : rrtmgp_set_gases_sw, &
53 : rrtmgp_set_cloud_lw, &
54 : rrtmgp_set_cloud_sw, &
55 : rrtmgp_set_aer_lw, &
56 : rrtmgp_set_aer_sw
57 :
58 :
59 : ! This value is to match the arbitrary small value used in RRTMG to decide
60 : ! when a quantity is effectively zero.
61 : real(r8), parameter :: tiny = 1.0e-80_r8
62 :
63 : ! Indices for copying data between cam and rrtmgp arrays
64 : integer :: ktopcam ! Index in CAM arrays of top level (layer or interface) at which
65 : ! RRTMGP is active.
66 : integer :: ktoprad ! Index in RRTMGP arrays of the layer or interface corresponding
67 : ! to CAM's top layer or interface
68 :
69 : ! wavenumber (cm^-1) boundaries of shortwave bands
70 : real(r8) :: sw_low_bounds(nswbands), sw_high_bounds(nswbands)
71 :
72 : ! Mapping from RRTMG shortwave bands to RRTMGP. Currently needed to continue using
73 : ! the SW optics datasets from RRTMG (even thought there is a slight mismatch in the
74 : ! band boundaries of the 2 bands that overlap with the LW bands).
75 : integer, parameter, dimension(14) :: rrtmg_to_rrtmgp_swbands = &
76 : [ 14, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 ]
77 :
78 : !==================================================================================================
79 : contains
80 : !==================================================================================================
81 :
82 1536 : subroutine rrtmgp_inputs_init(ktcam, ktrad)
83 :
84 : ! Note that this routine must be called after the calls to set_wavenumber_bands which set
85 : ! the sw/lw band boundaries in the radconstants module.
86 :
87 : integer, intent(in) :: ktcam
88 : integer, intent(in) :: ktrad
89 :
90 1536 : ktopcam = ktcam
91 1536 : ktoprad = ktrad
92 :
93 : ! Initialize the module data containing the SW band boundaries.
94 1536 : call get_sw_spectral_boundaries(sw_low_bounds, sw_high_bounds, 'cm^-1')
95 :
96 1536 : end subroutine rrtmgp_inputs_init
97 :
98 : !=========================================================================================
99 :
100 746136 : subroutine rrtmgp_set_state( &
101 : state, cam_in, ncol, nlay, nday, &
102 2238408 : idxday, coszrs, kdist_sw, t_sfc, emis_sfc, &
103 746136 : t_rad, pmid_rad, pint_rad, t_day, pmid_day, &
104 746136 : pint_day, coszrs_day, alb_dir, alb_dif)
105 :
106 : ! arguments
107 : type(physics_state), intent(in) :: state ! CAM physics state
108 : type(cam_in_t), intent(in) :: cam_in ! CAM import state
109 : integer, intent(in) :: ncol ! # cols in CAM chunk
110 : integer, intent(in) :: nlay ! # layers in rrtmgp grid
111 : integer, intent(in) :: nday ! # daylight columns
112 : integer, intent(in) :: idxday(:) ! chunk indicies of daylight columns
113 : real(r8), intent(in) :: coszrs(:) ! cosine of solar zenith angle
114 : class(ty_gas_optics_rrtmgp), intent(in) :: kdist_sw ! spectral information
115 :
116 : real(r8), intent(out) :: t_sfc(ncol) ! surface temperature [K]
117 : real(r8), intent(out) :: emis_sfc(nlwbands,ncol) ! emissivity at surface []
118 : real(r8), intent(out) :: t_rad(ncol,nlay) ! layer midpoint temperatures [K]
119 : real(r8), intent(out) :: pmid_rad(ncol,nlay) ! layer midpoint pressures [Pa]
120 : real(r8), intent(out) :: pint_rad(ncol,nlay+1) ! layer interface pressures [Pa]
121 : real(r8), intent(out) :: t_day(nday,nlay) ! layer midpoint temperatures [K]
122 : real(r8), intent(out) :: pmid_day(nday,nlay) ! layer midpoint pressure [Pa]
123 : real(r8), intent(out) :: pint_day(nday,nlay+1) ! layer interface pressures [Pa]
124 : real(r8), intent(out) :: coszrs_day(nday) ! cosine of solar zenith angle
125 : real(r8), intent(out) :: alb_dir(nswbands,nday) ! surface albedo, direct radiation
126 : real(r8), intent(out) :: alb_dif(nswbands,nday) ! surface albedo, diffuse radiation
127 :
128 : ! local variables
129 : integer :: i, k, iband
130 :
131 : real(r8) :: tref_min, tref_max
132 :
133 : character(len=*), parameter :: sub='rrtmgp_set_state'
134 : character(len=512) :: errmsg
135 : !--------------------------------------------------------------------------------
136 :
137 12458736 : t_sfc = sqrt(sqrt(cam_in%lwup(:ncol)/stebol)) ! Surface temp set based on longwave up flux.
138 :
139 : ! Set surface emissivity to 1.0.
140 : ! The land model *does* have its own surface emissivity, but is not spectrally resolved.
141 : ! The LW upward flux is calculated with that land emissivity, and the "radiative temperature"
142 : ! t_sfc is derived from that flux. We assume, therefore, that the emissivity is unity
143 : ! to be consistent with t_sfc.
144 199860336 : emis_sfc(:,:) = 1._r8
145 :
146 : ! Level ordering is the same for both CAM and RRTMGP (top to bottom)
147 1159408584 : t_rad(:,ktoprad:) = state%t(:ncol,ktopcam:)
148 1159408584 : pmid_rad(:,ktoprad:) = state%pmid(:ncol,ktopcam:)
149 1171867320 : pint_rad(:,ktoprad:) = state%pint(:ncol,ktopcam:)
150 :
151 : ! Add extra layer values if needed.
152 746136 : if (nlay == pverp) then
153 12458736 : t_rad(:,1) = state%t(:ncol,1)
154 : ! The top reference pressure from the RRTMGP coefficients datasets is 1.005183574463 Pa
155 : ! Set the top of the extra layer just below that.
156 12458736 : pint_rad(:,1) = 1.01_r8
157 :
158 : ! next interface down in LT will always be > 1Pa
159 : ! but in MT we apply adjustment to have it be 1.02 Pa if it was too high
160 12458736 : where (pint_rad(:,2) <= pint_rad(:,1)) pint_rad(:,2) = pint_rad(:,1)+0.01_r8
161 :
162 : ! set the highest pmid (in the "extra layer") to the midpoint (guarantees > 1Pa)
163 12458736 : pmid_rad(:,1) = pint_rad(:,1) + 0.5_r8 * (pint_rad(:,2) - pint_rad(:,1))
164 :
165 : ! For case of CAM MT, also ensure pint_rad(:,2) > pint_rad(:,1) & pmid_rad(:,2) > max(pmid_rad(:,1), min_pressure)
166 12458736 : where (pmid_rad(:,2) <= kdist_sw%get_press_min()) pmid_rad(:,2) = pint_rad(:,2) + 0.01_r8
167 : else
168 : ! nlay < pverp, thus the 1 Pa level is within a CAM layer. Assuming the top interface of
169 : ! this layer is at a pressure < 1 Pa, we need to adjust the top of this layer so that it
170 : ! is within the valid pressure range of RRTMGP (otherwise RRTMGP issues an error). Then
171 : ! set the midpoint pressure halfway between the interfaces.
172 0 : pint_rad(:,1) = 1.01_r8
173 0 : pmid_rad(:,1) = 0.5_r8 * (pint_rad(:,1) + pint_rad(:,2))
174 : end if
175 :
176 : ! Limit temperatures to be within the limits of RRTMGP validity.
177 746136 : tref_min = kdist_sw%get_temp_min()
178 746136 : tref_max = kdist_sw%get_temp_max()
179 1171867320 : t_rad = merge(t_rad, tref_min, t_rad > tref_min)
180 1171867320 : t_rad = merge(t_rad, tref_max, t_rad < tref_max)
181 12458736 : t_sfc = merge(t_sfc, tref_min, t_sfc > tref_min)
182 12458736 : t_sfc = merge(t_sfc, tref_max, t_sfc < tref_max)
183 :
184 : ! Construct arrays containing only daylight columns
185 6602436 : do i = 1, nday
186 556348500 : t_day(i,:) = t_rad(idxday(i),:)
187 556348500 : pmid_day(i,:) = pmid_rad(idxday(i),:)
188 562204800 : pint_day(i,:) = pint_rad(idxday(i),:)
189 6602436 : coszrs_day(i) = coszrs(idxday(i))
190 : end do
191 :
192 : ! Assign albedos to the daylight columns (from E3SM implementation)
193 : ! Albedos are imported from the surface models as broadband (visible, and near-IR),
194 : ! and we need to map these to appropriate narrower bands used in RRTMGP. Bands
195 : ! are categorized broadly as "visible/UV" or "infrared" based on wavenumber.
196 : ! Loop over bands, and determine for each band whether it is broadly in the
197 : ! visible or infrared part of the spectrum based on a dividing line of
198 : ! 0.7 micron, or 14286 cm^-1
199 11192040 : do iband = 1,nswbands
200 10445904 : if (is_visible(sw_low_bounds(iband)) .and. &
201 746136 : is_visible(sw_high_bounds(iband))) then
202 :
203 : ! Entire band is in the visible
204 26409744 : do i = 1, nday
205 23425200 : alb_dir(iband,i) = cam_in%asdir(idxday(i))
206 26409744 : alb_dif(iband,i) = cam_in%asdif(idxday(i))
207 : end do
208 :
209 7461360 : else if (.not.is_visible(sw_low_bounds(iband)) .and. &
210 : .not.is_visible(sw_high_bounds(iband))) then
211 : ! Entire band is in the longwave (near-infrared)
212 59421924 : do i = 1, nday
213 52706700 : alb_dir(iband,i) = cam_in%aldir(idxday(i))
214 59421924 : alb_dif(iband,i) = cam_in%aldif(idxday(i))
215 : end do
216 : else
217 : ! Band straddles the visible to near-infrared transition, so we take
218 : ! the albedo to be the average of the visible and near-infrared
219 : ! broadband albedos
220 6602436 : do i = 1, nday
221 5856300 : alb_dir(iband,i) = 0.5_r8 * (cam_in%aldir(idxday(i)) + cam_in%asdir(idxday(i)))
222 6602436 : alb_dif(iband,i) = 0.5_r8 * (cam_in%aldif(idxday(i)) + cam_in%asdif(idxday(i)))
223 : end do
224 : end if
225 : end do
226 :
227 : ! Strictly enforce albedo bounds
228 88590636 : where (alb_dir < 0)
229 : alb_dir = 0.0_r8
230 : end where
231 88590636 : where (alb_dir > 1)
232 : alb_dir = 1.0_r8
233 : end where
234 88590636 : where (alb_dif < 0)
235 : alb_dif = 0.0_r8
236 : end where
237 88590636 : where (alb_dif > 1)
238 : alb_dif = 1.0_r8
239 : end where
240 :
241 746136 : end subroutine rrtmgp_set_state
242 :
243 : !=========================================================================================
244 :
245 20891808 : pure logical function is_visible(wavenumber)
246 :
247 : ! Wavenumber is in the visible if it is above the visible threshold
248 : ! wavenumber, and in the infrared if it is below the threshold
249 : ! This function doesn't distinquish between visible and UV.
250 :
251 : ! wavenumber in inverse cm (cm^-1)
252 : real(r8), intent(in) :: wavenumber
253 :
254 : ! Set threshold between visible and infrared to 0.7 micron, or 14286 cm^-1
255 : real(r8), parameter :: visible_wavenumber_threshold = 14286._r8 ! cm^-1
256 :
257 20891808 : if (wavenumber > visible_wavenumber_threshold) then
258 : is_visible = .true.
259 : else
260 14176584 : is_visible = .false.
261 : end if
262 :
263 20891808 : end function is_visible
264 :
265 : !=========================================================================================
266 :
267 9083192 : function get_molar_mass_ratio(gas_name) result(massratio)
268 :
269 : ! return the molar mass ratio of dry air to gas based on gas_name
270 :
271 : character(len=*),intent(in) :: gas_name
272 : real(r8) :: massratio
273 :
274 : ! local variables
275 : real(r8), parameter :: amdw = 1.607793_r8 ! Molecular weight of dry air / water vapor
276 : real(r8), parameter :: amdc = 0.658114_r8 ! Molecular weight of dry air / carbon dioxide
277 : real(r8), parameter :: amdo = 0.603428_r8 ! Molecular weight of dry air / ozone
278 : real(r8), parameter :: amdm = 1.805423_r8 ! Molecular weight of dry air / methane
279 : real(r8), parameter :: amdn = 0.658090_r8 ! Molecular weight of dry air / nitrous oxide
280 : real(r8), parameter :: amdo2 = 0.905140_r8 ! Molecular weight of dry air / oxygen
281 : real(r8), parameter :: amdc1 = 0.210852_r8 ! Molecular weight of dry air / CFC11
282 : real(r8), parameter :: amdc2 = 0.239546_r8 ! Molecular weight of dry air / CFC12
283 :
284 : character(len=*), parameter :: sub='get_molar_mass_ratio'
285 : !----------------------------------------------------------------------------
286 :
287 18166384 : select case (trim(gas_name))
288 : case ('H2O')
289 1135399 : massratio = amdw
290 : case ('CO2')
291 1135399 : massratio = amdc
292 : case ('O3')
293 1135399 : massratio = amdo
294 : case ('CH4')
295 1135399 : massratio = amdm
296 : case ('N2O')
297 1135399 : massratio = amdn
298 : case ('O2')
299 1135399 : massratio = amdo2
300 : case ('CFC11')
301 1135399 : massratio = amdc1
302 : case ('CFC12')
303 1135399 : massratio = amdc2
304 : case default
305 18166384 : call endrun(sub//": Invalid gas: "//trim(gas_name))
306 : end select
307 :
308 9083192 : end function get_molar_mass_ratio
309 :
310 : !=========================================================================================
311 :
312 9083192 : subroutine rad_gas_get_vmr(icall, gas_name, state, pbuf, nlay, numactivecols, gas_concs, idxday)
313 :
314 : ! Set volume mixing ratio in gas_concs object.
315 : ! The gas_concs%set_vmr method copies data into internally allocated storage.
316 :
317 : integer, intent(in) :: icall ! index of climate/diagnostic radiation call
318 : character(len=*), intent(in) :: gas_name
319 : type(physics_state), target, intent(in) :: state
320 : type(physics_buffer_desc), pointer :: pbuf(:)
321 : integer, intent(in) :: nlay ! number of layers in radiation calculation
322 : integer, intent(in) :: numactivecols ! number of columns, ncol for LW, nday for SW
323 :
324 : type(ty_gas_concs), intent(inout) :: gas_concs ! the result is VRM inside gas_concs
325 :
326 : integer, optional, intent(in) :: idxday(:) ! indices of daylight columns in a chunk
327 :
328 : ! Local variables
329 18166384 : integer :: i, idx(numactivecols)
330 : integer :: istat
331 9083192 : real(r8), pointer :: gas_mmr(:,:)
332 9083192 : real(r8), allocatable :: gas_vmr(:,:)
333 9083192 : real(r8), allocatable :: mmr(:,:)
334 : real(r8) :: massratio
335 :
336 : ! For ozone profile above model
337 : real(r8) :: P_top, P_int, P_mid, alpha, beta, a, b, chi_mid, chi_0, chi_eff
338 :
339 : character(len=128) :: errmsg
340 : character(len=*), parameter :: sub = 'rad_gas_get_vmr'
341 : !----------------------------------------------------------------------------
342 :
343 : ! set the column indices; when idxday is provided (e.g. daylit columns) use them, otherwise just count.
344 149634392 : do i = 1, numactivecols
345 149634392 : if (present(idxday)) then
346 46850400 : idx(i) = idxday(i)
347 : else
348 93700800 : idx(i) = i
349 : end if
350 : end do
351 :
352 : ! gas_mmr points to a "chunk" in either the state or pbuf objects. That storage is
353 : ! dimensioned (pcols,pver).
354 9083192 : call rad_cnst_get_gas(icall, gas_name, state, pbuf, gas_mmr)
355 :
356 : ! Copy into storage for RRTMGP
357 36332768 : allocate(mmr(numactivecols, nlay), stat=istat)
358 9083192 : call alloc_err(istat, sub, 'mmr', numactivecols*nlay)
359 27249576 : allocate(gas_vmr(numactivecols, nlay), stat=istat)
360 9083192 : call alloc_err(istat, sub, 'gas_vmr', numactivecols*nlay)
361 :
362 149634392 : do i = 1, numactivecols
363 13220895992 : mmr(i,ktoprad:) = gas_mmr(idx(i),ktopcam:)
364 : end do
365 :
366 : ! If an extra layer is being used, copy mmr from the top layer of CAM to the extra layer.
367 9083192 : if (nlay == pverp) then
368 149634392 : mmr(:,1) = mmr(:,2)
369 : end if
370 :
371 : ! special case: H2O is specific humidity, not mixing ratio. Use r = q/(1-q):
372 9083192 : if (gas_name == 'H2O') then
373 1759339505 : mmr = mmr / (1._r8 - mmr)
374 : end if
375 :
376 : ! convert MMR to VMR, multipy by ratio of dry air molar mas to gas molar mass.
377 9083192 : massratio = get_molar_mass_ratio(gas_name)
378 14083799232 : gas_vmr = mmr * massratio
379 :
380 : ! special case: Setting O3 in the extra layer:
381 : !
382 : ! For the purpose of attenuating solar fluxes above the CAM model top, we assume that ozone
383 : ! mixing decreases linearly in each column from the value in the top layer of CAM to zero at
384 : ! the pressure level set by P_top. P_top has been set to 50 Pa (0.5 hPa) based on model tuning
385 : ! to produce temperatures at the top of CAM that are most consistent with WACCM at similar pressure levels.
386 :
387 9083192 : if ((gas_name == 'O3') .and. (nlay == pverp)) then
388 : P_top = 50.0_r8
389 18704299 : do i = 1, numactivecols
390 17568900 : P_int = state%pint(idx(i),1) ! pressure (Pa) at upper interface of CAM
391 17568900 : P_mid = state%pmid(idx(i),1) ! pressure (Pa) at midpoint of top layer of CAM
392 17568900 : alpha = log(P_int/P_top)
393 17568900 : beta = log(P_mid/P_int)/log(P_mid/P_top)
394 :
395 17568900 : a = ( (1._r8 + alpha) * exp(-alpha) - 1._r8 ) / alpha
396 17568900 : b = 1._r8 - exp(-alpha)
397 :
398 18704299 : if (alpha .gt. 0) then ! only apply where top level is below 80 km
399 0 : chi_mid = gas_vmr(i,1) ! molar mixing ratio of O3 at midpoint of top layer
400 0 : chi_0 = chi_mid / (1._r8 + beta)
401 0 : chi_eff = chi_0 * (a + b)
402 0 : gas_vmr(i,1) = chi_eff
403 : end if
404 : end do
405 : end if
406 :
407 9083192 : errmsg = gas_concs%set_vmr(gas_name, gas_vmr)
408 9083192 : if (len_trim(errmsg) > 0) then
409 0 : call endrun(sub//': ERROR, gas_concs%set_vmr: '//trim(errmsg))
410 : end if
411 :
412 9083192 : deallocate(gas_vmr)
413 9083192 : deallocate(mmr)
414 :
415 18166384 : end subroutine rad_gas_get_vmr
416 :
417 : !==================================================================================================
418 :
419 746136 : subroutine rrtmgp_set_gases_lw(icall, state, pbuf, nlay, gas_concs)
420 :
421 : ! Set gas vmr for the gases in the radconstants module's gaslist.
422 :
423 : ! The memory management for the gas_concs object is internal. The arrays passed to it
424 : ! are copied to the internally allocated memory. Each call to the set_vmr method checks
425 : ! whether the gas already has memory allocated, and if it does that memory is deallocated
426 : ! and new memory is allocated.
427 :
428 : ! arguments
429 : integer, intent(in) :: icall ! index of climate/diagnostic radiation call
430 : type(physics_state), target, intent(in) :: state
431 : type(physics_buffer_desc), pointer :: pbuf(:)
432 : integer, intent(in) :: nlay
433 : type(ty_gas_concs), intent(inout) :: gas_concs
434 :
435 : ! local variables
436 : integer :: i, ncol
437 : character(len=*), parameter :: sub = 'rrtmgp_set_gases_lw'
438 : !--------------------------------------------------------------------------------
439 :
440 746136 : ncol = state%ncol
441 6715224 : do i = 1, nradgas
442 6715224 : call rad_gas_get_vmr(icall, gaslist(i), state, pbuf, nlay, ncol, gas_concs)
443 : end do
444 746136 : end subroutine rrtmgp_set_gases_lw
445 :
446 : !==================================================================================================
447 :
448 389263 : subroutine rrtmgp_set_gases_sw( &
449 : icall, state, pbuf, nlay, nday, &
450 389263 : idxday, gas_concs)
451 :
452 : ! Return gas_concs with gas volume mixing ratio on DAYLIT columns.
453 : ! Set all gases in radconstants gaslist.
454 :
455 : ! arguments
456 : integer, intent(in) :: icall ! index of climate/diagnostic radiation call
457 : type(physics_state), target, intent(in) :: state
458 : type(physics_buffer_desc), pointer :: pbuf(:)
459 : integer, intent(in) :: nlay
460 : integer, intent(in) :: nday
461 : integer, intent(in) :: idxday(:)
462 : type(ty_gas_concs), intent(inout) :: gas_concs
463 :
464 : ! local variables
465 : integer :: i
466 : character(len=*), parameter :: sub = 'rrtmgp_set_gases_sw'
467 : !----------------------------------------------------------------------------
468 :
469 : ! use the optional argument idxday to specify which columns are sunlit
470 3503367 : do i = 1,nradgas
471 3503367 : call rad_gas_get_vmr(icall, gaslist(i), state, pbuf, nlay, nday, gas_concs, idxday=idxday)
472 : end do
473 :
474 389263 : end subroutine rrtmgp_set_gases_sw
475 :
476 : !==================================================================================================
477 :
478 746136 : subroutine rrtmgp_set_cloud_lw( &
479 : state, pbuf, ncol, nlay, nlaycam, &
480 : cld, cldfsnow, cldfgrau, cldfprime, graupel_in_rad, &
481 : kdist_lw, cloud_lw, cld_lw_abs_cloudsim, snow_lw_abs_cloudsim, grau_lw_abs_cloudsim)
482 :
483 : ! Compute combined cloud optical properties.
484 : ! Create MCICA stochastic arrays for cloud LW optical properties.
485 : ! Initialize optical properties object (cloud_lw) and load with MCICA columns.
486 :
487 : ! arguments
488 : type(physics_state), intent(in) :: state
489 : type(physics_buffer_desc), pointer :: pbuf(:)
490 : integer, intent(in) :: ncol ! number of columns in CAM chunk
491 : integer, intent(in) :: nlay ! number of layers in radiation calculation (may include "extra layer")
492 : integer, intent(in) :: nlaycam ! number of CAM layers in radiation calculation
493 : real(r8), pointer :: cld(:,:) ! cloud fraction (liq+ice)
494 : real(r8), pointer :: cldfsnow(:,:) ! cloud fraction of just "snow clouds"
495 : real(r8), pointer :: cldfgrau(:,:) ! cloud fraction of just "graupel clouds"
496 : real(r8), intent(in) :: cldfprime(pcols,pver) ! combined cloud fraction
497 :
498 : logical, intent(in) :: graupel_in_rad ! use graupel in radiation code
499 : class(ty_gas_optics_rrtmgp), intent(in) :: kdist_lw
500 : type(ty_optical_props_1scl), intent(out) :: cloud_lw
501 :
502 : ! Diagnostic outputs
503 : real(r8), intent(out) :: cld_lw_abs_cloudsim(pcols,pver) ! in-cloud liq+ice optical depth (for COSP)
504 : real(r8), intent(out) :: snow_lw_abs_cloudsim(pcols,pver)! in-cloud snow optical depth (for COSP)
505 : real(r8), intent(out) :: grau_lw_abs_cloudsim(pcols,pver)! in-cloud Graupel optical depth (for COSP)
506 :
507 : ! Local variables
508 :
509 : integer :: i, k
510 :
511 : ! cloud radiative parameters are "in cloud" not "in cell"
512 : real(r8) :: liq_lw_abs(nlwbands,pcols,pver) ! liquid absorption optics depth (LW)
513 : real(r8) :: ice_lw_abs(nlwbands,pcols,pver) ! ice absorption optics depth (LW)
514 : real(r8) :: cld_lw_abs(nlwbands,pcols,pver) ! cloud absorption optics depth (LW)
515 : real(r8) :: snow_lw_abs(nlwbands,pcols,pver) ! snow absorption optics depth (LW)
516 : real(r8) :: grau_lw_abs(nlwbands,pcols,pver) ! graupel absorption optics depth (LW)
517 : real(r8) :: c_cld_lw_abs(nlwbands,pcols,pver) ! combined cloud absorption optics depth (LW)
518 :
519 : ! Arrays for converting from CAM chunks to RRTMGP inputs.
520 1492272 : real(r8) :: cldf(ncol,nlaycam)
521 1492272 : real(r8) :: tauc(nlwbands,ncol,nlaycam)
522 746136 : real(r8) :: taucmcl(nlwgpts,ncol,nlaycam)
523 :
524 : character(len=128) :: errmsg
525 : character(len=*), parameter :: sub = 'rrtmgp_set_cloud_lw'
526 : !--------------------------------------------------------------------------------
527 :
528 : ! Combine the cloud optical properties. These calculations are done on CAM "chunks".
529 :
530 : ! gammadist liquid optics
531 746136 : call liquid_cloud_get_rad_props_lw(state, pbuf, liq_lw_abs)
532 : ! Mitchell ice optics
533 746136 : call ice_cloud_get_rad_props_lw(state, pbuf, ice_lw_abs)
534 :
535 18587757384 : cld_lw_abs(:,:ncol,:) = liq_lw_abs(:,:ncol,:) + ice_lw_abs(:,:ncol,:)
536 :
537 746136 : if (associated(cldfsnow)) then
538 : ! add in snow
539 746136 : call snow_cloud_get_rad_props_lw(state, pbuf, snow_lw_abs)
540 12458736 : do i = 1, ncol
541 1101730536 : do k = 1, pver
542 1100984400 : if (cldfprime(i,k) > 0._r8) then
543 0 : c_cld_lw_abs(:,i,k) = ( cldfsnow(i,k)*snow_lw_abs(:,i,k) &
544 3972169782 : + cld(i,k)*cld_lw_abs(:,i,k) )/cldfprime(i,k)
545 : else
546 14545450818 : c_cld_lw_abs(:,i,k) = 0._r8
547 : end if
548 : end do
549 : end do
550 : else
551 0 : c_cld_lw_abs(:,:ncol,:) = cld_lw_abs(:,:ncol,:)
552 : end if
553 :
554 : ! add in graupel
555 746136 : if (associated(cldfgrau) .and. graupel_in_rad) then
556 0 : call grau_cloud_get_rad_props_lw(state, pbuf, grau_lw_abs)
557 0 : do i = 1, ncol
558 0 : do k = 1, pver
559 0 : if (cldfprime(i,k) > 0._r8) then
560 0 : c_cld_lw_abs(:,i,k) = ( cldfgrau(i,k)*grau_lw_abs(:,i,k) &
561 0 : + cld(i,k)*c_cld_lw_abs(:,i,k) )/cldfprime(i,k)
562 : else
563 0 : c_cld_lw_abs(:,i,k) = 0._r8
564 : end if
565 : end do
566 : end do
567 : end if
568 :
569 : ! Cloud optics for COSP
570 1180387152 : cld_lw_abs_cloudsim = cld_lw_abs(idx_lw_cloudsim,:,:)
571 1180387152 : snow_lw_abs_cloudsim = snow_lw_abs(idx_lw_cloudsim,:,:)
572 1180387152 : grau_lw_abs_cloudsim = grau_lw_abs(idx_lw_cloudsim,:,:)
573 :
574 : ! Extract just the layers of CAM where RRTMGP does calculations.
575 :
576 : ! Subset "chunk" data so just the number of CAM layers in the
577 : ! radiation calculation are used by MCICA to produce subcolumns.
578 1159408584 : cldf = cldfprime(:ncol, ktopcam:)
579 18587757384 : tauc = c_cld_lw_abs(:, :ncol, ktopcam:)
580 :
581 : ! Enforce tauc >= 0.
582 18587757384 : tauc = merge(tauc, 0.0_r8, tauc > 0.0_r8)
583 :
584 : ! MCICA uses spectral data (on bands) to construct subcolumns (one per g-point)
585 : call mcica_subcol_lw( &
586 : kdist_lw, nlwbands, nlwgpts, ncol, nlaycam, &
587 746136 : nlwgpts, state%pmid, cldf, tauc, taucmcl )
588 :
589 746136 : errmsg =cloud_lw%alloc_1scl(ncol, nlay, kdist_lw)
590 746136 : if (len_trim(errmsg) > 0) then
591 0 : call endrun(sub//': ERROR: cloud_lw%alloc_1scalar: '//trim(errmsg))
592 : end if
593 :
594 : ! If there is an extra layer in the radiation then this initialization
595 : ! will provide zero optical depths there.
596 >14999*10^7 : cloud_lw%tau = 0.0_r8
597 :
598 : ! Set the properties on g-points.
599 96251544 : do i = 1, nlwgpts
600 >14840*10^7 : cloud_lw%tau(:,ktoprad:,i) = taucmcl(i,:,:)
601 : end do
602 :
603 : ! validate checks that: tau > 0
604 746136 : errmsg = cloud_lw%validate()
605 746136 : if (len_trim(errmsg) > 0) then
606 0 : call endrun(sub//': ERROR: cloud_lw%validate: '//trim(errmsg))
607 : end if
608 :
609 746136 : end subroutine rrtmgp_set_cloud_lw
610 :
611 : !==================================================================================================
612 :
613 746136 : subroutine rrtmgp_set_cloud_sw( &
614 : state, pbuf, nlay, nday, idxday, &
615 746136 : nnite, idxnite, pmid, cld, cldfsnow, &
616 : cldfgrau, cldfprime, graupel_in_rad, kdist_sw, cloud_sw, &
617 : tot_cld_vistau, tot_icld_vistau, liq_icld_vistau, ice_icld_vistau, snow_icld_vistau, &
618 : grau_icld_vistau, cld_tau_cloudsim, snow_tau_cloudsim, grau_tau_cloudsim)
619 :
620 : ! Compute combined cloud optical properties.
621 : ! Create MCICA stochastic arrays for cloud SW optical properties.
622 : ! Initialize optical properties object (cloud_sw) and load with MCICA columns.
623 :
624 : ! arguments
625 : type(physics_state), target, intent(in) :: state
626 : type(physics_buffer_desc), pointer :: pbuf(:)
627 : integer, intent(in) :: nlay ! number of layers in radiation calculation (may include "extra layer")
628 : integer, intent(in) :: nday ! number of daylight columns
629 : integer, intent(in) :: idxday(pcols) ! indices of daylight columns in the chunk
630 : integer, intent(in) :: nnite ! number of night columns
631 : integer, intent(in) :: idxnite(pcols) ! indices of night columns in the chunk
632 :
633 : real(r8), intent(in) :: pmid(nday,nlay)! pressure at layer midpoints (Pa) used to seed RNG.
634 :
635 : real(r8), pointer :: cld(:,:) ! cloud fraction (liq+ice)
636 : real(r8), pointer :: cldfsnow(:,:) ! cloud fraction of just "snow clouds"
637 : real(r8), pointer :: cldfgrau(:,:) ! cloud fraction of just "graupel clouds"
638 : real(r8), intent(in) :: cldfprime(pcols,pver) ! combined cloud fraction
639 :
640 : logical, intent(in) :: graupel_in_rad ! graupel in radiation code
641 : class(ty_gas_optics_rrtmgp), intent(in) :: kdist_sw ! shortwave gas optics object
642 : type(ty_optical_props_2str), intent(out) :: cloud_sw ! SW cloud optical properties object
643 :
644 : ! Diagnostic outputs
645 : real(r8), intent(out) :: tot_cld_vistau(pcols,pver) ! gbx total cloud optical depth
646 : real(r8), intent(out) :: tot_icld_vistau(pcols,pver) ! in-cld total cloud optical depth
647 : real(r8), intent(out) :: liq_icld_vistau(pcols,pver) ! in-cld liq cloud optical depth
648 : real(r8), intent(out) :: ice_icld_vistau(pcols,pver) ! in-cld ice cloud optical depth
649 : real(r8), intent(out) :: snow_icld_vistau(pcols,pver) ! snow in-cloud visible sw optical depth
650 : real(r8), intent(out) :: grau_icld_vistau(pcols,pver) ! Graupel in-cloud visible sw optical depth
651 : real(r8), intent(out) :: cld_tau_cloudsim(pcols,pver) ! in-cloud liq+ice optical depth (for COSP)
652 : real(r8), intent(out) :: snow_tau_cloudsim(pcols,pver)! in-cloud snow optical depth (for COSP)
653 : real(r8), intent(out) :: grau_tau_cloudsim(pcols,pver)! in-cloud Graupel optical depth (for COSP)
654 :
655 : ! Local variables
656 :
657 : integer :: i, k, ncol
658 : integer :: igpt, nver
659 : integer :: istat
660 : integer, parameter :: changeseed = 1
661 :
662 : ! cloud radiative parameters are "in cloud" not "in cell"
663 : real(r8) :: ice_tau (nswbands,pcols,pver) ! ice extinction optical depth
664 : real(r8) :: ice_tau_w (nswbands,pcols,pver) ! ice single scattering albedo * tau
665 : real(r8) :: ice_tau_w_g(nswbands,pcols,pver) ! ice asymmetry parameter * tau * w
666 : real(r8) :: liq_tau (nswbands,pcols,pver) ! liquid extinction optical depth
667 : real(r8) :: liq_tau_w (nswbands,pcols,pver) ! liquid single scattering albedo * tau
668 : real(r8) :: liq_tau_w_g(nswbands,pcols,pver) ! liquid asymmetry parameter * tau * w
669 : real(r8) :: cld_tau (nswbands,pcols,pver) ! cloud extinction optical depth
670 : real(r8) :: cld_tau_w (nswbands,pcols,pver) ! cloud single scattering albedo * tau
671 : real(r8) :: cld_tau_w_g(nswbands,pcols,pver) ! cloud asymmetry parameter * w * tau
672 : real(r8) :: snow_tau (nswbands,pcols,pver) ! snow extinction optical depth
673 : real(r8) :: snow_tau_w (nswbands,pcols,pver) ! snow single scattering albedo * tau
674 : real(r8) :: snow_tau_w_g(nswbands,pcols,pver) ! snow asymmetry parameter * tau * w
675 : real(r8) :: grau_tau (nswbands,pcols,pver) ! graupel extinction optical depth
676 : real(r8) :: grau_tau_w (nswbands,pcols,pver) ! graupel single scattering albedo * tau
677 : real(r8) :: grau_tau_w_g(nswbands,pcols,pver) ! graupel asymmetry parameter * tau * w
678 : real(r8) :: c_cld_tau (nswbands,pcols,pver) ! combined cloud extinction optical depth
679 : real(r8) :: c_cld_tau_w (nswbands,pcols,pver) ! combined cloud single scattering albedo * tau
680 : real(r8) :: c_cld_tau_w_g(nswbands,pcols,pver) ! combined cloud asymmetry parameter * w * tau
681 :
682 : ! RRTMGP does not use this property in its 2-stream calculations.
683 : real(r8) :: sw_tau_w_f(nswbands,pcols,pver) ! Forward scattered fraction * tau * w.
684 :
685 : ! Arrays for converting from CAM chunks to RRTMGP inputs.
686 746136 : real(r8), allocatable :: cldf(:,:)
687 746136 : real(r8), allocatable :: tauc(:,:,:)
688 746136 : real(r8), allocatable :: ssac(:,:,:)
689 746136 : real(r8), allocatable :: asmc(:,:,:)
690 746136 : real(r8), allocatable :: taucmcl(:,:,:)
691 746136 : real(r8), allocatable :: ssacmcl(:,:,:)
692 746136 : real(r8), allocatable :: asmcmcl(:,:,:)
693 746136 : real(r8), allocatable :: day_cld_tau(:,:,:)
694 746136 : real(r8), allocatable :: day_cld_tau_w(:,:,:)
695 746136 : real(r8), allocatable :: day_cld_tau_w_g(:,:,:)
696 :
697 : character(len=128) :: errmsg
698 : character(len=*), parameter :: sub = 'rrtmgp_set_cloud_sw'
699 : !--------------------------------------------------------------------------------
700 :
701 746136 : ncol = state%ncol
702 :
703 : ! Combine the cloud optical properties. These calculations are done on CAM "chunks".
704 :
705 : ! gammadist liquid optics
706 746136 : call get_liquid_optics_sw(state, pbuf, liq_tau, liq_tau_w, liq_tau_w_g, sw_tau_w_f)
707 : ! Mitchell ice optics
708 746136 : call get_ice_optics_sw(state, pbuf, ice_tau, ice_tau_w, ice_tau_w_g, sw_tau_w_f)
709 :
710 16409213784 : cld_tau(:,:ncol,:) = liq_tau(:,:ncol,:) + ice_tau(:,:ncol,:)
711 16409213784 : cld_tau_w(:,:ncol,:) = liq_tau_w(:,:ncol,:) + ice_tau_w(:,:ncol,:)
712 16409213784 : cld_tau_w_g(:,:ncol,:) = liq_tau_w_g(:,:ncol,:) + ice_tau_w_g(:,:ncol,:)
713 :
714 : ! add in snow
715 746136 : if (associated(cldfsnow)) then
716 746136 : call get_snow_optics_sw(state, pbuf, snow_tau, snow_tau_w, snow_tau_w_g, sw_tau_w_f)
717 12458736 : do i = 1, ncol
718 1101730536 : do k = 1, pver
719 1100984400 : if (cldfprime(i,k) > 0._r8) then
720 0 : c_cld_tau(:,i,k) = ( cldfsnow(i,k)*snow_tau(:,i,k) &
721 3504855690 : + cld(i,k)*cld_tau(:,i,k) )/cldfprime(i,k)
722 0 : c_cld_tau_w(:,i,k) = ( cldfsnow(i,k)*snow_tau_w(:,i,k) &
723 3504855690 : + cld(i,k)*cld_tau_w(:,i,k) )/cldfprime(i,k)
724 0 : c_cld_tau_w_g(:,i,k) = ( cldfsnow(i,k)*snow_tau_w_g(:,i,k) &
725 3504855690 : + cld(i,k)*cld_tau_w_g(:,i,k) )/cldfprime(i,k)
726 : else
727 12834221310 : c_cld_tau(:,i,k) = 0._r8
728 12834221310 : c_cld_tau_w(:,i,k) = 0._r8
729 12834221310 : c_cld_tau_w_g(:,i,k) = 0._r8
730 : end if
731 : end do
732 : end do
733 : else
734 0 : c_cld_tau(:,:ncol,:) = cld_tau(:,:ncol,:)
735 0 : c_cld_tau_w(:,:ncol,:) = cld_tau_w(:,:ncol,:)
736 0 : c_cld_tau_w_g(:,:ncol,:) = cld_tau_w_g(:,:ncol,:)
737 : end if
738 :
739 : ! add in graupel
740 746136 : if (associated(cldfgrau) .and. graupel_in_rad) then
741 0 : call get_grau_optics_sw(state, pbuf, grau_tau, grau_tau_w, grau_tau_w_g, sw_tau_w_f)
742 0 : do i = 1, ncol
743 0 : do k = 1, pver
744 0 : if (cldfprime(i,k) > 0._r8) then
745 0 : c_cld_tau(:,i,k) = ( cldfgrau(i,k)*grau_tau(:,i,k) &
746 0 : + cld(i,k)*c_cld_tau(:,i,k) )/cldfprime(i,k)
747 0 : c_cld_tau_w(:,i,k) = ( cldfgrau(i,k)*grau_tau_w(:,i,k) &
748 0 : + cld(i,k)*c_cld_tau_w(:,i,k) )/cldfprime(i,k)
749 0 : c_cld_tau_w_g(:,i,k) = ( cldfgrau(i,k)*grau_tau_w_g(:,i,k) &
750 0 : + cld(i,k)*c_cld_tau_w_g(:,i,k) )/cldfprime(i,k)
751 : else
752 0 : c_cld_tau(:,i,k) = 0._r8
753 0 : c_cld_tau_w(:,i,k) = 0._r8
754 0 : c_cld_tau_w_g(:,i,k) = 0._r8
755 : end if
756 : end do
757 : end do
758 : end if
759 :
760 : ! cloud optical properties need to be re-ordered from the RRTMG spectral bands
761 : ! (assumed in the optics datasets) to RRTMGP's
762 32817681432 : ice_tau(:,:ncol,:) = ice_tau(rrtmg_to_rrtmgp_swbands,:ncol,:)
763 32817681432 : liq_tau(:,:ncol,:) = liq_tau(rrtmg_to_rrtmgp_swbands,:ncol,:)
764 32817681432 : c_cld_tau(:,:ncol,:) = c_cld_tau(rrtmg_to_rrtmgp_swbands,:ncol,:)
765 32817681432 : c_cld_tau_w(:,:ncol,:) = c_cld_tau_w(rrtmg_to_rrtmgp_swbands,:ncol,:)
766 32817681432 : c_cld_tau_w_g(:,:ncol,:) = c_cld_tau_w_g(rrtmg_to_rrtmgp_swbands,:ncol,:)
767 746136 : if (associated(cldfsnow)) then
768 32817681432 : snow_tau(:,:ncol,:) = snow_tau(rrtmg_to_rrtmgp_swbands,:ncol,:)
769 : end if
770 746136 : if (associated(cldfgrau) .and. graupel_in_rad) then
771 0 : grau_tau(:,:ncol,:) = grau_tau(rrtmg_to_rrtmgp_swbands,:ncol,:)
772 : end if
773 :
774 : ! Set arrays for diagnostic output.
775 : ! cloud optical depth fields for the visible band
776 1159408584 : tot_icld_vistau(:ncol,:) = c_cld_tau(idx_sw_diag,:ncol,:)
777 1159408584 : liq_icld_vistau(:ncol,:) = liq_tau(idx_sw_diag,:ncol,:)
778 1159408584 : ice_icld_vistau(:ncol,:) = ice_tau(idx_sw_diag,:ncol,:)
779 746136 : if (associated(cldfsnow)) then
780 1159408584 : snow_icld_vistau(:ncol,:) = snow_tau(idx_sw_diag,:ncol,:)
781 : endif
782 746136 : if (associated(cldfgrau) .and. graupel_in_rad) then
783 0 : grau_icld_vistau(:ncol,:) = grau_tau(idx_sw_diag,:ncol,:)
784 : endif
785 :
786 : ! multiply by total cloud fraction to get gridbox value
787 1159408584 : tot_cld_vistau(:ncol,:) = c_cld_tau(idx_sw_diag,:ncol,:)*cldfprime(:ncol,:)
788 :
789 : ! overwrite night columns with fillvalue
790 6602436 : do i = 1, Nnite
791 550492200 : tot_cld_vistau(IdxNite(i),:) = fillvalue
792 550492200 : tot_icld_vistau(IdxNite(i),:) = fillvalue
793 550492200 : liq_icld_vistau(IdxNite(i),:) = fillvalue
794 550492200 : ice_icld_vistau(IdxNite(i),:) = fillvalue
795 5856300 : if (associated(cldfsnow)) then
796 550492200 : snow_icld_vistau(IdxNite(i),:) = fillvalue
797 : end if
798 6602436 : if (associated(cldfgrau) .and. graupel_in_rad) then
799 0 : grau_icld_vistau(IdxNite(i),:) = fillvalue
800 : end if
801 : end do
802 :
803 : ! Cloud optics for COSP
804 1180387152 : cld_tau_cloudsim = cld_tau(idx_sw_cloudsim,:,:)
805 1180387152 : snow_tau_cloudsim = snow_tau(idx_sw_cloudsim,:,:)
806 1180387152 : grau_tau_cloudsim = grau_tau(idx_sw_cloudsim,:,:)
807 :
808 : ! if no daylight columns the cloud_sw object isn't initialized
809 746136 : if (nday > 0) then
810 :
811 : ! number of CAM's layers in radiation calculation. Does not include the "extra layer".
812 389263 : nver = pver - ktopcam + 1
813 :
814 : allocate( &
815 : cldf(nday,nver), &
816 : day_cld_tau(nswbands,nday,nver), &
817 : day_cld_tau_w(nswbands,nday,nver), &
818 : day_cld_tau_w_g(nswbands,nday,nver), &
819 : tauc(nswbands,nday,nver), taucmcl(nswgpts,nday,nver), &
820 : ssac(nswbands,nday,nver), ssacmcl(nswgpts,nday,nver), &
821 10510101 : asmc(nswbands,nday,nver), asmcmcl(nswgpts,nday,nver), stat=istat)
822 389263 : call alloc_err(istat, sub, 'cldf,..,asmcmcl', 9*nswgpts*nday*nver)
823 :
824 : ! Subset "chunk" data so just the daylight columns, and the number of CAM layers in the
825 : ! radiation calculation are used by MCICA to produce subcolumns.
826 581615885 : cldf = cldfprime( idxday(1:nday), ktopcam:)
827 8206518485 : day_cld_tau = c_cld_tau( :, idxday(1:nday), ktopcam:)
828 8206518485 : day_cld_tau_w = c_cld_tau_w( :, idxday(1:nday), ktopcam:)
829 8206518485 : day_cld_tau_w_g = c_cld_tau_w_g(:, idxday(1:nday), ktopcam:)
830 :
831 : ! Compute the optical properties needed for the 2-stream calculations. These calculations
832 : ! are the same as the RRTMG version.
833 :
834 : ! set cloud optical depth, clip @ zero
835 8206518485 : tauc = merge(day_cld_tau, 0.0_r8, day_cld_tau > 0.0_r8)
836 : ! set value of asymmetry
837 8206518485 : asmc = merge(day_cld_tau_w_g / max(day_cld_tau_w, tiny), 0.0_r8, day_cld_tau_w > 0.0_r8)
838 : ! set value of single scattering albedo
839 8206518485 : ssac = merge(max(day_cld_tau_w, tiny) / max(tauc, tiny), 1.0_r8 , tauc > 0.0_r8)
840 : ! set asymmetry to zero when tauc = 0
841 8206518485 : asmc = merge(asmc, 0.0_r8, tauc > 0.0_r8)
842 :
843 : ! MCICA uses spectral data (on bands) to construct subcolumns (one per g-point)
844 : call mcica_subcol_sw( &
845 : kdist_sw, nswbands, nswgpts, nday, nlay, &
846 : nver, changeseed, pmid, cldf, tauc, &
847 389263 : ssac, asmc, taucmcl, ssacmcl, asmcmcl)
848 :
849 : ! Initialize object for SW cloud optical properties.
850 389263 : errmsg = cloud_sw%alloc_2str(nday, nlay, kdist_sw)
851 389263 : if (len_trim(errmsg) > 0) then
852 0 : call endrun(trim(sub)//': ERROR: cloud_sw%alloc_2str: '//trim(errmsg))
853 : end if
854 :
855 : ! If there is an extra layer in the radiation then this initialization
856 : ! will provide the optical properties there.
857 65797273983 : cloud_sw%tau = 0.0_r8
858 65797273983 : cloud_sw%ssa = 1.0_r8
859 65797273983 : cloud_sw%g = 0.0_r8
860 :
861 : ! Set the properties on g-points.
862 43986719 : do igpt = 1,nswgpts
863 65097381664 : cloud_sw%g (:,ktoprad:,igpt) = asmcmcl(igpt,:,:)
864 65097381664 : cloud_sw%ssa(:,ktoprad:,igpt) = ssacmcl(igpt,:,:)
865 65097770927 : cloud_sw%tau(:,ktoprad:,igpt) = taucmcl(igpt,:,:)
866 : end do
867 :
868 : ! validate checks that: tau > 0, ssa is in range [0,1], and g is in range [-1,1].
869 389263 : errmsg = cloud_sw%validate()
870 389263 : if (len_trim(errmsg) > 0) then
871 0 : call endrun(sub//': ERROR: cloud_sw%validate: '//trim(errmsg))
872 : end if
873 :
874 : ! delta scaling adjusts for forward scattering
875 389263 : errmsg = cloud_sw%delta_scale()
876 389263 : if (len_trim(errmsg) > 0) then
877 0 : call endrun(sub//': ERROR: cloud_sw%delta_scale: '//trim(errmsg))
878 : end if
879 :
880 : ! All information is in cloud_sw, now deallocate local vars.
881 0 : deallocate( &
882 0 : cldf, tauc, ssac, asmc, &
883 0 : taucmcl, ssacmcl, asmcmcl,&
884 389263 : day_cld_tau, day_cld_tau_w, day_cld_tau_w_g )
885 :
886 : end if
887 :
888 746136 : end subroutine rrtmgp_set_cloud_sw
889 :
890 : !==================================================================================================
891 :
892 746136 : subroutine rrtmgp_set_aer_lw(icall, state, pbuf, aer_lw)
893 :
894 : ! Load LW aerosol optical properties into the RRTMGP object.
895 :
896 : ! Arguments
897 : integer, intent(in) :: icall
898 : type(physics_state), target, intent(in) :: state
899 : type(physics_buffer_desc), pointer :: pbuf(:)
900 :
901 : type(ty_optical_props_1scl), intent(inout) :: aer_lw
902 :
903 : ! Local variables
904 : integer :: ncol
905 :
906 : ! Aerosol LW absorption optical depth
907 : real(r8) :: aer_lw_abs (pcols,pver,nlwbands)
908 :
909 : character(len=*), parameter :: sub = 'rrtmgp_set_aer_lw'
910 : character(len=128) :: errmsg
911 : !--------------------------------------------------------------------------------
912 :
913 746136 : ncol = state%ncol
914 :
915 : ! Get aerosol longwave optical properties.
916 746136 : call aer_rad_props_lw(icall, state, pbuf, aer_lw_abs)
917 :
918 : ! If there is an extra layer in the radiation then this initialization
919 : ! will provide zero optical depths there.
920 18750623256 : aer_lw%tau = 0.0_r8
921 :
922 18551283480 : aer_lw%tau(:ncol, ktoprad:, :) = aer_lw_abs(:ncol, ktopcam:, :)
923 :
924 746136 : errmsg = aer_lw%validate()
925 746136 : if (len_trim(errmsg) > 0) then
926 0 : call endrun(sub//': ERROR: aer_lw%validate: '//trim(errmsg))
927 : end if
928 746136 : end subroutine rrtmgp_set_aer_lw
929 :
930 : !==================================================================================================
931 :
932 746136 : subroutine rrtmgp_set_aer_sw( &
933 746136 : icall, state, pbuf, nday, idxday, nnite, idxnite, aer_sw)
934 :
935 : ! Load SW aerosol optical properties into the RRTMGP object.
936 :
937 : ! Arguments
938 : integer, intent(in) :: icall
939 : type(physics_state), target, intent(in) :: state
940 : type(physics_buffer_desc), pointer :: pbuf(:)
941 : integer, intent(in) :: nday
942 : integer, intent(in) :: idxday(:)
943 : integer, intent(in) :: nnite ! number of night columns
944 : integer, intent(in) :: idxnite(pcols) ! indices of night columns in the chunk
945 :
946 : type(ty_optical_props_2str), intent(inout) :: aer_sw
947 :
948 : ! local variables
949 : integer :: i
950 :
951 : ! The optical arrays dimensioned in the vertical as 0:pver.
952 : ! The index 0 is for the extra layer used in the radiation
953 : ! calculation. The index ktopcam assumes the CAM vertical indices are
954 : ! in the range 1:pver, so using this index correctly ignores vertical
955 : ! index 0. If an "extra" layer is used in the calculations, it is
956 : ! provided and set in the RRTMGP aerosol object aer_sw.
957 : real(r8) :: aer_tau (pcols,0:pver,nswbands) ! extinction optical depth
958 : real(r8) :: aer_tau_w (pcols,0:pver,nswbands) ! single scattering albedo * tau
959 : real(r8) :: aer_tau_w_g(pcols,0:pver,nswbands) ! asymmetry parameter * w * tau
960 : real(r8) :: aer_tau_w_f(pcols,0:pver,nswbands) ! forward scattered fraction * w * tau
961 : ! aer_tau_w_f is not used by RRTMGP.
962 : character(len=*), parameter :: sub = 'rrtmgp_set_aer_sw'
963 : !--------------------------------------------------------------------------------
964 :
965 : ! Get aerosol shortwave optical properties.
966 : ! Make outfld calls for aerosol optical property diagnostics.
967 : call aer_rad_props_sw( &
968 : icall, state, pbuf, nnite, idxnite, &
969 746136 : aer_tau, aer_tau_w, aer_tau_w_g, aer_tau_w_f)
970 :
971 : ! The aer_sw object is only initialized if nday > 0.
972 746136 : if (nday > 0) then
973 :
974 : ! aerosol optical properties need to be re-ordered from the RRTMG spectral bands
975 : ! (as assumed in the optics datasets) to the RRTMGP band order.
976 17428861562 : aer_tau(:,:,:) = aer_tau( :,:,rrtmg_to_rrtmgp_swbands)
977 17428861562 : aer_tau_w(:,:,:) = aer_tau_w( :,:,rrtmg_to_rrtmgp_swbands)
978 17428861562 : aer_tau_w_g(:,:,:) = aer_tau_w_g(:,:,rrtmg_to_rrtmgp_swbands)
979 :
980 : ! If there is an extra layer in the radiation then this initialization
981 : ! will provide default values.
982 8224999853 : aer_sw%tau = 0.0_r8
983 8224999853 : aer_sw%ssa = 1.0_r8
984 8224999853 : aer_sw%g = 0.0_r8
985 :
986 : ! CAM fields are products tau, tau*ssa, tau*ssa*asy
987 : ! Fields expected by RRTMGP are computed by
988 : ! aer_sw%tau = aer_tau
989 : ! aer_sw%ssa = aer_tau_w / aer_tau
990 : ! aer_sw%g = aer_tau_w_g / aer_taw_w
991 : ! aer_sw arrays have dimensions of (nday,nlay,nswbands)
992 :
993 6245563 : do i = 1, nday
994 : ! set aerosol optical depth, clip to zero
995 7712747100 : aer_sw%tau(i,ktoprad:,:) = max(aer_tau(idxday(i),ktopcam:,:), 0._r8)
996 : ! set value of single scattering albedo
997 5856300 : aer_sw%ssa(i,ktoprad:,:) = merge(aer_tau_w(idxday(i),ktopcam:,:)/aer_tau(idxday(i),ktopcam:,:), &
998 7718603400 : 1._r8, aer_tau(idxday(i),ktopcam:,:) > 0._r8)
999 : ! set value of asymmetry
1000 5856300 : aer_sw%g(i,ktoprad:,:) = merge(aer_tau_w_g(idxday(i),ktopcam:,:)/aer_tau_w(idxday(i),ktopcam:,:), &
1001 7718992663 : 0._r8, aer_tau_w(idxday(i),ktopcam:,:) > tiny)
1002 : end do
1003 :
1004 : ! impose limits on the components
1005 8224999853 : aer_sw%ssa = min(max(aer_sw%ssa, 0._r8), 1._r8)
1006 8224999853 : aer_sw%g = min(max(aer_sw%g, -1._r8), 1._r8)
1007 :
1008 : end if
1009 :
1010 746136 : end subroutine rrtmgp_set_aer_sw
1011 :
1012 : !==================================================================================================
1013 :
1014 : end module rrtmgp_inputs
|
