LCOV - code coverage report
Current view: | top level - chemistry/mozart - mo_waccm_hrates.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 15 | 181 | 8.3 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 1 | 2 | 50.0 %

              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

Generated by: LCOV version 1.14