Line data Source code
1 :
2 : module radheat
3 : !-----------------------------------------------------------------------
4 : !
5 : ! Purpose: Provide an interface to convert shortwave and longwave
6 : ! radiative heating terms into net heating.
7 : !
8 : ! This module provides a hook to allow incorporating additional
9 : ! radiative terms (eUV heating and nonLTE longwave cooling).
10 : !
11 : ! Original version: B.A. Boville
12 : ! Change weighting function for RRTMG: A J Conley
13 : !-----------------------------------------------------------------------
14 :
15 : ! Use a cubic polynomial over the domain from minimum pressure to maximum pressure
16 : ! Cubic polynomial is chosen so that derivative is zero at minimum and maximum pressures
17 : ! and is monotonically increasing from zero at minimum pressure to one at maximum pressure
18 :
19 : use shr_kind_mod, only: r8 => shr_kind_r8
20 : use spmd_utils, only: masterproc
21 : use ppgrid, only: pcols, pver
22 : use physics_types, only: physics_state, physics_ptend, physics_ptend_init
23 : use physconst, only: gravit
24 : use air_composition, only: cpairv
25 : use perf_mod
26 : use cam_logfile, only: iulog
27 :
28 : implicit none
29 : private
30 : save
31 :
32 : ! Public interfaces
33 : public &
34 : radheat_readnl, &!
35 : radheat_register, &!
36 : radheat_init, &!
37 : radheat_timestep_init, &!
38 : radheat_tend ! return net radiative heating
39 :
40 : public :: radheat_disable_waccm ! disable waccm heating in the upper atm
41 :
42 : ! Namelist variables
43 : logical :: nlte_use_mo = .true. ! Determines which constituents are used from NLTE calculations
44 : ! = .true. uses prognostic constituents
45 : ! = .false. uses constituents from prescribed dataset waccm_forcing_file
46 : logical :: nlte_limit_co2 = .false. ! if true apply upper limit to co2 in the Fomichev scheme
47 : logical :: nlte_use_aliarms = .false. ! If true, use ALI-ARMS for the cooling rate calculation
48 : integer :: nlte_aliarms_every_X = 1 ! Call aliarms every X times radiation is called
49 :
50 : ! Private variables for merging heating rates
51 : real(r8):: qrs_wt(pver) ! merge weight for cam solar heating
52 : real(r8):: qrl_wt(pver) ! merge weight for cam long wave heating
53 :
54 : logical :: waccm_heating
55 : logical :: waccm_heating_on = .true.
56 :
57 : ! sw merge region
58 : ! highest altitude (lowest pressure) of merge region (Pa)
59 : real(r8) :: min_pressure_sw= 5._r8
60 : ! lowest altitude (lowest pressure) of merge region (Pa)
61 : real(r8) :: max_pressure_sw=50._r8
62 :
63 : ! lw merge region
64 : ! highest altitude (lowest pressure) of merge region (Pa)
65 : real(r8) :: min_pressure_lw= 5._r8
66 : ! lowest altitude (highest pressure) of merge region (Pa)
67 : real(r8) :: max_pressure_lw=50._r8
68 :
69 : integer :: ntop_qrs_cam ! top level for pure cam solar heating
70 :
71 : !===============================================================================
72 : contains
73 : !===============================================================================
74 :
75 1536 : subroutine radheat_readnl(nlfile)
76 :
77 : use namelist_utils, only: find_group_name
78 : use units, only: getunit, freeunit
79 : use cam_abortutils, only: endrun
80 : use spmd_utils, only : mpicom, masterprocid, mpi_logical, mpi_integer
81 :
82 : use waccm_forcing, only: waccm_forcing_readnl
83 :
84 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
85 :
86 : ! Local variables
87 : integer :: unitn, ierr
88 : character(len=*), parameter :: subname = 'radheat_readnl'
89 :
90 : namelist /radheat_nl/ nlte_use_mo, nlte_limit_co2, nlte_use_aliarms,nlte_aliarms_every_X
91 :
92 1536 : if (masterproc) then
93 2 : unitn = getunit()
94 2 : open( unitn, file=trim(nlfile), status='old' )
95 2 : call find_group_name(unitn, 'radheat_nl', status=ierr)
96 2 : if (ierr == 0) then
97 0 : read(unitn, radheat_nl, iostat=ierr)
98 0 : if (ierr /= 0) then
99 0 : call endrun(subname // ':: ERROR reading namelist')
100 : end if
101 : end if
102 2 : close(unitn)
103 2 : call freeunit(unitn)
104 :
105 : end if
106 :
107 1536 : call mpi_bcast (nlte_use_mo, 1, mpi_logical, masterprocid, mpicom, ierr)
108 1536 : if (ierr /= 0) call endrun("radheat_readnl: FATAL: mpi_bcast: nlte_use_mo")
109 1536 : call mpi_bcast (nlte_limit_co2, 1, mpi_logical, masterprocid, mpicom, ierr)
110 1536 : if (ierr /= 0) call endrun("radheat_readnl: FATAL: mpi_bcast: nlte_limit_co2")
111 1536 : call mpi_bcast (nlte_use_aliarms, 1, mpi_logical, masterprocid, mpicom, ierr)
112 1536 : if (ierr /= 0) call endrun("radheat_readnl: FATAL: mpi_bcast: nlte_use_aliarms")
113 1536 : call mpi_bcast (nlte_aliarms_every_X, 1, mpi_integer, masterprocid, mpicom, ierr)
114 1536 : if (ierr /= 0) call endrun("radheat_readnl: FATAL: mpi_bcast: nlte_aliarms_every_X")
115 :
116 : ! Have waccm_forcing read its namelist as well.
117 1536 : call waccm_forcing_readnl(nlfile)
118 :
119 1536 : end subroutine radheat_readnl
120 :
121 : !================================================================================================
122 :
123 1536 : subroutine radheat_register
124 :
125 1536 : use nlte_lw, only : nlte_register
126 :
127 : ! only ALI-ARMS has pbuf fields to register
128 1536 : if (nlte_use_aliarms) then
129 0 : call nlte_register()
130 : end if
131 :
132 1536 : end subroutine radheat_register
133 :
134 : !================================================================================================
135 :
136 1536 : subroutine radheat_init(pref_mid)
137 :
138 1536 : use nlte_lw, only: nlte_init
139 : use cam_history, only: add_default, addfld
140 : use phys_control, only: phys_getopts
141 : use physics_buffer, only : physics_buffer_desc
142 :
143 : ! args
144 :
145 : real(r8), intent(in) :: pref_mid(pver) ! mid point reference pressure
146 :
147 : ! local vars
148 :
149 : real(r8) :: delta_merge_sw ! range of merge region
150 : real(r8) :: midpoint_sw ! midpoint of merge region
151 : real(r8) :: delta_merge_lw ! range of merge region
152 : real(r8) :: midpoint_lw ! midpoint of merge region
153 : real(r8) :: psh(pver) ! pressure scale height
154 : integer :: k
155 : logical :: camrt
156 :
157 : character(len=16) :: rad_pkg
158 : logical :: history_scwaccm_forcing
159 : logical :: history_waccm
160 :
161 : !-----------------------------------------------------------------------
162 :
163 :
164 : call phys_getopts(radiation_scheme_out=rad_pkg, &
165 : history_waccm_out=history_waccm, &
166 1536 : history_scwaccm_forcing_out=history_scwaccm_forcing)
167 1536 : camrt = rad_pkg == 'CAMRT' .or. rad_pkg == 'camrt'
168 :
169 : ! set max/min pressures for merging regions.
170 :
171 1536 : if (camrt) then
172 0 : min_pressure_sw = 1e5_r8*exp(-10._r8)
173 0 : max_pressure_sw = 1e5_r8*exp(-9._r8)
174 0 : min_pressure_lw = 1e5_r8*exp(-10._r8)
175 0 : max_pressure_lw = 1e5_r8*exp(-8.57_r8)
176 : else
177 1536 : min_pressure_sw = 5._r8
178 1536 : max_pressure_sw = 50._r8
179 1536 : min_pressure_lw = 5._r8
180 1536 : max_pressure_lw = 50._r8
181 : endif
182 :
183 1536 : delta_merge_sw = max_pressure_sw - min_pressure_sw
184 1536 : delta_merge_lw = max_pressure_lw - min_pressure_lw
185 :
186 1536 : midpoint_sw = (max_pressure_sw + min_pressure_sw)/2._r8
187 1536 : midpoint_lw = (max_pressure_lw + min_pressure_lw)/2._r8
188 :
189 109056 : do k=1,pver
190 :
191 : ! pressure scale heights for camrt merging (waccm4)
192 107520 : psh(k)=log(1e5_r8/pref_mid(k))
193 :
194 107520 : if ( pref_mid(k) .le. min_pressure_sw ) then
195 29184 : qrs_wt(k) = 0._r8
196 78336 : else if( pref_mid(k) .ge. max_pressure_sw) then
197 70656 : qrs_wt(k) = 1._r8
198 : else
199 7680 : if (camrt) then
200 : ! camrt
201 0 : qrs_wt(k) = 1._r8 - tanh( (psh(k) - 9._r8)/.75_r8 )
202 : else
203 : ! rrtmg
204 : qrs_wt(k) = 0.5_r8 + 1.5_r8*((pref_mid(k)-midpoint_sw)/delta_merge_sw) &
205 7680 : - 2._r8*((pref_mid(k)-midpoint_sw)/delta_merge_sw)**3._r8
206 : endif
207 : endif
208 :
209 109056 : if ( pref_mid(k) .le. min_pressure_lw ) then
210 29184 : qrl_wt(k)= 0._r8
211 78336 : else if( pref_mid(k) .ge. max_pressure_lw) then
212 70656 : qrl_wt(k)= 1._r8
213 : else
214 7680 : if (camrt) then
215 : ! camrt
216 0 : qrl_wt(k) = 1._r8 - tanh( (psh(k) - 8.57_r8) / 0.71_r8 )
217 : else
218 : ! rrtmg
219 : qrl_wt(k) = 0.5_r8 + 1.5_r8*((pref_mid(k)-midpoint_lw)/delta_merge_lw) &
220 7680 : - 2._r8*((pref_mid(k)-midpoint_lw)/delta_merge_lw)**3._r8
221 : endif
222 : endif
223 :
224 : end do
225 :
226 : ! determine upppermost level that is purely solar heating (no MLT chem heationg)
227 1536 : ntop_qrs_cam = 0
228 109056 : do k=pver,1,-1
229 109056 : if (qrs_wt(k)==1._r8) ntop_qrs_cam = k
230 : enddo
231 :
232 1536 : if (masterproc) then
233 2 : write(iulog,*) 'RADHEAT_INIT: pref_mid', pref_mid(:)
234 2 : write(iulog,*) 'RADHEAT_INIT: QRS_WT ', qrs_wt(:)
235 2 : write(iulog,*) 'RADHEAT_INIT: QRL_WT ', qrl_wt(:)
236 : end if
237 :
238 : ! WACCM heating if top-most layer is above merge region
239 1536 : waccm_heating = (pref_mid(1) .le. min_pressure_sw)
240 :
241 1536 : if (masterproc) then
242 2 : write(iulog,*) 'WACCM Heating is computed (true/false): ',waccm_heating
243 : end if
244 :
245 1536 : if (waccm_heating) then
246 1536 : call nlte_init(pref_mid, max_pressure_lw, nlte_use_mo, nlte_limit_co2, nlte_use_aliarms,nlte_aliarms_every_X)
247 : endif
248 :
249 : ! Add history variables to master field list
250 3072 : call addfld ('QRL_TOT',(/ 'lev' /), 'A','K/s','Merged LW heating: QRL+QRLNLTE')
251 3072 : call addfld ('QRS_TOT',(/ 'lev' /), 'A','K/s','Merged SW heating: QRS+QCP+QRS_EUV+QRS_CO2NIR+QRS_AUR+QTHERMAL')
252 :
253 3072 : call addfld ('QRS_TOT_24_COS',(/ 'lev' /), 'A','K/s','SW heating 24hr. cos coeff.')
254 3072 : call addfld ('QRS_TOT_24_SIN',(/ 'lev' /), 'A','K/s','SW heating 24hr. sin coeff.')
255 3072 : call addfld ('QRS_TOT_12_COS',(/ 'lev' /), 'A','K/s','SW heating 12hr. cos coeff.')
256 3072 : call addfld ('QRS_TOT_12_SIN',(/ 'lev' /), 'A','K/s','SW heating 12hr. sin coeff.')
257 3072 : call addfld ('QRS_TOT_08_COS',(/ 'lev' /), 'A','K/s','SW heating 8hr. cos coeff.')
258 3072 : call addfld ('QRS_TOT_08_SIN',(/ 'lev' /), 'A','K/s','SW heating 8hr. sin coeff.')
259 :
260 : ! Add default history variables to files
261 1536 : if (history_waccm) then
262 1536 : call add_default ('QRL_TOT', 1, ' ')
263 1536 : call add_default ('QRS_TOT', 1, ' ')
264 : end if
265 1536 : if (history_scwaccm_forcing) then
266 0 : call add_default ('QRS_TOT', 8, ' ')
267 : end if
268 :
269 1536 : end subroutine radheat_init
270 :
271 : !================================================================================================
272 :
273 32256 : subroutine radheat_timestep_init (state, pbuf2d)
274 :
275 1536 : use nlte_lw, only: nlte_timestep_init
276 : use physics_types,only : physics_state
277 : use ppgrid, only : begchunk, endchunk
278 : use physics_buffer, only : physics_buffer_desc
279 :
280 : type(physics_state), intent(in):: state(begchunk:endchunk)
281 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
282 :
283 :
284 16128 : if (waccm_heating) then
285 16128 : call nlte_timestep_init (state, pbuf2d)
286 : endif
287 :
288 16128 : end subroutine radheat_timestep_init
289 :
290 : !================================================================================================
291 :
292 17660160 : subroutine radheat_tend(state, pbuf, ptend, qrl, qrs, fsns, &
293 : fsnt, flns, flnt, asdir, net_flx)
294 : !-----------------------------------------------------------------------
295 : ! Compute net radiative heating from qrs and qrl, and the associated net
296 : ! boundary flux.
297 : !
298 : ! This routine provides the waccm hook for computing nonLTE cooling and
299 : ! eUV heating.
300 : !-----------------------------------------------------------------------
301 :
302 16128 : use cam_history, only: outfld
303 : use nlte_lw, only: nlte_tend
304 : use mo_waccm_hrates, only: waccm_hrates, has_hrates
305 : use waccm_forcing, only: get_solar
306 :
307 : use physics_buffer, only : physics_buffer_desc
308 : use tidal_diag, only: get_tidal_coeffs
309 :
310 : ! Arguments
311 : type(physics_state), intent(in) :: state ! Physics state variables
312 :
313 : type(physics_buffer_desc), pointer :: pbuf(:)
314 : type(physics_ptend), intent(out) :: ptend ! indivdual parameterization tendencie
315 : real(r8), intent(in) :: qrl(pcols,pver) ! longwave heating
316 : real(r8), intent(in) :: qrs(pcols,pver) ! shortwave heating
317 : real(r8), intent(in) :: fsns(pcols) ! Surface solar absorbed flux
318 : real(r8), intent(in) :: fsnt(pcols) ! Net column abs solar flux at model top
319 : real(r8), intent(in) :: flns(pcols) ! Srf longwave cooling (up-down) flux
320 : real(r8), intent(in) :: flnt(pcols) ! Net outgoing lw flux at model top
321 : real(r8), intent(in) :: asdir(pcols) ! shortwave, direct albedo
322 : real(r8), intent(out) :: net_flx(pcols)
323 :
324 : ! Local variables
325 : integer :: i, k
326 : integer :: ncol ! number of atmospheric columns
327 : integer :: lchnk ! chunk identifier
328 : real(r8) :: qrl_mrg(pcols,pver) ! merged LW heating
329 : real(r8) :: qrl_mlt(pcols,pver) ! M/LT longwave heating rates
330 : real(r8) :: qrs_mrg(pcols,pver) ! merged SW heating
331 : real(r8) :: qrs_mlt(pcols,pver) ! M/LT solar heating rates
332 : real(r8) :: qout(pcols,pver) ! temp for outfld call
333 : real(r8) :: dcoef(6) ! for tidal component of heating
334 :
335 : !-----------------------------------------------------------------------
336 :
337 80640 : ncol = state%ncol
338 80640 : lchnk = state%lchnk
339 :
340 80640 : call physics_ptend_init(ptend, state%psetcols, 'radheat', ls=.true.)
341 :
342 : ! WACCM interactive heating rate
343 80640 : if (waccm_heating.and.waccm_heating_on) then
344 80640 : call t_startf( 'hrates' )
345 80640 : if (has_hrates) then
346 80640 : call waccm_hrates(ncol, state, asdir, ntop_qrs_cam-1, qrs_mlt, pbuf)
347 : else
348 0 : call get_solar(ncol, lchnk, qrs_mlt)
349 : endif
350 80640 : call t_stopf( 'hrates' )
351 : else
352 0 : qrs_mlt(:,:) = 0._r8
353 : endif
354 :
355 : ! Merge cam solar heating for lower atmosphere with M/LT heating
356 80640 : call merge_qrs (ncol, qrs, qrs_mlt, qrs_mrg, cpairv(:,:,lchnk))
357 87010560 : qout(:ncol,:) = qrs_mrg(:ncol,:)/cpairv(:ncol,:,lchnk)
358 80640 : call outfld ('QRS_TOT', qout, pcols, lchnk)
359 :
360 : ! Output tidal coefficients of total SW heating
361 80640 : call get_tidal_coeffs( dcoef )
362 87010560 : call outfld( 'QRS_TOT_24_SIN', qout(:ncol,:)*dcoef(1), ncol, lchnk )
363 87010560 : call outfld( 'QRS_TOT_24_COS', qout(:ncol,:)*dcoef(2), ncol, lchnk )
364 87010560 : call outfld( 'QRS_TOT_12_SIN', qout(:ncol,:)*dcoef(3), ncol, lchnk )
365 87010560 : call outfld( 'QRS_TOT_12_COS', qout(:ncol,:)*dcoef(4), ncol, lchnk )
366 87010560 : call outfld( 'QRS_TOT_08_SIN', qout(:ncol,:)*dcoef(5), ncol, lchnk )
367 87010560 : call outfld( 'QRS_TOT_08_COS', qout(:ncol,:)*dcoef(6), ncol, lchnk )
368 :
369 80640 : if (waccm_heating.and.waccm_heating_on) then
370 80640 : call t_startf( 'nltedrv' )
371 80640 : call nlte_tend(state, pbuf, qrl_mlt)
372 80640 : call t_stopf( 'nltedrv' )
373 : else
374 0 : qrl_mlt(:,:) = 0._r8
375 : endif
376 :
377 : ! Merge cam long wave heating for lower atmosphere with M/LT (nlte) heating
378 80640 : call merge_qrl (ncol, qrl, qrl_mlt, qrl_mrg)
379 87010560 : qout(:ncol,:) = qrl_mrg(:ncol,:)/cpairv(:ncol,:,lchnk)
380 80640 : call outfld ('QRL_TOT', qout, pcols, lchnk)
381 :
382 87010560 : ptend%s(:ncol,:) = qrs_mrg(:ncol,:) + qrl_mrg(:ncol,:)
383 :
384 80640 : net_flx = 0._r8
385 5725440 : do k = 1, pver
386 87010560 : do i = 1, ncol
387 86929920 : net_flx(i) = net_flx(i) + ptend%s(i,k)*state%pdel(i,k)/gravit
388 : end do
389 : end do
390 :
391 80640 : end subroutine radheat_tend
392 :
393 : !================================================================================================
394 0 : subroutine radheat_disable_waccm()
395 0 : waccm_heating_on = .false.
396 80640 : end subroutine radheat_disable_waccm
397 : !================================================================================================
398 :
399 80640 : subroutine merge_qrs (ncol, hcam, hmlt, hmrg, cpair)
400 : !
401 : ! Merges short wave heating rates
402 : !
403 : implicit none
404 :
405 : !-----------------Input arguments-----------------------------------
406 : integer ncol
407 :
408 : real(r8), intent(in) :: hmlt(pcols,pver) ! Upper atmosphere heating rates
409 : real(r8), intent(in) :: hcam(pcols,pver) ! CAM heating rate
410 : real(r8), intent(out) :: hmrg(pcols,pver) ! merged heating rates
411 : real(r8), intent(in) :: cpair(pcols,pver) ! Specific heat of dry air
412 :
413 : !-----------------Local workspace------------------------------------
414 :
415 : integer k
416 :
417 5725440 : do k = 1, pver
418 87010560 : hmrg(:ncol,k) = qrs_wt(k)*hcam(:ncol,k) + (1._r8 - qrs_wt(k))*cpair(:ncol,k)*hmlt(:ncol,k)
419 : end do
420 :
421 80640 : end subroutine merge_qrs
422 :
423 : !==================================================================================================
424 :
425 80640 : subroutine merge_qrl (ncol, hcam, hmlt, hmrg)
426 : !
427 : ! Merges long wave heating rates
428 : !
429 : !-----------------Input arguments-----------------------------------
430 : integer ncol
431 :
432 : real(r8), intent(in) :: hmlt(pcols,pver) ! Upper atmosphere heating rates
433 : real(r8), intent(in) :: hcam(pcols,pver) ! CAM heating rate
434 : real(r8), intent(out) :: hmrg(pcols,pver) ! merged heating rates
435 :
436 : !-----------------Local workspace------------------------------------
437 :
438 : integer k
439 :
440 : !--------------------------------------------------------------------
441 :
442 5725440 : do k = 1, pver
443 87010560 : hmrg(:ncol,k) = qrl_wt(k) * hcam(:ncol,k) + (1._r8-qrl_wt(k)) * hmlt(:ncol,k)
444 : end do
445 :
446 80640 : end subroutine merge_qrl
447 :
448 : end module radheat
|
