Line data Source code
1 :
2 : module mo_waccm_hrates
3 :
4 : use shr_kind_mod, only : r8 => shr_kind_r8
5 : use cam_logfile, only : iulog
6 : use physics_buffer, only : pbuf_get_index, pbuf_get_field
7 :
8 : implicit none
9 :
10 : save
11 :
12 : real(r8), parameter :: secpday = 86400._r8
13 : real(r8), parameter :: daypsec = 1._r8/secpday
14 : real(r8), parameter :: aur_therm = 807._r8
15 : real(r8), parameter :: jkcal = 4184._r8
16 : real(r8), parameter :: aur_heat_eff = .05_r8
17 : real(r8), parameter :: aur_hconst = 1.e3_r8*jkcal*aur_therm*aur_heat_eff
18 :
19 : real(r8) :: max_zen_angle
20 :
21 : private
22 : public :: waccm_hrates, init_hrates, has_hrates
23 :
24 : integer :: id_co2, id_o2, id_o3, id_o2_1d, id_o2_1s, id_o1d, id_h2o, id_o, id_h
25 : logical :: has_hrates
26 : integer :: ele_temp_ndx, ion_temp_ndx
27 :
28 : contains
29 :
30 1536 : subroutine init_hrates( )
31 : use mo_chem_utls, only : get_spc_ndx
32 : use cam_history, only : addfld
33 : use ref_pres, only : ptop_ref, psurf_ref
34 :
35 :
36 : implicit none
37 :
38 : integer :: ids(9), err
39 : character(len=128) :: attr ! netcdf variable attribute
40 :
41 1536 : id_co2 = get_spc_ndx( 'CO2' )
42 1536 : id_o2 = get_spc_ndx( 'O2' )
43 1536 : id_o3 = get_spc_ndx( 'O3' )
44 1536 : id_o2_1d = get_spc_ndx( 'O2_1D' )
45 1536 : id_o2_1s = get_spc_ndx( 'O2_1S' )
46 1536 : id_o1d = get_spc_ndx( 'O1D' )
47 1536 : id_h2o = get_spc_ndx( 'H2O' )
48 1536 : id_o = get_spc_ndx( 'O' )
49 1536 : id_h = get_spc_ndx( 'H' )
50 :
51 15360 : ids = (/ id_co2, id_o2, id_o3, id_o2_1d, id_o2_1s, id_o1d, id_h2o, id_o, id_h /)
52 :
53 3072 : has_hrates = all( ids(:) > 0 ) .and. ptop_ref < 0.0004_r8 * psurf_ref
54 :
55 1536 : if (.not. has_hrates) return
56 :
57 0 : call addfld( 'CPAIR', (/ 'lev' /), 'I', 'J/K/kg', 'specific heat cap air' )
58 0 : call addfld( 'QRS_AUR', (/ 'lev' /), 'I', 'K/s', 'total auroral heating rate' )
59 0 : call addfld( 'QRS_CO2NIR', (/ 'lev' /), 'I', 'K/s', 'co2 nir heating rate' )
60 0 : call addfld( 'QTHERMAL', (/ 'lev' /), 'I', 'K/s', 'non-euv photolysis heating rate' )
61 0 : call addfld( 'QRS_MLT', (/ 'lev' /), 'I', 'K/s', 'Total heating rate (unmerged with tropospheric RT heating)' )
62 :
63 0 : attr = 'O2 + hv -> O1D + O3P solar heating rate < 200nm'
64 0 : call addfld( 'QRS_SO2A', (/ 'lev' /), 'I', 'K/s', trim(attr) )
65 0 : attr = 'O2 + hv -> O3P + O3P solar heating rate < 200nm'
66 0 : call addfld( 'QRS_SO2B', (/ 'lev' /), 'I', 'K/s', trim(attr) )
67 0 : attr = 'O3 + hv -> O1D + O2_1S solar heating rate < 200nm'
68 0 : call addfld( 'QRS_SO3A', (/ 'lev' /), 'I', 'K/s', trim(attr) )
69 0 : attr = 'O3 + hv -> O3P + O2 solar heating rate < 200nm'
70 0 : call addfld( 'QRS_SO3B', (/ 'lev' /), 'I', 'K/s', trim(attr) )
71 0 : attr = 'O2 + hv -> 2*O3P solar heating rate > 200nm'
72 0 : call addfld( 'QRS_LO2B', (/ 'lev' /), 'I', 'K/s', trim(attr) )
73 0 : attr = 'O3 + hv -> O1D + O2_1S solar heating rate > 200nm'
74 0 : call addfld( 'QRS_LO3A', (/ 'lev' /), 'I', 'K/s', trim(attr) )
75 0 : attr = 'O3 + hv -> O3P + O2 solar heating rate > 200nm'
76 0 : call addfld( 'QRS_LO3B', (/ 'lev' /), 'I', 'K/s', trim(attr) )
77 0 : attr = 'Total O3 solar heating > 200nm'
78 0 : call addfld( 'QRS_LO3', (/ 'lev' /), 'I', 'K/s', trim(attr) )
79 0 : attr = 'total euv heating rate'
80 0 : call addfld( 'QRS_EUV', (/ 'lev' /), 'I', 'K/s', trim(attr) )
81 0 : attr = 'total jo2 euv photolysis rate'
82 0 : call addfld( 'JO2_EUV', (/ 'lev' /), 'I', '/s', trim(attr) )
83 :
84 0 : ele_temp_ndx = pbuf_get_index('TElec',errcode=err)! electron temperature index
85 0 : ion_temp_ndx = pbuf_get_index('TIon',errcode=err) ! ion temperature index
86 :
87 1536 : end subroutine init_hrates
88 :
89 0 : subroutine waccm_hrates(ncol, state, asdir, bot_mlt_lev, qrs_tot, pbuf )
90 : !-----------------------------------------------------------------------
91 : ! ... computes the short wavelength heating rates
92 : !-----------------------------------------------------------------------
93 :
94 1536 : use chem_mods, only : nabscol, nfs, gas_pcnst, rxntot, indexm
95 : use ppgrid, only : pcols, pver
96 : use physconst, only : rga
97 : use air_composition, only : mbarv, cpairv
98 : use constituents, only : pcnst
99 : use mo_gas_phase_chemdr,only: map2chm
100 : use mo_photo, only : set_ub_col, setcol
101 : use mo_jlong, only : jlong
102 : use mo_jshort, only : jshort
103 : use mo_jeuv, only : heuv
104 : use mo_cph, only : cph
105 : use mo_heatnirco2, only : heatnirco2
106 : use mo_airglow, only : airglow
107 : use mo_aurora, only : aurora
108 : use mo_setrxt, only : setrxt_hrates
109 : use mo_adjrxt, only : adjrxt
110 : use mo_usrrxt, only : usrrxt_hrates
111 : use mo_setinv, only : setinv
112 : use mo_mass_xforms, only : mmr2vmr
113 : use physics_types, only : physics_state
114 : use phys_grid, only : get_rlat_all_p, get_rlon_all_p
115 : use mo_mean_mass, only : set_mean_mass
116 : use set_cp, only : calc_cp
117 : use cam_history, only : outfld
118 : use shr_orb_mod, only : shr_orb_decl
119 : use time_manager, only : get_curr_calday
120 : use cam_control_mod, only : lambm0, eccen, mvelpp, obliqr
121 : use mo_constants, only : r2d, n2min
122 : use short_lived_species,only: get_short_lived_species
123 : use physics_buffer, only : physics_buffer_desc
124 : use phys_control, only : waccmx_is
125 : use orbit, only : zenith
126 :
127 : !-----------------------------------------------------------------------
128 : ! ... dummy arguments
129 : !-----------------------------------------------------------------------
130 : integer, intent(in) :: ncol ! number columns in chunk
131 : type(physics_state),target, intent(in) :: state ! physics state structure
132 : real(r8), intent(in) :: asdir(pcols) ! shortwave, direct albedo
133 : integer, intent(in) :: bot_mlt_lev ! lowest model level where MLT heating is needed
134 : real(r8), intent(out) :: qrs_tot(pcols,pver) ! total heating (K/s)
135 : type(physics_buffer_desc), pointer :: pbuf(:)
136 :
137 : !-----------------------------------------------------------------------
138 : ! ... local variables
139 : !-----------------------------------------------------------------------
140 : integer :: lchnk ! chunk index
141 : real(r8), parameter :: m2km = 1.e-3_r8
142 : real(r8), parameter :: Pa2mb = 1.e-2_r8
143 :
144 : integer :: i, k, m, n
145 : integer :: kbot_hrates
146 : real(r8) :: esfact
147 : real(r8) :: sza ! solar zenith angle (degrees)
148 0 : real(r8) :: invariants(ncol,pver,nfs)
149 0 : real(r8) :: col_dens(ncol,pver,max(1,nabscol)) ! column densities (molecules/cm^2)
150 0 : real(r8) :: col_delta(ncol,0:pver,max(1,nabscol)) ! layer column densities (molecules/cm^2)
151 0 : real(r8) :: vmr(ncol,pver,gas_pcnst) ! xported species (vmr)
152 0 : real(r8) :: reaction_rates(ncol,pver,rxntot) ! reaction rates
153 : real(r8) :: mmr(pcols,pver,gas_pcnst) ! chem working concentrations (kg/kg)
154 0 : real(r8) :: h2ovmr(ncol,pver) ! water vapor concentration (mol/mol)
155 0 : real(r8) :: mbar(ncol,pver) ! mean wet atmospheric mass (kg/mole)
156 0 : real(r8) :: zmid(ncol,pver) ! midpoint geopotential (km)
157 0 : real(r8) :: cpair(ncol,pver) ! specific heat capacity (J/K/kg)
158 0 : real(r8) :: cphrate(ncol,pver) ! chemical pot heat rate (K/s)
159 0 : real(r8) :: aghrate(ncol,pver) ! airglow heat rate (K/s)
160 : real(r8) :: qrs_col(pver,4) ! column thermal heating < 200nm
161 : real(r8) :: qrl_col(pver,4) ! column thermal heating > 200nm
162 0 : real(r8) :: qrs(ncol,pver,4) ! chunk thermal heating < 200nm
163 0 : real(r8) :: qrl(ncol,pver,4) ! chunk thermal heating > 200nm
164 : real(r8) :: euv_hrate_col(pver) ! column euv thermal heating rate
165 : real(r8) :: co2_hrate_col(pver) ! column co2 nir heating rate
166 0 : real(r8) :: euv_hrate(ncol,pver) ! chunk euv thermal heating rate
167 0 : real(r8) :: aur_hrate(ncol,pver) ! chunk auroral heating rate
168 0 : real(r8) :: co2_hrate(ncol,pver) ! chunk co2 nir heating rate
169 : real(r8) :: colo3(pver) ! vertical o3 column density
170 : real(r8) :: zarg(pver) ! vertical height array
171 : real(r8) :: parg(pver) ! vertical pressure array (hPa)
172 : real(r8) :: tline(pver) ! vertical temperature array
173 : real(r8) :: eff_alb(pver) ! albedo
174 : real(r8) :: mw(pver) ! atms molecular weight
175 : real(r8) :: n2_line(pver) ! n2 density (mol/mol)
176 : real(r8) :: o_line(pver) ! o density (mol/mol)
177 : real(r8) :: o2_line(pver) ! o2 density (mol/mol)
178 : real(r8) :: o3_line(pver) ! o3 density (mol/mol)
179 : real(r8) :: co2_line(pver) ! co2 density (mol/mol)
180 : real(r8) :: scco2(pver) ! co2 slant column concentration (molec/cm^2)
181 : real(r8) :: scco2i(pver) ! co2 slant column concentration (molec/cm^2)
182 : real(r8) :: occ(pver) ! o density (molecules/cm^3)
183 : real(r8) :: o2cc(pver) ! o2 density (molecules/cm^3)
184 : real(r8) :: co2cc(pver) ! co2 density (molecules/cm^3)
185 : real(r8) :: n2cc(pver) ! n2 density (molecules/cm^3)
186 : real(r8) :: o3cc(pver) ! o3 density (molecules/cm^3)
187 : real(r8) :: cparg(pver) ! specific heat capacity
188 0 : real(r8) :: zen_angle(ncol) ! solar zenith angles (radians)
189 0 : real(r8) :: zsurf(ncol) ! surface height (m)
190 0 : real(r8) :: rlats(ncol) ! chunk latitudes (radians)
191 0 : real(r8) :: rlons(ncol) ! chunk longitudes (radians)
192 : real(r8) :: calday ! day of year
193 : real(r8) :: delta ! solar declination (radians)
194 : logical :: do_diag
195 :
196 0 : real(r8), pointer :: ele_temp_fld(:,:) ! electron temperature pointer
197 0 : real(r8), pointer :: ion_temp_fld(:,:) ! ion temperature pointer
198 :
199 0 : if ( ele_temp_ndx>0 .and. ion_temp_ndx>0 ) then
200 0 : call pbuf_get_field(pbuf, ele_temp_ndx, ele_temp_fld)
201 0 : call pbuf_get_field(pbuf, ion_temp_ndx, ion_temp_fld)
202 : else
203 0 : ele_temp_fld => state%t
204 0 : ion_temp_fld => state%t
205 : endif
206 :
207 0 : qrs_tot(:ncol,:) = 0._r8
208 0 : if (.not. has_hrates) return
209 :
210 : !-------------------------------------------------------------------------
211 : ! ... set maximum zenith angle - higher value for higher top model
212 : !-------------------------------------------------------------------------
213 0 : if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then
214 0 : max_zen_angle = 116._r8
215 : else
216 0 : max_zen_angle = 97.01_r8 ! degrees
217 : endif
218 :
219 : !-----------------------------------------------------------------------
220 : ! ... get chunk latitudes and longitudes
221 : !-----------------------------------------------------------------------
222 0 : lchnk = state%lchnk
223 :
224 0 : call get_rlat_all_p( lchnk, ncol, rlats )
225 0 : call get_rlon_all_p( lchnk, ncol, rlons )
226 :
227 : !-----------------------------------------------------------------------
228 : ! ... set lower limit for heating rates which is now dictated by radheat module
229 : !-----------------------------------------------------------------------
230 0 : kbot_hrates = bot_mlt_lev
231 0 : kbot_hrates = min( kbot_hrates,pver )
232 : ! write(iulog,*) 'hrates: kbot_hrates = ',kbot_hrates
233 :
234 : !-----------------------------------------------------------------------
235 : ! ... calculate cosine of zenith angle then cast back to angle
236 : !-----------------------------------------------------------------------
237 0 : calday = get_curr_calday()
238 0 : call zenith( calday, rlats, rlons, zen_angle, ncol )
239 0 : zen_angle(:) = acos( zen_angle(:) )
240 :
241 : !-----------------------------------------------------------------------
242 : ! ... map incoming concentrations to working array
243 : !-----------------------------------------------------------------------
244 0 : do m = 1,pcnst
245 0 : n = map2chm(m)
246 0 : if( n > 0 ) then
247 0 : do k = 1,pver
248 0 : mmr(:ncol,k,n) = state%q(:ncol,k,m)
249 : end do
250 : end if
251 : end do
252 0 : call get_short_lived_species( mmr, lchnk, ncol, pbuf )
253 :
254 : !-----------------------------------------------------------------------
255 : ! ... set atmosphere mean mass
256 : !-----------------------------------------------------------------------
257 0 : if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then
258 0 : do k = 1,pver
259 0 : mbar(:ncol,k) = mbarv(:ncol,k,lchnk)
260 : enddo
261 : else
262 0 : call set_mean_mass( ncol, mmr, mbar )
263 : endif
264 : !
265 : !-----------------------------------------------------------------------
266 : ! ... xform from mmr to vmr
267 : !-----------------------------------------------------------------------
268 0 : call mmr2vmr( mmr(:ncol,:,:), vmr(:ncol,:,:), mbar(:ncol,:), ncol )
269 : !-----------------------------------------------------------------------
270 : ! ... xform water vapor from mmr to vmr
271 : !-----------------------------------------------------------------------
272 0 : do k = 1,pver
273 0 : h2ovmr(:ncol,k) = vmr(:ncol,k,id_h2o)
274 : end do
275 : !-----------------------------------------------------------------------
276 : ! ... xform geopotential height from m to km
277 : ! and pressure from Pa to mb
278 : !-----------------------------------------------------------------------
279 0 : zsurf(:ncol) = rga * state%phis(:ncol)
280 0 : do k = 1,pver
281 0 : zmid(:ncol,k) = m2km * (state%zm(:ncol,k) + zsurf(:ncol))
282 : end do
283 :
284 : !-----------------------------------------------------------------------
285 : ! ... set the "invariants"
286 : !-----------------------------------------------------------------------
287 0 : call setinv( invariants, state%t, h2ovmr, vmr, state%pmid, ncol, lchnk, pbuf )
288 :
289 : !-----------------------------------------------------------------------
290 : ! ... set the column densities at the upper boundary
291 : !-----------------------------------------------------------------------
292 0 : call set_ub_col( col_delta, vmr, invariants, state%pint(:,1), state%pdel, ncol, lchnk )
293 :
294 : !-----------------------------------------------------------------------
295 : ! ... set rates for "tabular" and user specified reactions
296 : !-----------------------------------------------------------------------
297 0 : do m = 1,rxntot
298 0 : do k = 1,pver
299 0 : reaction_rates(:,k,m) = 0._r8
300 : end do
301 : end do
302 0 : call setrxt_hrates( reaction_rates, state%t, invariants(1,1,indexm), ncol, kbot_hrates )
303 : call usrrxt_hrates( reaction_rates, state%t, ele_temp_fld, ion_temp_fld, &
304 0 : h2ovmr, invariants(:,:,indexm), ncol, kbot_hrates )
305 0 : call adjrxt( reaction_rates, invariants, invariants(1,1,indexm), ncol,pver )
306 :
307 : !-----------------------------------------------------------------------
308 : ! ... set cp array
309 : !-----------------------------------------------------------------------
310 0 : if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then
311 0 : do k = 1, pver
312 0 : cpair(:ncol,k) = cpairv(:ncol,k,lchnk)
313 : enddo
314 : else
315 0 : call calc_cp( ncol, vmr, cpair )
316 : endif
317 :
318 0 : call outfld( 'CPAIR', cpair, ncol, lchnk )
319 : #ifdef HRATES_DEBUG
320 : write(iulog,*) ' '
321 : write(iulog,*) '---------------------------------------------'
322 : write(iulog,*) 'waccm_hrates: cp at lchnk = ',lchnk
323 : write(iulog,'(1p,5g15.7)') cpair(1,:)
324 : write(iulog,*) '---------------------------------------------'
325 : write(iulog,*) ' '
326 : #endif
327 :
328 : !-----------------------------------------------------------------------
329 : ! ... set the earth-sun distance factor
330 : !-----------------------------------------------------------------------
331 : call shr_orb_decl( calday, eccen, mvelpp, lambm0, obliqr , &
332 0 : delta, esfact )
333 : !-----------------------------------------------------------------------
334 : ! ... set the column densities
335 : !-----------------------------------------------------------------------
336 0 : call setcol( col_delta, col_dens )
337 : !-----------------------------------------------------------------------
338 : ! ... compute the thermal heating rates
339 : !-----------------------------------------------------------------------
340 0 : do m = 1,4
341 0 : do k = 1,pver
342 0 : qrs(:,k,m) = 0._r8
343 0 : qrl(:,k,m) = 0._r8
344 : end do
345 : end do
346 0 : do k = 1,pver
347 0 : euv_hrate(:,k) = 0._r8
348 0 : co2_hrate(:,k) = 0._r8
349 : end do
350 : column_loop : &
351 0 : do i = 1,ncol
352 0 : sza = zen_angle(i)*r2d
353 0 : if( sza < max_zen_angle ) then
354 0 : zarg(:) = zmid(i,:)
355 0 : parg(:) = Pa2mb*state%pmid(i,:)
356 0 : colo3(:) = col_dens(i,:,1)
357 0 : tline(:) = state%t(i,:)
358 0 : eff_alb(:) = asdir(i)
359 0 : o_line(:) = vmr(i,:,id_o)
360 0 : o2_line(:) = vmr(i,:,id_o2)
361 0 : co2_line(:) = vmr(i,:,id_co2)
362 0 : n2_line(:) = 1._r8 - (o_line(:) + o2_line(:) + vmr(i,:,id_h))
363 0 : where( n2_line(:) < n2min )
364 : n2_line = n2min
365 : end where
366 0 : o3_line(:) = vmr(i,:,id_o3)
367 0 : occ(:) = o_line(:) * invariants(i,:,indexm)
368 0 : o2cc(:) = o2_line(:) * invariants(i,:,indexm)
369 0 : co2cc(:) = co2_line(:) * invariants(i,:,indexm)
370 0 : n2cc(:) = n2_line(:) * invariants(i,:,indexm)
371 0 : o3cc(:) = o3_line(:) * invariants(i,:,indexm)
372 0 : mw(:) = mbar(i,:)
373 0 : cparg(:) = cpair(i,:)
374 0 : do_diag = .false.
375 : call jshort( pver, sza, o2_line, o3_line, o2cc, &
376 : o3cc, tline, zarg, mw, qrs_col, &
377 0 : cparg, lchnk, i, co2cc, scco2, do_diag )
378 : call jlong( pver, sza, eff_alb, parg, tline, &
379 : mw, o2_line, o3_line, colo3, qrl_col, &
380 0 : cparg, kbot_hrates )
381 0 : do m = 1,4
382 0 : qrs(i,pver:1:-1,m) = qrs_col(:,m) * esfact
383 : end do
384 0 : do m = 2,4
385 0 : qrl(i,:,m) = qrl_col(:,m) * esfact
386 : end do
387 : call heuv( pver, sza, occ, o2cc, n2cc, &
388 : o_line, o2_line, n2_line, cparg, mw, &
389 0 : zarg, euv_hrate_col, kbot_hrates )
390 0 : euv_hrate(i,:) = euv_hrate_col(:) * esfact
391 0 : scco2i(1:pver) = scco2(pver:1:-1)
392 0 : call heatnirco2( co2_line, scco2i, state%pmid(i,:kbot_hrates), co2_hrate_col, kbot_hrates, &
393 0 : zarg, sza )
394 : #ifdef HRATES_DEBUG
395 : write(iulog,*) '==================================='
396 : write(iulog,*) 'hrates: diagnostics for heatco2nir'
397 : write(iulog,*) 'hrates: co2_line'
398 : write(iulog,'(1p,5g15.7)') co2_line(:)
399 : write(iulog,*) 'hrates: scco2'
400 : write(iulog,'(1p,5g15.7)') scco2i(:)
401 : write(iulog,*) 'hrates: co2_hrate'
402 : write(iulog,'(1p,5g15.7)') co2_hrate_col(:)
403 : write(iulog,*) '==================================='
404 : #endif
405 0 : co2_hrate(i,:kbot_hrates) = co2_hrate_col(:kbot_hrates) * esfact * daypsec
406 : end if
407 : end do column_loop
408 :
409 :
410 0 : call outfld( 'QRS_SO2A', qrs(:,:,1), ncol, lchnk )
411 0 : call outfld( 'QRS_SO2B', qrs(:,:,2), ncol, lchnk )
412 0 : call outfld( 'QRS_SO3A', qrs(:,:,3), ncol, lchnk )
413 0 : call outfld( 'QRS_SO3B', qrs(:,:,4), ncol, lchnk )
414 0 : call outfld( 'QRS_LO2B', qrl(:,:,2), ncol, lchnk )
415 0 : call outfld( 'QRS_LO3A', qrl(:,:,3), ncol, lchnk )
416 0 : call outfld( 'QRS_LO3B', qrl(:,:,4), ncol, lchnk )
417 0 : call outfld( 'QRS_LO3', qrl(:,:,3)+qrl(:,:,4), ncol, lchnk )
418 0 : call outfld( 'QRS_EUV', euv_hrate(:,:), ncol, lchnk )
419 0 : call outfld( 'QRS_CO2NIR', co2_hrate(:,:), ncol, lchnk )
420 :
421 : !-----------------------------------------------------------------------
422 : ! ... chemical pot heating rate
423 : !-----------------------------------------------------------------------
424 : call cph( cphrate, vmr, reaction_rates, cpair, mbar, &
425 0 : kbot_hrates, ncol, lchnk )
426 :
427 : !-----------------------------------------------------------------------
428 : ! ... auroral ion production
429 : !-----------------------------------------------------------------------
430 : call aurora( state%t, mbar, rlats, &
431 : aur_hrate, cpair, state%pmid, lchnk, calday, &
432 0 : ncol, rlons, pbuf )
433 0 : do k = 1,pver
434 0 : aur_hrate(:,k) = aur_hrate(:,k)/invariants(:,k,indexm)
435 : end do
436 0 : call outfld( 'QRS_AUR', aur_hrate(:,:), ncol, lchnk )
437 :
438 : !-----------------------------------------------------------------------
439 : ! ... airglow heating rate
440 : !-----------------------------------------------------------------------
441 0 : call airglow( aghrate, vmr(1,1,id_o2_1s), vmr(1,1,id_o2_1d), vmr(1,1,id_o1d), reaction_rates, cpair, &
442 0 : ncol, lchnk )
443 :
444 : !-----------------------------------------------------------------------
445 : ! ... form total heating rate
446 : !-----------------------------------------------------------------------
447 0 : do k = 1,kbot_hrates
448 0 : qrs_tot(:ncol,k) = qrs(:,k,1) + qrs(:,k,2) + qrs(:,k,3) + qrs(:,k,4) &
449 0 : + qrl(:,k,1) + qrl(:,k,2) + qrl(:,k,3) + qrl(:,k,4)
450 : end do
451 0 : call outfld( 'QTHERMAL', qrs_tot, pcols, lchnk )
452 0 : do k = 1,kbot_hrates
453 0 : qrs_tot(:ncol,k) = qrs_tot(:ncol,k) &
454 0 : + cphrate(:,k) + euv_hrate(:,k) + aur_hrate(:,k) + co2_hrate(:,k)
455 : end do
456 0 : call outfld( 'QRS_MLT', qrs_tot, pcols, lchnk )
457 :
458 0 : end subroutine waccm_hrates
459 :
460 : end module mo_waccm_hrates
|
