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 30960 : subroutine rrtmgp_set_state( &
101 : state, cam_in, ncol, nlay, nday, &
102 92880 : idxday, coszrs, kdist_sw, t_sfc, emis_sfc, &
103 30960 : t_rad, pmid_rad, pint_rad, t_day, pmid_day, &
104 30960 : 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 516960 : 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 8292960 : emis_sfc(:,:) = 1._r8
145 :
146 : ! Level ordering is the same for both CAM and RRTMGP (top to bottom)
147 48108240 : t_rad(:,ktoprad:) = state%t(:ncol,ktopcam:)
148 48108240 : pmid_rad(:,ktoprad:) = state%pmid(:ncol,ktopcam:)
149 48625200 : pint_rad(:,ktoprad:) = state%pint(:ncol,ktopcam:)
150 :
151 : ! Add extra layer values if needed.
152 30960 : if (nlay == pverp) then
153 516960 : 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 516960 : 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 516960 : 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 516960 : 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 516960 : 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 30960 : tref_min = kdist_sw%get_temp_min()
178 30960 : tref_max = kdist_sw%get_temp_max()
179 48625200 : t_rad = merge(t_rad, tref_min, t_rad > tref_min)
180 48625200 : t_rad = merge(t_rad, tref_max, t_rad < tref_max)
181 :
182 : ! Construct arrays containing only daylight columns
183 273960 : do i = 1, nday
184 23085000 : t_day(i,:) = t_rad(idxday(i),:)
185 23085000 : pmid_day(i,:) = pmid_rad(idxday(i),:)
186 23328000 : pint_day(i,:) = pint_rad(idxday(i),:)
187 273960 : coszrs_day(i) = coszrs(idxday(i))
188 : end do
189 :
190 : ! Assign albedos to the daylight columns (from E3SM implementation)
191 : ! Albedos are imported from the surface models as broadband (visible, and near-IR),
192 : ! and we need to map these to appropriate narrower bands used in RRTMGP. Bands
193 : ! are categorized broadly as "visible/UV" or "infrared" based on wavenumber.
194 : ! Loop over bands, and determine for each band whether it is broadly in the
195 : ! visible or infrared part of the spectrum based on a dividing line of
196 : ! 0.7 micron, or 14286 cm^-1
197 464400 : do iband = 1,nswbands
198 433440 : if (is_visible(sw_low_bounds(iband)) .and. &
199 30960 : is_visible(sw_high_bounds(iband))) then
200 :
201 : ! Entire band is in the visible
202 1095840 : do i = 1, nday
203 972000 : alb_dir(iband,i) = cam_in%asdir(idxday(i))
204 1095840 : alb_dif(iband,i) = cam_in%asdif(idxday(i))
205 : end do
206 :
207 309600 : else if (.not.is_visible(sw_low_bounds(iband)) .and. &
208 : .not.is_visible(sw_high_bounds(iband))) then
209 : ! Entire band is in the longwave (near-infrared)
210 2465640 : do i = 1, nday
211 2187000 : alb_dir(iband,i) = cam_in%aldir(idxday(i))
212 2465640 : alb_dif(iband,i) = cam_in%aldif(idxday(i))
213 : end do
214 : else
215 : ! Band straddles the visible to near-infrared transition, so we take
216 : ! the albedo to be the average of the visible and near-infrared
217 : ! broadband albedos
218 273960 : do i = 1, nday
219 243000 : alb_dir(iband,i) = 0.5_r8 * (cam_in%aldir(idxday(i)) + cam_in%asdir(idxday(i)))
220 273960 : alb_dif(iband,i) = 0.5_r8 * (cam_in%aldif(idxday(i)) + cam_in%asdif(idxday(i)))
221 : end do
222 : end if
223 : end do
224 :
225 : ! Strictly enforce albedo bounds
226 3675960 : where (alb_dir < 0)
227 : alb_dir = 0.0_r8
228 : end where
229 3675960 : where (alb_dir > 1)
230 : alb_dir = 1.0_r8
231 : end where
232 3675960 : where (alb_dif < 0)
233 : alb_dif = 0.0_r8
234 : end where
235 3675960 : where (alb_dif > 1)
236 : alb_dif = 1.0_r8
237 : end where
238 :
239 30960 : end subroutine rrtmgp_set_state
240 :
241 : !=========================================================================================
242 :
243 866880 : pure logical function is_visible(wavenumber)
244 :
245 : ! Wavenumber is in the visible if it is above the visible threshold
246 : ! wavenumber, and in the infrared if it is below the threshold
247 : ! This function doesn't distinquish between visible and UV.
248 :
249 : ! wavenumber in inverse cm (cm^-1)
250 : real(r8), intent(in) :: wavenumber
251 :
252 : ! Set threshold between visible and infrared to 0.7 micron, or 14286 cm^-1
253 : real(r8), parameter :: visible_wavenumber_threshold = 14286._r8 ! cm^-1
254 :
255 866880 : if (wavenumber > visible_wavenumber_threshold) then
256 : is_visible = .true.
257 : else
258 588240 : is_visible = .false.
259 : end if
260 :
261 866880 : end function is_visible
262 :
263 : !=========================================================================================
264 :
265 376888 : function get_molar_mass_ratio(gas_name) result(massratio)
266 :
267 : ! return the molar mass ratio of dry air to gas based on gas_name
268 :
269 : character(len=*),intent(in) :: gas_name
270 : real(r8) :: massratio
271 :
272 : ! local variables
273 : real(r8), parameter :: amdw = 1.607793_r8 ! Molecular weight of dry air / water vapor
274 : real(r8), parameter :: amdc = 0.658114_r8 ! Molecular weight of dry air / carbon dioxide
275 : real(r8), parameter :: amdo = 0.603428_r8 ! Molecular weight of dry air / ozone
276 : real(r8), parameter :: amdm = 1.805423_r8 ! Molecular weight of dry air / methane
277 : real(r8), parameter :: amdn = 0.658090_r8 ! Molecular weight of dry air / nitrous oxide
278 : real(r8), parameter :: amdo2 = 0.905140_r8 ! Molecular weight of dry air / oxygen
279 : real(r8), parameter :: amdc1 = 0.210852_r8 ! Molecular weight of dry air / CFC11
280 : real(r8), parameter :: amdc2 = 0.239546_r8 ! Molecular weight of dry air / CFC12
281 :
282 : character(len=*), parameter :: sub='get_molar_mass_ratio'
283 : !----------------------------------------------------------------------------
284 :
285 753776 : select case (trim(gas_name))
286 : case ('H2O')
287 47111 : massratio = amdw
288 : case ('CO2')
289 47111 : massratio = amdc
290 : case ('O3')
291 47111 : massratio = amdo
292 : case ('CH4')
293 47111 : massratio = amdm
294 : case ('N2O')
295 47111 : massratio = amdn
296 : case ('O2')
297 47111 : massratio = amdo2
298 : case ('CFC11')
299 47111 : massratio = amdc1
300 : case ('CFC12')
301 47111 : massratio = amdc2
302 : case default
303 753776 : call endrun(sub//": Invalid gas: "//trim(gas_name))
304 : end select
305 :
306 376888 : end function get_molar_mass_ratio
307 :
308 : !=========================================================================================
309 :
310 376888 : subroutine rad_gas_get_vmr(icall, gas_name, state, pbuf, nlay, numactivecols, gas_concs, idxday)
311 :
312 : ! Set volume mixing ratio in gas_concs object.
313 : ! The gas_concs%set_vmr method copies data into internally allocated storage.
314 :
315 : integer, intent(in) :: icall ! index of climate/diagnostic radiation call
316 : character(len=*), intent(in) :: gas_name
317 : type(physics_state), target, intent(in) :: state
318 : type(physics_buffer_desc), pointer :: pbuf(:)
319 : integer, intent(in) :: nlay ! number of layers in radiation calculation
320 : integer, intent(in) :: numactivecols ! number of columns, ncol for LW, nday for SW
321 :
322 : type(ty_gas_concs), intent(inout) :: gas_concs ! the result is VRM inside gas_concs
323 :
324 : integer, optional, intent(in) :: idxday(:) ! indices of daylight columns in a chunk
325 :
326 : ! Local variables
327 753776 : integer :: i, idx(numactivecols)
328 : integer :: istat
329 376888 : real(r8), pointer :: gas_mmr(:,:)
330 376888 : real(r8), allocatable :: gas_vmr(:,:)
331 376888 : real(r8), allocatable :: mmr(:,:)
332 : real(r8) :: massratio
333 :
334 : ! For ozone profile above model
335 : real(r8) :: P_top, P_int, P_mid, alpha, beta, a, b, chi_mid, chi_0, chi_eff
336 :
337 : character(len=128) :: errmsg
338 : character(len=*), parameter :: sub = 'rad_gas_get_vmr'
339 : !----------------------------------------------------------------------------
340 :
341 : ! set the column indices; when idxday is provided (e.g. daylit columns) use them, otherwise just count.
342 6208888 : do i = 1, numactivecols
343 6208888 : if (present(idxday)) then
344 1944000 : idx(i) = idxday(i)
345 : else
346 3888000 : idx(i) = i
347 : end if
348 : end do
349 :
350 : ! gas_mmr points to a "chunk" in either the state or pbuf objects. That storage is
351 : ! dimensioned (pcols,pver).
352 376888 : call rad_cnst_get_gas(icall, gas_name, state, pbuf, gas_mmr)
353 :
354 : ! Copy into storage for RRTMGP
355 1507552 : allocate(mmr(numactivecols, nlay), stat=istat)
356 376888 : call alloc_err(istat, sub, 'mmr', numactivecols*nlay)
357 1130664 : allocate(gas_vmr(numactivecols, nlay), stat=istat)
358 376888 : call alloc_err(istat, sub, 'gas_vmr', numactivecols*nlay)
359 :
360 6208888 : do i = 1, numactivecols
361 548584888 : mmr(i,ktoprad:) = gas_mmr(idx(i),ktopcam:)
362 : end do
363 :
364 : ! If an extra layer is being used, copy mmr from the top layer of CAM to the extra layer.
365 376888 : if (nlay == pverp) then
366 6208888 : mmr(:,1) = mmr(:,2)
367 : end if
368 :
369 : ! special case: H2O is specific humidity, not mixing ratio. Use r = q/(1-q):
370 376888 : if (gas_name == 'H2O') then
371 73001545 : mmr = mmr / (1._r8 - mmr)
372 : end if
373 :
374 : ! convert MMR to VMR, multipy by ratio of dry air molar mas to gas molar mass.
375 376888 : massratio = get_molar_mass_ratio(gas_name)
376 584389248 : gas_vmr = mmr * massratio
377 :
378 : ! special case: Setting O3 in the extra layer:
379 : !
380 : ! For the purpose of attenuating solar fluxes above the CAM model top, we assume that ozone
381 : ! mixing decreases linearly in each column from the value in the top layer of CAM to zero at
382 : ! the pressure level set by P_top. P_top has been set to 50 Pa (0.5 hPa) based on model tuning
383 : ! to produce temperatures at the top of CAM that are most consistent with WACCM at similar pressure levels.
384 :
385 376888 : if ((gas_name == 'O3') .and. (nlay == pverp)) then
386 : P_top = 50.0_r8
387 776111 : do i = 1, numactivecols
388 729000 : P_int = state%pint(idx(i),1) ! pressure (Pa) at upper interface of CAM
389 729000 : P_mid = state%pmid(idx(i),1) ! pressure (Pa) at midpoint of top layer of CAM
390 729000 : alpha = log(P_int/P_top)
391 729000 : beta = log(P_mid/P_int)/log(P_mid/P_top)
392 :
393 729000 : a = ( (1._r8 + alpha) * exp(-alpha) - 1._r8 ) / alpha
394 729000 : b = 1._r8 - exp(-alpha)
395 :
396 776111 : if (alpha .gt. 0) then ! only apply where top level is below 80 km
397 0 : chi_mid = gas_vmr(i,1) ! molar mixing ratio of O3 at midpoint of top layer
398 0 : chi_0 = chi_mid / (1._r8 + beta)
399 0 : chi_eff = chi_0 * (a + b)
400 0 : gas_vmr(i,1) = chi_eff
401 : end if
402 : end do
403 : end if
404 :
405 376888 : errmsg = gas_concs%set_vmr(gas_name, gas_vmr)
406 376888 : if (len_trim(errmsg) > 0) then
407 0 : call endrun(sub//': ERROR, gas_concs%set_vmr: '//trim(errmsg))
408 : end if
409 :
410 376888 : deallocate(gas_vmr)
411 376888 : deallocate(mmr)
412 :
413 753776 : end subroutine rad_gas_get_vmr
414 :
415 : !==================================================================================================
416 :
417 30960 : subroutine rrtmgp_set_gases_lw(icall, state, pbuf, nlay, gas_concs)
418 :
419 : ! Set gas vmr for the gases in the radconstants module's gaslist.
420 :
421 : ! The memory management for the gas_concs object is internal. The arrays passed to it
422 : ! are copied to the internally allocated memory. Each call to the set_vmr method checks
423 : ! whether the gas already has memory allocated, and if it does that memory is deallocated
424 : ! and new memory is allocated.
425 :
426 : ! arguments
427 : integer, intent(in) :: icall ! index of climate/diagnostic radiation call
428 : type(physics_state), target, intent(in) :: state
429 : type(physics_buffer_desc), pointer :: pbuf(:)
430 : integer, intent(in) :: nlay
431 : type(ty_gas_concs), intent(inout) :: gas_concs
432 :
433 : ! local variables
434 : integer :: i, ncol
435 : character(len=*), parameter :: sub = 'rrtmgp_set_gases_lw'
436 : !--------------------------------------------------------------------------------
437 :
438 30960 : ncol = state%ncol
439 278640 : do i = 1, nradgas
440 278640 : call rad_gas_get_vmr(icall, gaslist(i), state, pbuf, nlay, ncol, gas_concs)
441 : end do
442 30960 : end subroutine rrtmgp_set_gases_lw
443 :
444 : !==================================================================================================
445 :
446 16151 : subroutine rrtmgp_set_gases_sw( &
447 : icall, state, pbuf, nlay, nday, &
448 16151 : idxday, gas_concs)
449 :
450 : ! Return gas_concs with gas volume mixing ratio on DAYLIT columns.
451 : ! Set all gases in radconstants gaslist.
452 :
453 : ! arguments
454 : integer, intent(in) :: icall ! index of climate/diagnostic radiation call
455 : type(physics_state), target, intent(in) :: state
456 : type(physics_buffer_desc), pointer :: pbuf(:)
457 : integer, intent(in) :: nlay
458 : integer, intent(in) :: nday
459 : integer, intent(in) :: idxday(:)
460 : type(ty_gas_concs), intent(inout) :: gas_concs
461 :
462 : ! local variables
463 : integer :: i
464 : character(len=*), parameter :: sub = 'rrtmgp_set_gases_sw'
465 : !----------------------------------------------------------------------------
466 :
467 : ! use the optional argument idxday to specify which columns are sunlit
468 145359 : do i = 1,nradgas
469 145359 : call rad_gas_get_vmr(icall, gaslist(i), state, pbuf, nlay, nday, gas_concs, idxday=idxday)
470 : end do
471 :
472 16151 : end subroutine rrtmgp_set_gases_sw
473 :
474 : !==================================================================================================
475 :
476 30960 : subroutine rrtmgp_set_cloud_lw( &
477 : state, pbuf, ncol, nlay, nlaycam, &
478 : cld, cldfsnow, cldfgrau, cldfprime, graupel_in_rad, &
479 : kdist_lw, cloud_lw, cld_lw_abs_cloudsim, snow_lw_abs_cloudsim, grau_lw_abs_cloudsim)
480 :
481 : ! Compute combined cloud optical properties.
482 : ! Create MCICA stochastic arrays for cloud LW optical properties.
483 : ! Initialize optical properties object (cloud_lw) and load with MCICA columns.
484 :
485 : ! arguments
486 : type(physics_state), intent(in) :: state
487 : type(physics_buffer_desc), pointer :: pbuf(:)
488 : integer, intent(in) :: ncol ! number of columns in CAM chunk
489 : integer, intent(in) :: nlay ! number of layers in radiation calculation (may include "extra layer")
490 : integer, intent(in) :: nlaycam ! number of CAM layers in radiation calculation
491 : real(r8), pointer :: cld(:,:) ! cloud fraction (liq+ice)
492 : real(r8), pointer :: cldfsnow(:,:) ! cloud fraction of just "snow clouds"
493 : real(r8), pointer :: cldfgrau(:,:) ! cloud fraction of just "graupel clouds"
494 : real(r8), intent(in) :: cldfprime(pcols,pver) ! combined cloud fraction
495 :
496 : logical, intent(in) :: graupel_in_rad ! use graupel in radiation code
497 : class(ty_gas_optics_rrtmgp), intent(in) :: kdist_lw
498 : type(ty_optical_props_1scl), intent(out) :: cloud_lw
499 :
500 : ! Diagnostic outputs
501 : real(r8), intent(out) :: cld_lw_abs_cloudsim(pcols,pver) ! in-cloud liq+ice optical depth (for COSP)
502 : real(r8), intent(out) :: snow_lw_abs_cloudsim(pcols,pver)! in-cloud snow optical depth (for COSP)
503 : real(r8), intent(out) :: grau_lw_abs_cloudsim(pcols,pver)! in-cloud Graupel optical depth (for COSP)
504 :
505 : ! Local variables
506 :
507 : integer :: i, k
508 :
509 : ! cloud radiative parameters are "in cloud" not "in cell"
510 : real(r8) :: liq_lw_abs(nlwbands,pcols,pver) ! liquid absorption optics depth (LW)
511 : real(r8) :: ice_lw_abs(nlwbands,pcols,pver) ! ice absorption optics depth (LW)
512 : real(r8) :: cld_lw_abs(nlwbands,pcols,pver) ! cloud absorption optics depth (LW)
513 : real(r8) :: snow_lw_abs(nlwbands,pcols,pver) ! snow absorption optics depth (LW)
514 : real(r8) :: grau_lw_abs(nlwbands,pcols,pver) ! graupel absorption optics depth (LW)
515 : real(r8) :: c_cld_lw_abs(nlwbands,pcols,pver) ! combined cloud absorption optics depth (LW)
516 :
517 : ! Arrays for converting from CAM chunks to RRTMGP inputs.
518 61920 : real(r8) :: cldf(ncol,nlaycam)
519 61920 : real(r8) :: tauc(nlwbands,ncol,nlaycam)
520 30960 : real(r8) :: taucmcl(nlwgpts,ncol,nlaycam)
521 :
522 : character(len=128) :: errmsg
523 : character(len=*), parameter :: sub = 'rrtmgp_set_cloud_lw'
524 : !--------------------------------------------------------------------------------
525 :
526 : ! Combine the cloud optical properties. These calculations are done on CAM "chunks".
527 :
528 : ! gammadist liquid optics
529 30960 : call liquid_cloud_get_rad_props_lw(state, pbuf, liq_lw_abs)
530 : ! Mitchell ice optics
531 30960 : call ice_cloud_get_rad_props_lw(state, pbuf, ice_lw_abs)
532 :
533 771276240 : cld_lw_abs(:,:ncol,:) = liq_lw_abs(:,:ncol,:) + ice_lw_abs(:,:ncol,:)
534 :
535 30960 : if (associated(cldfsnow)) then
536 : ! add in snow
537 30960 : call snow_cloud_get_rad_props_lw(state, pbuf, snow_lw_abs)
538 516960 : do i = 1, ncol
539 45714960 : do k = 1, pver
540 45684000 : if (cldfprime(i,k) > 0._r8) then
541 0 : c_cld_lw_abs(:,i,k) = ( cldfsnow(i,k)*snow_lw_abs(:,i,k) &
542 143353027 : + cld(i,k)*cld_lw_abs(:,i,k) )/cldfprime(i,k)
543 : else
544 625012973 : c_cld_lw_abs(:,i,k) = 0._r8
545 : end if
546 : end do
547 : end do
548 : else
549 0 : c_cld_lw_abs(:,:ncol,:) = cld_lw_abs(:,:ncol,:)
550 : end if
551 :
552 : ! add in graupel
553 30960 : if (associated(cldfgrau) .and. graupel_in_rad) then
554 0 : call grau_cloud_get_rad_props_lw(state, pbuf, grau_lw_abs)
555 0 : do i = 1, ncol
556 0 : do k = 1, pver
557 0 : if (cldfprime(i,k) > 0._r8) then
558 0 : c_cld_lw_abs(:,i,k) = ( cldfgrau(i,k)*grau_lw_abs(:,i,k) &
559 0 : + cld(i,k)*c_cld_lw_abs(:,i,k) )/cldfprime(i,k)
560 : else
561 0 : c_cld_lw_abs(:,i,k) = 0._r8
562 : end if
563 : end do
564 : end do
565 : end if
566 :
567 : ! Cloud optics for COSP
568 48978720 : cld_lw_abs_cloudsim = cld_lw_abs(idx_lw_cloudsim,:,:)
569 48978720 : snow_lw_abs_cloudsim = snow_lw_abs(idx_lw_cloudsim,:,:)
570 48978720 : grau_lw_abs_cloudsim = grau_lw_abs(idx_lw_cloudsim,:,:)
571 :
572 : ! Extract just the layers of CAM where RRTMGP does calculations.
573 :
574 : ! Subset "chunk" data so just the number of CAM layers in the
575 : ! radiation calculation are used by MCICA to produce subcolumns.
576 48108240 : cldf = cldfprime(:ncol, ktopcam:)
577 771276240 : tauc = c_cld_lw_abs(:, :ncol, ktopcam:)
578 :
579 : ! Enforce tauc >= 0.
580 771276240 : tauc = merge(tauc, 0.0_r8, tauc > 0.0_r8)
581 :
582 : ! MCICA uses spectral data (on bands) to construct subcolumns (one per g-point)
583 : call mcica_subcol_lw( &
584 : kdist_lw, nlwbands, nlwgpts, ncol, nlaycam, &
585 30960 : nlwgpts, state%pmid, cldf, tauc, taucmcl )
586 :
587 30960 : errmsg =cloud_lw%alloc_1scl(ncol, nlay, kdist_lw)
588 30960 : if (len_trim(errmsg) > 0) then
589 0 : call endrun(sub//': ERROR: cloud_lw%alloc_1scalar: '//trim(errmsg))
590 : end if
591 :
592 : ! If there is an extra layer in the radiation then this initialization
593 : ! will provide zero optical depths there.
594 6224056560 : cloud_lw%tau = 0.0_r8
595 :
596 : ! Set the properties on g-points.
597 3993840 : do i = 1, nlwgpts
598 6157885680 : cloud_lw%tau(:,ktoprad:,i) = taucmcl(i,:,:)
599 : end do
600 :
601 : ! validate checks that: tau > 0
602 30960 : errmsg = cloud_lw%validate()
603 30960 : if (len_trim(errmsg) > 0) then
604 0 : call endrun(sub//': ERROR: cloud_lw%validate: '//trim(errmsg))
605 : end if
606 :
607 30960 : end subroutine rrtmgp_set_cloud_lw
608 :
609 : !==================================================================================================
610 :
611 30960 : subroutine rrtmgp_set_cloud_sw( &
612 : state, pbuf, nlay, nday, idxday, &
613 30960 : nnite, idxnite, pmid, cld, cldfsnow, &
614 : cldfgrau, cldfprime, graupel_in_rad, kdist_sw, cloud_sw, &
615 : tot_cld_vistau, tot_icld_vistau, liq_icld_vistau, ice_icld_vistau, snow_icld_vistau, &
616 : grau_icld_vistau, cld_tau_cloudsim, snow_tau_cloudsim, grau_tau_cloudsim)
617 :
618 : ! Compute combined cloud optical properties.
619 : ! Create MCICA stochastic arrays for cloud SW optical properties.
620 : ! Initialize optical properties object (cloud_sw) and load with MCICA columns.
621 :
622 : ! arguments
623 : type(physics_state), target, intent(in) :: state
624 : type(physics_buffer_desc), pointer :: pbuf(:)
625 : integer, intent(in) :: nlay ! number of layers in radiation calculation (may include "extra layer")
626 : integer, intent(in) :: nday ! number of daylight columns
627 : integer, intent(in) :: idxday(pcols) ! indices of daylight columns in the chunk
628 : integer, intent(in) :: nnite ! number of night columns
629 : integer, intent(in) :: idxnite(pcols) ! indices of night columns in the chunk
630 :
631 : real(r8), intent(in) :: pmid(nday,nlay)! pressure at layer midpoints (Pa) used to seed RNG.
632 :
633 : real(r8), pointer :: cld(:,:) ! cloud fraction (liq+ice)
634 : real(r8), pointer :: cldfsnow(:,:) ! cloud fraction of just "snow clouds"
635 : real(r8), pointer :: cldfgrau(:,:) ! cloud fraction of just "graupel clouds"
636 : real(r8), intent(in) :: cldfprime(pcols,pver) ! combined cloud fraction
637 :
638 : logical, intent(in) :: graupel_in_rad ! graupel in radiation code
639 : class(ty_gas_optics_rrtmgp), intent(in) :: kdist_sw ! shortwave gas optics object
640 : type(ty_optical_props_2str), intent(out) :: cloud_sw ! SW cloud optical properties object
641 :
642 : ! Diagnostic outputs
643 : real(r8), intent(out) :: tot_cld_vistau(pcols,pver) ! gbx total cloud optical depth
644 : real(r8), intent(out) :: tot_icld_vistau(pcols,pver) ! in-cld total cloud optical depth
645 : real(r8), intent(out) :: liq_icld_vistau(pcols,pver) ! in-cld liq cloud optical depth
646 : real(r8), intent(out) :: ice_icld_vistau(pcols,pver) ! in-cld ice cloud optical depth
647 : real(r8), intent(out) :: snow_icld_vistau(pcols,pver) ! snow in-cloud visible sw optical depth
648 : real(r8), intent(out) :: grau_icld_vistau(pcols,pver) ! Graupel in-cloud visible sw optical depth
649 : real(r8), intent(out) :: cld_tau_cloudsim(pcols,pver) ! in-cloud liq+ice optical depth (for COSP)
650 : real(r8), intent(out) :: snow_tau_cloudsim(pcols,pver)! in-cloud snow optical depth (for COSP)
651 : real(r8), intent(out) :: grau_tau_cloudsim(pcols,pver)! in-cloud Graupel optical depth (for COSP)
652 :
653 : ! Local variables
654 :
655 : integer :: i, k, ncol
656 : integer :: igpt, nver
657 : integer :: istat
658 : integer, parameter :: changeseed = 1
659 :
660 : ! cloud radiative parameters are "in cloud" not "in cell"
661 : real(r8) :: ice_tau (nswbands,pcols,pver) ! ice extinction optical depth
662 : real(r8) :: ice_tau_w (nswbands,pcols,pver) ! ice single scattering albedo * tau
663 : real(r8) :: ice_tau_w_g(nswbands,pcols,pver) ! ice asymmetry parameter * tau * w
664 : real(r8) :: liq_tau (nswbands,pcols,pver) ! liquid extinction optical depth
665 : real(r8) :: liq_tau_w (nswbands,pcols,pver) ! liquid single scattering albedo * tau
666 : real(r8) :: liq_tau_w_g(nswbands,pcols,pver) ! liquid asymmetry parameter * tau * w
667 : real(r8) :: cld_tau (nswbands,pcols,pver) ! cloud extinction optical depth
668 : real(r8) :: cld_tau_w (nswbands,pcols,pver) ! cloud single scattering albedo * tau
669 : real(r8) :: cld_tau_w_g(nswbands,pcols,pver) ! cloud asymmetry parameter * w * tau
670 : real(r8) :: snow_tau (nswbands,pcols,pver) ! snow extinction optical depth
671 : real(r8) :: snow_tau_w (nswbands,pcols,pver) ! snow single scattering albedo * tau
672 : real(r8) :: snow_tau_w_g(nswbands,pcols,pver) ! snow asymmetry parameter * tau * w
673 : real(r8) :: grau_tau (nswbands,pcols,pver) ! graupel extinction optical depth
674 : real(r8) :: grau_tau_w (nswbands,pcols,pver) ! graupel single scattering albedo * tau
675 : real(r8) :: grau_tau_w_g(nswbands,pcols,pver) ! graupel asymmetry parameter * tau * w
676 : real(r8) :: c_cld_tau (nswbands,pcols,pver) ! combined cloud extinction optical depth
677 : real(r8) :: c_cld_tau_w (nswbands,pcols,pver) ! combined cloud single scattering albedo * tau
678 : real(r8) :: c_cld_tau_w_g(nswbands,pcols,pver) ! combined cloud asymmetry parameter * w * tau
679 :
680 : ! RRTMGP does not use this property in its 2-stream calculations.
681 : real(r8) :: sw_tau_w_f(nswbands,pcols,pver) ! Forward scattered fraction * tau * w.
682 :
683 : ! Arrays for converting from CAM chunks to RRTMGP inputs.
684 30960 : real(r8), allocatable :: cldf(:,:)
685 30960 : real(r8), allocatable :: tauc(:,:,:)
686 30960 : real(r8), allocatable :: ssac(:,:,:)
687 30960 : real(r8), allocatable :: asmc(:,:,:)
688 30960 : real(r8), allocatable :: taucmcl(:,:,:)
689 30960 : real(r8), allocatable :: ssacmcl(:,:,:)
690 30960 : real(r8), allocatable :: asmcmcl(:,:,:)
691 30960 : real(r8), allocatable :: day_cld_tau(:,:,:)
692 30960 : real(r8), allocatable :: day_cld_tau_w(:,:,:)
693 30960 : real(r8), allocatable :: day_cld_tau_w_g(:,:,:)
694 :
695 : character(len=128) :: errmsg
696 : character(len=*), parameter :: sub = 'rrtmgp_set_cloud_sw'
697 : !--------------------------------------------------------------------------------
698 :
699 30960 : ncol = state%ncol
700 :
701 : ! Combine the cloud optical properties. These calculations are done on CAM "chunks".
702 :
703 : ! gammadist liquid optics
704 30960 : call get_liquid_optics_sw(state, pbuf, liq_tau, liq_tau_w, liq_tau_w_g, sw_tau_w_f)
705 : ! Mitchell ice optics
706 30960 : call get_ice_optics_sw(state, pbuf, ice_tau, ice_tau_w, ice_tau_w_g, sw_tau_w_f)
707 :
708 680880240 : cld_tau(:,:ncol,:) = liq_tau(:,:ncol,:) + ice_tau(:,:ncol,:)
709 680880240 : cld_tau_w(:,:ncol,:) = liq_tau_w(:,:ncol,:) + ice_tau_w(:,:ncol,:)
710 680880240 : cld_tau_w_g(:,:ncol,:) = liq_tau_w_g(:,:ncol,:) + ice_tau_w_g(:,:ncol,:)
711 :
712 : ! add in snow
713 30960 : if (associated(cldfsnow)) then
714 30960 : call get_snow_optics_sw(state, pbuf, snow_tau, snow_tau_w, snow_tau_w_g, sw_tau_w_f)
715 516960 : do i = 1, ncol
716 45714960 : do k = 1, pver
717 45684000 : if (cldfprime(i,k) > 0._r8) then
718 0 : c_cld_tau(:,i,k) = ( cldfsnow(i,k)*snow_tau(:,i,k) &
719 126487965 : + cld(i,k)*cld_tau(:,i,k) )/cldfprime(i,k)
720 0 : c_cld_tau_w(:,i,k) = ( cldfsnow(i,k)*snow_tau_w(:,i,k) &
721 126487965 : + cld(i,k)*cld_tau_w(:,i,k) )/cldfprime(i,k)
722 0 : c_cld_tau_w_g(:,i,k) = ( cldfsnow(i,k)*snow_tau_w_g(:,i,k) &
723 126487965 : + cld(i,k)*cld_tau_w_g(:,i,k) )/cldfprime(i,k)
724 : else
725 551482035 : c_cld_tau(:,i,k) = 0._r8
726 551482035 : c_cld_tau_w(:,i,k) = 0._r8
727 551482035 : c_cld_tau_w_g(:,i,k) = 0._r8
728 : end if
729 : end do
730 : end do
731 : else
732 0 : c_cld_tau(:,:ncol,:) = cld_tau(:,:ncol,:)
733 0 : c_cld_tau_w(:,:ncol,:) = cld_tau_w(:,:ncol,:)
734 0 : c_cld_tau_w_g(:,:ncol,:) = cld_tau_w_g(:,:ncol,:)
735 : end if
736 :
737 : ! add in graupel
738 30960 : if (associated(cldfgrau) .and. graupel_in_rad) then
739 0 : call get_grau_optics_sw(state, pbuf, grau_tau, grau_tau_w, grau_tau_w_g, sw_tau_w_f)
740 0 : do i = 1, ncol
741 0 : do k = 1, pver
742 0 : if (cldfprime(i,k) > 0._r8) then
743 0 : c_cld_tau(:,i,k) = ( cldfgrau(i,k)*grau_tau(:,i,k) &
744 0 : + cld(i,k)*c_cld_tau(:,i,k) )/cldfprime(i,k)
745 0 : c_cld_tau_w(:,i,k) = ( cldfgrau(i,k)*grau_tau_w(:,i,k) &
746 0 : + cld(i,k)*c_cld_tau_w(:,i,k) )/cldfprime(i,k)
747 0 : c_cld_tau_w_g(:,i,k) = ( cldfgrau(i,k)*grau_tau_w_g(:,i,k) &
748 0 : + cld(i,k)*c_cld_tau_w_g(:,i,k) )/cldfprime(i,k)
749 : else
750 0 : c_cld_tau(:,i,k) = 0._r8
751 0 : c_cld_tau_w(:,i,k) = 0._r8
752 0 : c_cld_tau_w_g(:,i,k) = 0._r8
753 : end if
754 : end do
755 : end do
756 : end if
757 :
758 : ! cloud optical properties need to be re-ordered from the RRTMG spectral bands
759 : ! (assumed in the optics datasets) to RRTMGP's
760 1361729520 : ice_tau(:,:ncol,:) = ice_tau(rrtmg_to_rrtmgp_swbands,:ncol,:)
761 1361729520 : liq_tau(:,:ncol,:) = liq_tau(rrtmg_to_rrtmgp_swbands,:ncol,:)
762 1361729520 : c_cld_tau(:,:ncol,:) = c_cld_tau(rrtmg_to_rrtmgp_swbands,:ncol,:)
763 1361729520 : c_cld_tau_w(:,:ncol,:) = c_cld_tau_w(rrtmg_to_rrtmgp_swbands,:ncol,:)
764 1361729520 : c_cld_tau_w_g(:,:ncol,:) = c_cld_tau_w_g(rrtmg_to_rrtmgp_swbands,:ncol,:)
765 30960 : if (associated(cldfsnow)) then
766 1361729520 : snow_tau(:,:ncol,:) = snow_tau(rrtmg_to_rrtmgp_swbands,:ncol,:)
767 : end if
768 30960 : if (associated(cldfgrau) .and. graupel_in_rad) then
769 0 : grau_tau(:,:ncol,:) = grau_tau(rrtmg_to_rrtmgp_swbands,:ncol,:)
770 : end if
771 :
772 : ! Set arrays for diagnostic output.
773 : ! cloud optical depth fields for the visible band
774 48108240 : tot_icld_vistau(:ncol,:) = c_cld_tau(idx_sw_diag,:ncol,:)
775 48108240 : liq_icld_vistau(:ncol,:) = liq_tau(idx_sw_diag,:ncol,:)
776 48108240 : ice_icld_vistau(:ncol,:) = ice_tau(idx_sw_diag,:ncol,:)
777 30960 : if (associated(cldfsnow)) then
778 48108240 : snow_icld_vistau(:ncol,:) = snow_tau(idx_sw_diag,:ncol,:)
779 : endif
780 30960 : if (associated(cldfgrau) .and. graupel_in_rad) then
781 0 : grau_icld_vistau(:ncol,:) = grau_tau(idx_sw_diag,:ncol,:)
782 : endif
783 :
784 : ! multiply by total cloud fraction to get gridbox value
785 48108240 : tot_cld_vistau(:ncol,:) = c_cld_tau(idx_sw_diag,:ncol,:)*cldfprime(:ncol,:)
786 :
787 : ! overwrite night columns with fillvalue
788 273960 : do i = 1, Nnite
789 22842000 : tot_cld_vistau(IdxNite(i),:) = fillvalue
790 22842000 : tot_icld_vistau(IdxNite(i),:) = fillvalue
791 22842000 : liq_icld_vistau(IdxNite(i),:) = fillvalue
792 22842000 : ice_icld_vistau(IdxNite(i),:) = fillvalue
793 243000 : if (associated(cldfsnow)) then
794 22842000 : snow_icld_vistau(IdxNite(i),:) = fillvalue
795 : end if
796 273960 : if (associated(cldfgrau) .and. graupel_in_rad) then
797 0 : grau_icld_vistau(IdxNite(i),:) = fillvalue
798 : end if
799 : end do
800 :
801 : ! Cloud optics for COSP
802 48978720 : cld_tau_cloudsim = cld_tau(idx_sw_cloudsim,:,:)
803 48978720 : snow_tau_cloudsim = snow_tau(idx_sw_cloudsim,:,:)
804 48978720 : grau_tau_cloudsim = grau_tau(idx_sw_cloudsim,:,:)
805 :
806 : ! if no daylight columns the cloud_sw object isn't initialized
807 30960 : if (nday > 0) then
808 :
809 : ! number of CAM's layers in radiation calculation. Does not include the "extra layer".
810 16151 : nver = pver - ktopcam + 1
811 :
812 : allocate( &
813 : cldf(nday,nver), &
814 : day_cld_tau(nswbands,nday,nver), &
815 : day_cld_tau_w(nswbands,nday,nver), &
816 : day_cld_tau_w_g(nswbands,nday,nver), &
817 : tauc(nswbands,nday,nver), taucmcl(nswgpts,nday,nver), &
818 : ssac(nswbands,nday,nver), ssacmcl(nswgpts,nday,nver), &
819 436077 : asmc(nswbands,nday,nver), asmcmcl(nswgpts,nday,nver), stat=istat)
820 16151 : call alloc_err(istat, sub, 'cldf,..,asmcmcl', 9*nswgpts*nday*nver)
821 :
822 : ! Subset "chunk" data so just the daylight columns, and the number of CAM layers in the
823 : ! radiation calculation are used by MCICA to produce subcolumns.
824 24133345 : cldf = cldfprime( idxday(1:nday), ktopcam:)
825 340519345 : day_cld_tau = c_cld_tau( :, idxday(1:nday), ktopcam:)
826 340519345 : day_cld_tau_w = c_cld_tau_w( :, idxday(1:nday), ktopcam:)
827 340519345 : day_cld_tau_w_g = c_cld_tau_w_g(:, idxday(1:nday), ktopcam:)
828 :
829 : ! Compute the optical properties needed for the 2-stream calculations. These calculations
830 : ! are the same as the RRTMG version.
831 :
832 : ! set cloud optical depth, clip @ zero
833 340519345 : tauc = merge(day_cld_tau, 0.0_r8, day_cld_tau > 0.0_r8)
834 : ! set value of asymmetry
835 340519345 : asmc = merge(day_cld_tau_w_g / max(day_cld_tau_w, tiny), 0.0_r8, day_cld_tau_w > 0.0_r8)
836 : ! set value of single scattering albedo
837 340519345 : ssac = merge(max(day_cld_tau_w, tiny) / max(tauc, tiny), 1.0_r8 , tauc > 0.0_r8)
838 : ! set asymmetry to zero when tauc = 0
839 340519345 : asmc = merge(asmc, 0.0_r8, tauc > 0.0_r8)
840 :
841 : ! MCICA uses spectral data (on bands) to construct subcolumns (one per g-point)
842 : call mcica_subcol_sw( &
843 : kdist_sw, nswbands, nswgpts, nday, nlay, &
844 : nver, changeseed, pmid, cldf, tauc, &
845 16151 : ssac, asmc, taucmcl, ssacmcl, asmcmcl)
846 :
847 : ! Initialize object for SW cloud optical properties.
848 16151 : errmsg = cloud_sw%alloc_2str(nday, nlay, kdist_sw)
849 16151 : if (len_trim(errmsg) > 0) then
850 0 : call endrun(trim(sub)//': ERROR: cloud_sw%alloc_2str: '//trim(errmsg))
851 : end if
852 :
853 : ! If there is an extra layer in the radiation then this initialization
854 : ! will provide the optical properties there.
855 2730166791 : cloud_sw%tau = 0.0_r8
856 2730166791 : cloud_sw%ssa = 1.0_r8
857 2730166791 : cloud_sw%g = 0.0_r8
858 :
859 : ! Set the properties on g-points.
860 1825063 : do igpt = 1,nswgpts
861 2701125728 : cloud_sw%g (:,ktoprad:,igpt) = asmcmcl(igpt,:,:)
862 2701125728 : cloud_sw%ssa(:,ktoprad:,igpt) = ssacmcl(igpt,:,:)
863 2701141879 : cloud_sw%tau(:,ktoprad:,igpt) = taucmcl(igpt,:,:)
864 : end do
865 :
866 : ! validate checks that: tau > 0, ssa is in range [0,1], and g is in range [-1,1].
867 16151 : errmsg = cloud_sw%validate()
868 16151 : if (len_trim(errmsg) > 0) then
869 0 : call endrun(sub//': ERROR: cloud_sw%validate: '//trim(errmsg))
870 : end if
871 :
872 : ! delta scaling adjusts for forward scattering
873 16151 : errmsg = cloud_sw%delta_scale()
874 16151 : if (len_trim(errmsg) > 0) then
875 0 : call endrun(sub//': ERROR: cloud_sw%delta_scale: '//trim(errmsg))
876 : end if
877 :
878 : ! All information is in cloud_sw, now deallocate local vars.
879 0 : deallocate( &
880 0 : cldf, tauc, ssac, asmc, &
881 0 : taucmcl, ssacmcl, asmcmcl,&
882 16151 : day_cld_tau, day_cld_tau_w, day_cld_tau_w_g )
883 :
884 : end if
885 :
886 30960 : end subroutine rrtmgp_set_cloud_sw
887 :
888 : !==================================================================================================
889 :
890 30960 : subroutine rrtmgp_set_aer_lw(icall, state, pbuf, aer_lw)
891 :
892 : ! Load LW aerosol optical properties into the RRTMGP object.
893 :
894 : ! Arguments
895 : integer, intent(in) :: icall
896 : type(physics_state), target, intent(in) :: state
897 : type(physics_buffer_desc), pointer :: pbuf(:)
898 :
899 : type(ty_optical_props_1scl), intent(inout) :: aer_lw
900 :
901 : ! Local variables
902 : integer :: ncol
903 :
904 : ! Aerosol LW absorption optical depth
905 : real(r8) :: aer_lw_abs (pcols,pver,nlwbands)
906 :
907 : character(len=*), parameter :: sub = 'rrtmgp_set_aer_lw'
908 : character(len=128) :: errmsg
909 : !--------------------------------------------------------------------------------
910 :
911 30960 : ncol = state%ncol
912 :
913 : ! Get aerosol longwave optical properties.
914 30960 : call aer_rad_props_lw(icall, state, pbuf, aer_lw_abs)
915 :
916 : ! If there is an extra layer in the radiation then this initialization
917 : ! will provide zero optical depths there.
918 778034160 : aer_lw%tau = 0.0_r8
919 :
920 769762800 : aer_lw%tau(:ncol, ktoprad:, :) = aer_lw_abs(:ncol, ktopcam:, :)
921 :
922 30960 : errmsg = aer_lw%validate()
923 30960 : if (len_trim(errmsg) > 0) then
924 0 : call endrun(sub//': ERROR: aer_lw%validate: '//trim(errmsg))
925 : end if
926 30960 : end subroutine rrtmgp_set_aer_lw
927 :
928 : !==================================================================================================
929 :
930 30960 : subroutine rrtmgp_set_aer_sw( &
931 30960 : icall, state, pbuf, nday, idxday, nnite, idxnite, aer_sw)
932 :
933 : ! Load SW aerosol optical properties into the RRTMGP object.
934 :
935 : ! Arguments
936 : integer, intent(in) :: icall
937 : type(physics_state), target, intent(in) :: state
938 : type(physics_buffer_desc), pointer :: pbuf(:)
939 : integer, intent(in) :: nday
940 : integer, intent(in) :: idxday(:)
941 : integer, intent(in) :: nnite ! number of night columns
942 : integer, intent(in) :: idxnite(pcols) ! indices of night columns in the chunk
943 :
944 : type(ty_optical_props_2str), intent(inout) :: aer_sw
945 :
946 : ! local variables
947 : integer :: i
948 :
949 : ! The optical arrays dimensioned in the vertical as 0:pver.
950 : ! The index 0 is for the extra layer used in the radiation
951 : ! calculation. The index ktopcam assumes the CAM vertical indices are
952 : ! in the range 1:pver, so using this index correctly ignores vertical
953 : ! index 0. If an "extra" layer is used in the calculations, it is
954 : ! provided and set in the RRTMGP aerosol object aer_sw.
955 : real(r8) :: aer_tau (pcols,0:pver,nswbands) ! extinction optical depth
956 : real(r8) :: aer_tau_w (pcols,0:pver,nswbands) ! single scattering albedo * tau
957 : real(r8) :: aer_tau_w_g(pcols,0:pver,nswbands) ! asymmetry parameter * w * tau
958 : real(r8) :: aer_tau_w_f(pcols,0:pver,nswbands) ! forward scattered fraction * w * tau
959 : ! aer_tau_w_f is not used by RRTMGP.
960 : character(len=*), parameter :: sub = 'rrtmgp_set_aer_sw'
961 : !--------------------------------------------------------------------------------
962 :
963 : ! Get aerosol shortwave optical properties.
964 : ! Make outfld calls for aerosol optical property diagnostics.
965 : call aer_rad_props_sw( &
966 : icall, state, pbuf, nnite, idxnite, &
967 30960 : aer_tau, aer_tau_w, aer_tau_w_g, aer_tau_w_f)
968 :
969 : ! The aer_sw object is only initialized if nday > 0.
970 30960 : if (nday > 0) then
971 :
972 : ! aerosol optical properties need to be re-ordered from the RRTMG spectral bands
973 : ! (as assumed in the optics datasets) to the RRTMGP band order.
974 723144874 : aer_tau(:,:,:) = aer_tau( :,:,rrtmg_to_rrtmgp_swbands)
975 723144874 : aer_tau_w(:,:,:) = aer_tau_w( :,:,rrtmg_to_rrtmgp_swbands)
976 723144874 : aer_tau_w_g(:,:,:) = aer_tau_w_g(:,:,rrtmg_to_rrtmgp_swbands)
977 :
978 : ! If there is an extra layer in the radiation then this initialization
979 : ! will provide default values.
980 341284981 : aer_sw%tau = 0.0_r8
981 341284981 : aer_sw%ssa = 1.0_r8
982 341284981 : aer_sw%g = 0.0_r8
983 :
984 : ! CAM fields are products tau, tau*ssa, tau*ssa*asy
985 : ! Fields expected by RRTMGP are computed by
986 : ! aer_sw%tau = aer_tau
987 : ! aer_sw%ssa = aer_tau_w / aer_tau
988 : ! aer_sw%g = aer_tau_w_g / aer_taw_w
989 : ! aer_sw arrays have dimensions of (nday,nlay,nswbands)
990 :
991 259151 : do i = 1, nday
992 : ! set aerosol optical depth, clip to zero
993 320031000 : aer_sw%tau(i,ktoprad:,:) = max(aer_tau(idxday(i),ktopcam:,:), 0._r8)
994 : ! set value of single scattering albedo
995 243000 : aer_sw%ssa(i,ktoprad:,:) = merge(aer_tau_w(idxday(i),ktopcam:,:)/aer_tau(idxday(i),ktopcam:,:), &
996 320274000 : 1._r8, aer_tau(idxday(i),ktopcam:,:) > 0._r8)
997 : ! set value of asymmetry
998 243000 : aer_sw%g(i,ktoprad:,:) = merge(aer_tau_w_g(idxday(i),ktopcam:,:)/aer_tau_w(idxday(i),ktopcam:,:), &
999 320290151 : 0._r8, aer_tau_w(idxday(i),ktopcam:,:) > tiny)
1000 : end do
1001 :
1002 : ! impose limits on the components
1003 341284981 : aer_sw%ssa = min(max(aer_sw%ssa, 0._r8), 1._r8)
1004 341284981 : aer_sw%g = min(max(aer_sw%g, -1._r8), 1._r8)
1005 :
1006 : end if
1007 :
1008 30960 : end subroutine rrtmgp_set_aer_sw
1009 :
1010 : !==================================================================================================
1011 :
1012 : end module rrtmgp_inputs
|
