LCOV - code coverage report
Current view: | top level - chemistry/mozart - mo_jeuv.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 10 | 267 | 3.7 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 1 | 3 | 33.3 %

              Line data    Source code
           1             :        module mo_jeuv
           2             :
           3             :        use shr_kind_mod,     only : r8 => shr_kind_r8
           4             :        use cam_abortutils,   only : endrun
           5             :        use cam_logfile,      only : iulog
           6             :        use euvac,            only : euvac_etf
           7             :        use solar_euv_data,   only : solar_euv_data_etf, solar_euv_data_active
           8             :        use spmd_utils,       only : masterproc, masterprocid, mpicom, mpi_real8
           9             :
          10             :        implicit none
          11             :
          12             :        private
          13             :        public :: jeuv_init, jeuv, heuv, neuv
          14             :        public :: nIonRates
          15             :
          16             :        save
          17             :
          18             : !------------------------------------------------------------------------------
          19             : !      define EUV photolysis cross sections, branching ratios,
          20             : !      wavelength parameters,etc
          21             : !------------------------------------------------------------------------------
          22             :        integer, parameter :: neuv = 26               ! number of photolysis/ionization reactions
          23             :        integer, parameter :: nIonRates = 11          ! number of photo-ionizations rates needed for waccmx
          24             :        integer, parameter :: nmaj = 3                ! number of major neutral species (O,O2,N2)
          25             :        integer, parameter :: nspecies = 5            ! number of neutral species(O,N,O2,N2,CO2)
          26             :        integer, parameter :: nstat = 6               ! maximum number of ionization/dissociation
          27             :                                                      ! /excitation states for each speies
          28             :        integer, parameter :: lmax = 23               ! number of wavelength bins in EUV
          29             :
          30             :        real(r8), parameter :: heat_eff_fac = .08_r8  ! heating efficiency factor
          31             :        ! Solar EUV direct heating was increased from 5% to 8% to bring it closer to the
          32             :        ! TIE-GCM value (5% applied twice for a total of 10%) -- Hanli Liu & Stan Solomon
          33             :
          34             :        real(r8), parameter :: hc = 6.62608e-34_r8 * 2.9979e8_r8 / 1.e-9_r8
          35             :
          36             :        real(r8) :: sigabs(lmax,nspecies)               ! absorption cross sections of major species
          37             :        real(r8) :: branch_p(lmax,nstat,nmaj) = 0._r8   ! branching ratios for photoionization/dissociation
          38             :        real(r8) :: branch_e(lmax,nstat,nmaj) = 0._r8   ! branching ratios for photoelectron ionization/dissociation/excitation
          39             :        real(r8) :: energy(lmax)                        ! solar energy
          40             :
          41             :        real(r8), pointer :: solar_etf(:)
          42             :        logical :: do_heating(13)
          43             :
          44             :        contains
          45             :
          46        1536 :        subroutine jeuv_init (photon_file, electron_file, indexer)
          47             : !==============================================================================
          48             : !   Purpose:
          49             : !      read tabulated data:
          50             : !           (1) thermosphere neutral species' absorption cross sections,
          51             : !               photoionization/dissociation branching ratios
          52             : !           (2) read photoelectron enhancement factor, photoelectron ionization/
          53             : !               dissociation/excitation branching ratios
          54             : !           (3) read solar flux
          55             : !==============================================================================
          56             :
          57             :        use units,         only : getunit, freeunit
          58             :        use ioFileMod,     only : getfil
          59             :        use mo_chem_utls,  only : get_rxt_ndx
          60             :
          61             :        implicit none
          62             :
          63             :        character(len=*), intent(in) :: photon_file
          64             :        character(len=*), intent(in) :: electron_file
          65             :        integer, optional,intent(inout) :: indexer(:)
          66             :
          67             : !------------------------------------------------------------------------------
          68             : !       local variables
          69             : !------------------------------------------------------------------------------
          70             :         integer  :: i, j, m                      ! loop indicies
          71             :         integer  :: unit                         ! fortran i/o unit number
          72             :         integer  :: istat                        ! file op status
          73             :         real(r8) :: waves(lmax)                  ! wavelength array for short bound of bins (A)
          74             :         real(r8) :: wavel(lmax)                  ! wavelength array for long bound of bins (A)
          75             :         real(r8) :: wc(lmax)                     ! wavelength bin center (nm)
          76             :         character(len=200) :: str,fmt            ! string for comments in data file
          77             :         character(len=256) :: locfn
          78             :
          79             :         integer :: jeuv_1_ndx, ierr, jndx
          80             :         character(len=2) :: mstring
          81             :         character(len=7) :: jstring
          82             :         logical :: do_jeuv
          83             :
          84        1536 :         do_jeuv=.false.
          85        1536 :         do_heating=.false.
          86             :
          87       41472 :         do m = 1,neuv
          88       39936 :            write ( mstring, '(I2)' ) m
          89       39936 :            jstring = 'jeuv_'//trim(adjustl(mstring))
          90       39936 :            jndx = get_rxt_ndx( jstring )
          91       39936 :            if (jndx>0) then
          92           0 :              do_jeuv=.true.
          93           0 :              indexer(m) = jndx
          94             :            endif
          95       41472 :            if (jndx>0 .and. m<14) then
          96           0 :              do_heating(m) = .true.
          97             :            endif
          98             :         enddo
          99             :
         100        1536 :         if (.not.do_jeuv) return
         101             :
         102           0 :         if (solar_euv_data_active) then
         103           0 :            solar_etf => solar_euv_data_etf
         104             :         else
         105           0 :            solar_etf => euvac_etf
         106             :         endif
         107             :
         108           0 :         if ( size(solar_etf) /= lmax ) then
         109           0 :            write(iulog,*) 'jeuv_init: the size of solar_etf is incorrect '
         110           0 :            write(iulog,*) ' ... size(solar_etf) = ',size(solar_etf)
         111           0 :            write(iulog,*) ' .............. lmax = ',lmax
         112           0 :            call endrun('jeuv_init: the size of solar_etf is incorrect ')
         113             :         endif
         114             :
         115             : !------------------------------------------------------------------------------
         116             : !       read from data file the absorption cross sections for neutral species,
         117             : !       braching ratios for photoionization/dissociation, and braching ratios
         118             : !       for photoelectron ionization/dissociation/excitation
         119             : !------------------------------------------------------------------------------
         120           0 :         master: if (masterproc) then ! read in ascii data only on master task then b-cast
         121             :            !------------------------------------------------------------------------------
         122             :            ! read neutral species' absorption cross section and
         123             :            ! photoionization/dissociation branching ratio
         124             :            !------------------------------------------------------------------------------
         125           0 :            unit = getunit()
         126           0 :            call getfil( photon_file, locfn, 0 )
         127           0 :            open( unit, file = trim(locfn), status='UNKNOWN', iostat=istat )
         128           0 :            if( istat /= 0 ) then
         129           0 :               write(iulog,*) 'jeuv_init: failed to open ',trim(locfn),'; error = ',istat
         130           0 :               call endrun
         131             :            end if
         132             :            !------------------------------------------------------------------------------
         133             :            ! read O
         134             :            !------------------------------------------------------------------------------
         135           0 :            read(unit,*,iostat=istat) str
         136           0 :            if( istat /= 0 ) then
         137           0 :               write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         138           0 :               call endrun
         139             :            end if
         140           0 :            read(unit,*,iostat=istat) str
         141           0 :            if( istat /= 0 ) then
         142           0 :               write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         143           0 :               call endrun
         144             :            end if
         145           0 :            fmt = "(f7.2,2x,f7.2,2x,f9.3,4(2x,f7.3))"
         146           0 :            do i = 1,lmax
         147           0 :               read(unit,fmt,iostat=istat) waves(i), wavel(i), sigabs(i,1), (branch_p(i,j,1),j=1,4)
         148           0 :               if( istat /= 0 ) then
         149           0 :                  write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         150           0 :                  call endrun
         151             :               end if
         152             :            end do
         153             :            !------------------------------------------------------------------------------
         154             :            ! read O2
         155             :            !------------------------------------------------------------------------------
         156           0 :            read(unit,*,iostat=istat) str
         157           0 :            if( istat /= 0 ) then
         158           0 :               write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         159           0 :               call endrun
         160             :            end if
         161           0 :            read(unit,*,iostat=istat) str
         162           0 :            if( istat /= 0 ) then
         163           0 :               write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         164           0 :               call endrun
         165             :            end if
         166           0 :            fmt = "(f7.2,2x,f7.2,2x,f9.3,5(2x,f7.3))"
         167           0 :            do i = 1,lmax
         168           0 :               read(unit,fmt,iostat=istat) waves(i), wavel(i), sigabs(i,2), (branch_p(i,j,2),j=1,5)
         169           0 :               if( istat /= 0 ) then
         170           0 :                  write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         171           0 :                  call endrun
         172             :               end if
         173             :            end do
         174             :            !------------------------------------------------------------------------------
         175             :            ! read N2
         176             :            !------------------------------------------------------------------------------
         177           0 :            read(unit,*,iostat=istat) str
         178           0 :            if( istat /= 0 ) then
         179           0 :               write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         180           0 :               call endrun
         181             :            end if
         182           0 :            read(unit,*,iostat=istat) str
         183           0 :            if( istat /= 0 ) then
         184           0 :               write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         185           0 :               call endrun
         186             :            end if
         187           0 :            do i = 1,lmax
         188           0 :               read(unit,fmt,iostat=istat) waves(i), wavel(i), sigabs(i,3), (branch_p(i,j,3),j=1,5)
         189           0 :               if( istat /= 0 ) then
         190           0 :                  write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         191           0 :                  call endrun
         192             :               end if
         193             :            end do
         194             :            !------------------------------------------------------------------------------
         195             :            ! read N
         196             :            !------------------------------------------------------------------------------
         197           0 :            read(unit,*,iostat=istat) str
         198           0 :            if( istat /= 0 ) then
         199           0 :               write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         200           0 :               call endrun
         201             :            end if
         202           0 :            read(unit,*,iostat=istat) str
         203           0 :            if( istat /= 0 ) then
         204           0 :               write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         205           0 :               call endrun
         206             :            end if
         207           0 :            fmt = "(f7.2,2x,f7.2,2x,f9.3)"
         208           0 :            do i = 1,lmax
         209           0 :               read(unit,fmt,iostat=istat) waves(i), wavel(i), sigabs(i,4)
         210           0 :               if( istat /= 0 ) then
         211           0 :                  write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         212           0 :                  call endrun
         213             :               end if
         214             :            end do
         215             :
         216             :            !------------------------------------------------------------------------------
         217             :            ! read CO2
         218             :            !------------------------------------------------------------------------------
         219           0 :            read(unit,*,iostat=istat) str
         220           0 :            if( istat /= 0 ) then
         221           0 :               write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         222           0 :               call endrun
         223             :            end if
         224           0 :            read(unit,*,iostat=istat) str
         225           0 :            if( istat /= 0 ) then
         226           0 :               write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         227           0 :               call endrun
         228             :            end if
         229           0 :            fmt = "(f7.2,2x,f7.2,2x,f9.3)"
         230           0 :            do i = 1,lmax
         231           0 :               read(unit,fmt,iostat=istat) waves(i), wavel(i), sigabs(i,5)
         232           0 :               if( istat /= 0 ) then
         233           0 :                  write(iulog,*) 'jeuv_init: failed to read CO2 data ',trim(locfn),'; error = ',istat
         234           0 :                  call endrun
         235             :               end if
         236             :            end do
         237             :
         238           0 :            close( unit )
         239             :
         240             :            !------------------------------------------------------------------------------
         241             :            ! unit transfer for absorption cross sections
         242             :            ! from Megabarns to cm^2
         243             :            !------------------------------------------------------------------------------
         244           0 :            sigabs(:,:) = sigabs(:,:)*1.e-18_r8
         245             :
         246             :            !------------------------------------------------------------------------------
         247             :            ! read photoelectron ionization/dissociation/excitation branching ratio
         248             :            !------------------------------------------------------------------------------
         249           0 :            call getfil( electron_file, locfn, 0 )
         250           0 :            open( unit, file = trim(locfn), status='UNKNOWN', iostat=istat )
         251           0 :            if( istat /= 0 ) then
         252           0 :               write(iulog,*) 'jeuv_init: failed to open ',trim(locfn),'; error = ',istat
         253           0 :               call endrun
         254             :            end if
         255             :            !------------------------------------------------------------------------------
         256             :            ! read O
         257             :            !------------------------------------------------------------------------------
         258           0 :            read(unit,*,iostat=istat) str
         259           0 :            if( istat /= 0 ) then
         260           0 :               write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         261           0 :               call endrun
         262             :            end if
         263           0 :            read(unit,*,iostat=istat) str
         264           0 :            if( istat /= 0 ) then
         265           0 :               write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         266           0 :               call endrun
         267             :            end if
         268           0 :            fmt="(f7.2,2x,f7.2,5(2x,f8.3))"
         269           0 :            do i = 1,lmax
         270           0 :               read(unit,fmt,iostat=istat) waves(i), wavel(i), (branch_e(i,j,1),j=1,5)
         271           0 :               if( istat /= 0 ) then
         272           0 :                  write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         273           0 :                  call endrun
         274             :               end if
         275             :            end do
         276             :            !------------------------------------------------------------------------------
         277             :            ! read O2
         278             :            !------------------------------------------------------------------------------
         279           0 :            read(unit,*,iostat=istat) str
         280           0 :            if( istat /= 0 ) then
         281           0 :               write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         282           0 :               call endrun
         283             :            end if
         284           0 :            read(unit,*,iostat=istat) str
         285           0 :            if( istat /= 0 ) then
         286           0 :               write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         287           0 :               call endrun
         288             :            end if
         289           0 :            fmt = "(f7.2,2x,f7.2,6(2x,f8.3))"
         290           0 :            do i = 1,lmax
         291           0 :               read(unit,fmt,iostat=istat) waves(i), wavel(i), (branch_e(i,j,2),j=1,6)
         292           0 :               if( istat /= 0 ) then
         293           0 :                  write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         294           0 :                  call endrun
         295             :               end if
         296             :            end do
         297             :            !------------------------------------------------------------------------------
         298             :            ! read N2
         299             :            !------------------------------------------------------------------------------
         300           0 :            read(unit,*,iostat=istat) str
         301           0 :            if( istat /= 0 ) then
         302           0 :               write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         303           0 :               call endrun
         304             :            end if
         305           0 :            read(unit,*,iostat=istat) str
         306           0 :            if( istat /= 0 ) then
         307           0 :               write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         308           0 :               call endrun
         309             :            end if
         310           0 :            do i = 1,lmax
         311           0 :               read(unit,fmt,iostat=istat) waves(i), wavel(i), (branch_e(i,j,3),j=1,6)
         312           0 :               if( istat /= 0 ) then
         313           0 :                  write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
         314           0 :                  call endrun
         315             :               end if
         316             :            end do
         317             :
         318           0 :            close( unit )
         319             :
         320           0 :            call freeunit( unit )
         321             :         end if master
         322             :
         323           0 :         call mpi_bcast(sigabs,   lmax*nspecies,   mpi_real8, masterprocid, mpicom, ierr)
         324           0 :         call mpi_bcast(branch_p, lmax*nstat*nmaj, mpi_real8, masterprocid, mpicom, ierr)
         325           0 :         call mpi_bcast(branch_e, lmax*nstat*nmaj, mpi_real8, masterprocid, mpicom, ierr)
         326           0 :         call mpi_bcast(waves,    lmax,            mpi_real8, masterprocid, mpicom, ierr)
         327           0 :         call mpi_bcast(wavel,    lmax,            mpi_real8, masterprocid, mpicom, ierr)
         328             :
         329           0 :         wc(:)     = .1_r8*(waves(:) + wavel(:))*0.5_r8          ! A to nm
         330           0 :         energy(:) = heat_eff_fac*hc/wc(:)
         331             :
         332             :         end subroutine jeuv_init
         333             :
         334           0 :         subroutine jeuv( nlev, zen, occ, o2cc, n2cc, &
         335           0 :                          zkm, euv_prates)
         336             : !==============================================================================
         337             : !   Purpose:
         338             : !      Calculate euv photolysis/ionization rates (in s-1)
         339             : !==============================================================================
         340             : !   Arguments:
         341             : !       nlev: number of model vertical levels
         342             : !       zen:  solar zenith angle in degree
         343             : !       occ:  atmoic oxygen number density (#/cm3)
         344             : !       o2cc: molecular oxygen number density (#/cm3)
         345             : !       n2cc: molecular nitrogen number density (#/cm3)
         346             : !       zkm:  altitude of model levels in KM
         347             : !       euv_prates: array for EUV photolysis/ionization rates
         348             : !==============================================================================
         349             : !   Approach:
         350             : !       call sphers
         351             : !            input: zenith angle
         352             : !            output: dsdh and nid used in slant column routine
         353             : !       call slant_col
         354             : !            input: dsdh and nid, etc
         355             : !            output: slant column density
         356             : !       calculate photon production rates and photoelectron production rates
         357             : !==============================================================================
         358             : !  Photolysis/ionization considered in EUV calculation
         359             : !
         360             : !  O + hv --> O+ (4S) + e*                        ! J1
         361             : !  O + hv --> O+ (2D) + e*                        ! J2
         362             : !  O + hv --> O+ (2P) + e*                        ! J3
         363             : !  N (4S) + hv --> N+ + e*                        ! J4
         364             : !  O2 + hv --> O2+ + e*                           ! J5
         365             : !  N2 + hv --> N2+ + e*                           ! J6
         366             : !  O2 + hv --> O + O+(4S) + e*                    ! J7
         367             : !  O2 + hv --> O + O+(2D) + e*                    ! J8
         368             : !  O2 + hv --> O + O+(2P) + e*                    ! J9
         369             : !  N2 + hv --> N (4S) + N+ + e*                   ! J10
         370             : !  N2 + hv --> N (2D) + N+ + e*                   ! J11
         371             : !  O2 + hv --> O (3P) + O (3P)                    ! J12
         372             : !  N2 + hv --> N (4S) + N (2D)                    ! J13
         373             : !
         374             : !  O + e* --> O+ (4S) + e* + e                    ! J14
         375             : !  O + e* --> O+ (2D) + e*+ e                     ! J15
         376             : !  O + e* --> O+ (2P) + e*+ e                     ! J16
         377             : !  O2 + e* --> O2+ + e*+ e                        ! J17
         378             : !  N2 + e*--> N2+ + e*+ e                         ! J18
         379             : !  O2 + e* --> O + O+(4S) + e*+ e                 ! J19
         380             : !  O2 + e*--> O + O+(2D) + e*+ e                  ! J20
         381             : !  O2 + e* --> O + O+(2P) + e*+ e                 ! J21
         382             : !  N2 + e* --> N (4S) + N+ + e*+ e                ! J22
         383             : !  N2 + e* --> N (2D) + N+ + e*+ e                ! J23
         384             : !  O2 + e* --> O (3P) + O (3P) + e*               ! J24
         385             : !  N2 + e* --> N (4S) + N (2D) + e*               ! J25
         386             : !==============================================================================
         387             :
         388             :         use mo_jshort,   only : sphers, slant_col
         389             :
         390             :         implicit none
         391             :
         392             : !------------------------------------------------------------------------------
         393             : !       dummy arguments
         394             : !------------------------------------------------------------------------------
         395             :         integer, intent(in)     :: nlev                        ! model vertical levels
         396             :         real(r8), intent(in)    :: zen                         ! Zenith  angle in degree
         397             :         real(r8), intent(in)    :: occ(nlev)                   ! atmic oxygen number density (#/cm3)
         398             :         real(r8), intent(in)    :: o2cc(nlev)                  ! Molecular oxygen number density (#/cm3)
         399             :         real(r8), intent(in)    :: n2cc(nlev)                  ! molecular nitrogen number density(#/cm3)
         400             :         real(r8), intent(in)    :: zkm(nlev)                   ! Altitude, km,from top to bottom
         401             :         real(r8), intent(out)   :: euv_prates(:,:)             ! EUV photolysis/ionization rates (s-1)
         402             :
         403             : !------------------------------------------------------------------------------
         404             : !       local variables
         405             : !------------------------------------------------------------------------------
         406             :         real(r8), parameter :: km2cm = 1.e5_r8
         407             :         integer  :: l, k, m, n              ! loop indecies
         408             :         real(r8) :: tau(lmax)               ! wavelength dependant optical depth
         409           0 :         real(r8) :: delz(nlev)              ! layer thickness (cm)
         410           0 :         real(r8) :: scol(nlev,nmaj)         ! major species' (O,O2,N2) Slant Column Density
         411           0 :         real(r8) :: nflux(nlev,lmax)
         412             :         real(r8) :: wrk(nmaj)               ! temporary array for photoabsorption rate
         413           0 :         real(r8) :: absorp(nlev,lmax)       ! temporary array for photoabsorption rate
         414           0 :         real(r8) :: ioniz(nlev,lmax)        ! temporary array for photoionization rate
         415           0 :         real(r8) :: dsdh(0:nlev,nlev)       ! Slant path of direct beam through each layer
         416             :                                             ! crossed  when travelling from the top of the atmosphere
         417             :                                             ! to layer i; dsdh(i,j), i = 0..NZ-1, j = 1..NZ-1
         418           0 :         integer :: nid(0:nlev)              ! Number of layers crossed by the direct
         419             :                                             ! beam when travelling from the top of the
         420             :                                             ! atmosphere to layer i; NID(i), i = 0..NZ-1
         421           0 :         real(r8) :: p_photon(nlev,nstat,nspecies)   !  photoionization/dissociation rates(s-1) (O,O2,N2,N)
         422           0 :         real(r8) :: p_electron(nlev,nstat,nmaj) !  photoelectron ionization/dissociation rates(s-1) (O,O2,N2)
         423             :
         424           0 :         real(r8) :: prates(nlev,neuv)
         425             : !------------------------------------------------------------------------------
         426             : ! zero arrays
         427             : !------------------------------------------------------------------------------
         428           0 :         p_photon(:,:,:)   = 0._r8
         429           0 :         p_electron(:,:,:) = 0._r8
         430             :
         431           0 :         call sphers( nlev, zkm, zen, dsdh, nid )
         432           0 :         delz(1:nlev-1) = km2cm*(zkm(1:nlev-1) - zkm(2:nlev))
         433           0 :         call slant_col( nlev, delz, dsdh, nid, occ, scol )
         434           0 :         call slant_col( nlev, delz, dsdh, nid, o2cc, scol(1,2) )
         435           0 :         call slant_col( nlev, delz, dsdh, nid, n2cc, scol(1,3) )
         436             :
         437             : !------------------------------------------------------------------------------
         438             : ! The calculation is in the order from model bottom to model top
         439             : ! because scol array is in this order.
         440             : !------------------------------------------------------------------------------
         441           0 :         do k = 1,nlev
         442           0 :            wrk(:) = scol(k,:)
         443           0 :            tau(:) = matmul( sigabs(:,:nmaj),wrk )
         444           0 :            where( tau(:) < 20._r8 )
         445           0 :               nflux(k,:) = solar_etf(:) * exp( -tau(:) )
         446             :            elsewhere
         447             :               nflux(k,:) = 0._r8
         448             :            endwhere
         449             :         end do
         450             :
         451             : !------------------------------------------------------------------------------
         452             : ! remember occ,o2cc and n2cc is from top to bottom
         453             : !------------------------------------------------------------------------------
         454           0 :         do m = 1,nspecies
         455           0 :            do l = 1,lmax
         456           0 :               absorp(:,l) = sigabs(l,m) * nflux(:,l)
         457             :            end do
         458           0 :            if( m <= nmaj ) then
         459           0 :               do l = 1,lmax
         460           0 :                  ioniz(:,l) = absorp(:,l) * branch_p(l,1,m)
         461             :               end do
         462           0 :               do n = 1,nstat
         463           0 :                  p_photon(:,n,m)   = matmul( absorp,branch_p(:,n,m) )
         464           0 :                  p_electron(:,n,m) = matmul( ioniz,branch_e(:,n,m) )
         465             :               end do
         466             :            else
         467           0 :               p_photon(:,1,m) = matmul( nflux,sigabs(:,m) )
         468             :            end if
         469             :         end do
         470             :
         471             : !------------------------------------------------------------------------------
         472             : ! set photolysis/ionization rate for each reaction
         473             : !------------------------------------------------------------------------------
         474           0 :        prates(:,1)  = p_photon(:,2,1)
         475           0 :        prates(:,2)  = p_photon(:,3,1)
         476           0 :        prates(:,3)  = p_photon(:,4,1)
         477           0 :        prates(:,4)  = p_photon(:,1,4)
         478           0 :        prates(:,5)  = p_photon(:,2,2) + p_photon(:,3,2)
         479           0 :        prates(:,6)  = p_photon(:,2,3) + p_photon(:,3,3)
         480           0 :        prates(:,7)  = .54_r8*p_photon(:,4,2)
         481           0 :        prates(:,8)  = .24_r8*p_photon(:,4,2)
         482           0 :        prates(:,9)  = .22_r8*p_photon(:,4,2)
         483           0 :        prates(:,10) = .2_r8*p_photon(:,4,3)
         484           0 :        prates(:,11) = .8_r8*p_photon(:,4,3)
         485           0 :        prates(:,12) = p_photon(:,5,2)
         486           0 :        prates(:,13) = p_photon(:,5,3)
         487           0 :        prates(:,14) = p_electron(:,2,1)
         488           0 :        prates(:,15) = p_electron(:,3,1)
         489           0 :        prates(:,16) = p_electron(:,4,1)
         490           0 :        prates(:,17) = p_electron(:,2,2) + p_electron(:,3,2)
         491           0 :        prates(:,18) = p_electron(:,2,3) + p_electron(:,3,3)
         492           0 :        prates(:,19) = .54_r8*p_electron(:,4,2)
         493           0 :        prates(:,20) = .24_r8*p_electron(:,4,2)
         494           0 :        prates(:,21) = .22_r8*p_electron(:,4,2)
         495           0 :        prates(:,22) = .2_r8*p_electron(:,4,3)
         496           0 :        prates(:,23) = .8_r8*p_electron(:,4,3)
         497           0 :        prates(:,24) = p_electron(:,5,2)
         498           0 :        prates(:,25) = p_electron(:,5,3)
         499           0 :        prates(:,26) = p_photon(:,1,5)
         500             :
         501           0 :        do m = 1,neuv
         502           0 :           euv_prates(:,m) = prates(nlev:1:-1,m)
         503             :        enddo
         504             :
         505           0 :        end subroutine jeuv
         506             :
         507           0 :        subroutine heuv( nlev, zen, occ, o2cc, n2cc, &
         508           0 :                         o_vmr, o2_vmr, n2_vmr, cparg, mw, &
         509           0 :                         zkm, euv_hrates, kbot )
         510             : !==============================================================================
         511             : !   Purpose:
         512             : !      Calculate euv photolysis heating rates
         513             : !==============================================================================
         514             : !   Arguments:
         515             : !       nlev: number of model vertical levels
         516             : !       zen:  solar zenith angle in degree
         517             : !       occ:  atmoic oxygen number density (#/cm3)
         518             : !       o2cc: molecular oxygen number density (#/cm3)
         519             : !       n2cc: molecular nitrogen number density (#/cm3)
         520             : !       zkm:  altitude of model levels in KM
         521             : !       euv_prates: array for EUV photolysis/ionization rates
         522             : !==============================================================================
         523             : !   Approach:
         524             : !       call sphers
         525             : !            input: zenith angle
         526             : !            output: dsdh and nid used in slant column routine
         527             : !       call slant_col
         528             : !            input: dsdh and nid, etc
         529             : !            output: slant column density
         530             : !       calculate photon production rates and photoelectron production rates
         531             : !==============================================================================
         532             : !  Photolysis/ionization considered in EUV heating rate calculation
         533             : !  O + hv --> O+ (4S) + e*                        ! J1
         534             : !  O + hv --> O+ (2D) + e*                        ! J2
         535             : !  O + hv --> O+ (2P) + e*                        ! J3
         536             : !  N (4S) + hv --> N+ + e*                        ! J4
         537             : !  O2 + hv --> O2+ + e*                           ! J5
         538             : !  N2 + hv --> N2+ + e*                           ! J6
         539             : !  O2 + hv --> O + O+(4S) + e*                    ! J7
         540             : !  O2 + hv --> O + O+(2D) + e*                    ! J8
         541             : !  O2 + hv --> O + O+(2P) + e*                    ! J9
         542             : !  N2 + hv --> N (4S) + N+ + e*                   ! J10
         543             : !  N2 + hv --> N (2D) + N+ + e*                   ! J11
         544             : !  O2 + hv --> O (3P) + O (3P)                    ! J12
         545             : !  N2 + hv --> N (4S) + N (2D)                    ! J13
         546             : !==============================================================================
         547             :
         548             :         use mo_jshort,     only : sphers, slant_col
         549             :         use physconst,     only : avogad
         550             :
         551             :         implicit none
         552             :
         553             : !------------------------------------------------------------------------------
         554             : !       dummy arguments
         555             : !------------------------------------------------------------------------------
         556             :         integer, intent(in)     :: nlev                        ! model vertical levels
         557             :         integer, intent(in)     :: kbot                        ! heating vertical levels
         558             :         real(r8), intent(in)    :: zen                         ! zenith  angle (degrees)
         559             :         real(r8), intent(in)    :: occ(nlev)                   ! atomic oxygen number density (molecules/cm^3)
         560             :         real(r8), intent(in)    :: o2cc(nlev)                  ! molecular oxygen concentration (molecules/cm^3)
         561             :         real(r8), intent(in)    :: n2cc(nlev)                  ! molecular nitrogen concentration (molecules/cm^3)
         562             :         real(r8), intent(in)    :: o_vmr(nlev)                 ! atomic oxygen concentration (mol/mol)
         563             :         real(r8), intent(in)    :: o2_vmr(nlev)                ! molecular oxygen concentration (mol/mol)
         564             :         real(r8), intent(in)    :: n2_vmr(nlev)                ! molecular nitrogen concentration (mol/mol)
         565             :         real(r8), intent(in)    :: zkm(nlev)                   ! midpoint geopotential (km)
         566             :         real(r8), intent(in)    :: cparg(nlev)                 ! specific heat capacity
         567             :         real(r8), intent(in)    :: mw(nlev)                    ! atm mean mass
         568             :         real(r8), intent(out)   :: euv_hrates(:)               ! euv heating rates
         569             :
         570             : !------------------------------------------------------------------------------
         571             : !       local variables
         572             : !------------------------------------------------------------------------------
         573             :         real(r8), parameter :: km2cm = 1.e5_r8
         574             :         integer  :: l, k, k1, m, n          ! indicies
         575             :         real(r8) :: tau(lmax)               ! wavelength dependant optical depth
         576           0 :         real(r8) :: delz(nlev)              ! layer thickness (cm)
         577           0 :         real(r8) :: hfactor(kbot)
         578           0 :         real(r8) :: scol(nlev,nmaj)         ! major species' (O,O2,N2) Slant Column Density
         579           0 :         real(r8) :: nflux(kbot,lmax)
         580           0 :         real(r8) :: prates(kbot,13)         ! working photorates array
         581             :         real(r8) :: wrk(nmaj)               ! temporary array for photoabsorption rate
         582           0 :         real(r8) :: absorp(kbot,lmax)       ! temporary array for photoabsorption rate
         583           0 :         real(r8) :: dsdh(0:nlev,nlev)       ! Slant path of direct beam through each layer
         584             :                                             ! crossed  when travelling from the top of the atmosphere
         585             :                                             ! to layer i; dsdh(i,j), i = 0..NZ-1, j = 1..NZ-1
         586           0 :         integer :: nid(0:nlev)              ! Number of layers crossed by the direct
         587             :                                             ! beam when travelling from the top of the
         588             :                                             ! atmosphere to layer i; NID(i), i = 0..NZ-1
         589           0 :         real(r8) :: p_photon(kbot,nstat,nspecies)   !  photoionization/dissociation rates(s-1) (O,O2,N2,N)
         590             :
         591             : !------------------------------------------------------------------------------
         592             : ! zero arrays
         593             : !------------------------------------------------------------------------------
         594           0 :         p_photon(:,:,:)   = 0._r8
         595             :
         596           0 :         call sphers( nlev, zkm, zen, dsdh, nid )
         597           0 :         delz(1:nlev-1) = km2cm*(zkm(1:nlev-1) - zkm(2:nlev))
         598           0 :         call slant_col( nlev, delz, dsdh, nid, occ, scol )
         599           0 :         call slant_col( nlev, delz, dsdh, nid, o2cc, scol(1,2) )
         600           0 :         call slant_col( nlev, delz, dsdh, nid, n2cc, scol(1,3) )
         601             :
         602             : !------------------------------------------------------------------------------
         603             : ! The calculation is in the order from model bottom to model top
         604             : ! because scol array is in this order.
         605             : !------------------------------------------------------------------------------
         606           0 :         do k = 1,kbot
         607           0 :            k1 = nlev - kbot + k
         608           0 :            wrk(:) = scol(k1,:)
         609           0 :            tau(:) = matmul( sigabs(:,:nmaj),wrk )
         610           0 :            where( tau(:) < 20._r8 )
         611           0 :               nflux(k,:) = energy(:) * solar_etf(:) * exp( -tau(:) )
         612             :            elsewhere
         613           0 :               nflux(k,:) = 0._r8
         614             :            endwhere
         615             :         end do
         616             :
         617             : !------------------------------------------------------------------------------
         618             : ! remember occ,o2cc and n2cc is from top to bottom
         619             : !------------------------------------------------------------------------------
         620           0 :         do m = 1,nspecies
         621           0 :            do l = 1,lmax
         622           0 :               absorp(:,l) = sigabs(l,m) * nflux(:,l)
         623             :            end do
         624           0 :            if( m <= nmaj ) then
         625           0 :               do n = 1,nstat
         626           0 :                  p_photon(:,n,m) = matmul( absorp,branch_p(:,n,m) )
         627             :               end do
         628             :            else
         629           0 :               p_photon(:,1,m) = matmul( nflux,sigabs(:,m) )
         630             :            end if
         631             :         end do
         632             :
         633             : !------------------------------------------------------------------------------
         634             : ! use photolysis rates to compute heating rates
         635             : ! only EUV photolysis reactions in the mechanism contribute to heating
         636             : !------------------------------------------------------------------------------
         637           0 :        prates = 0._r8
         638             :
         639           0 :        if (do_heating(1)) prates(:,1)  = p_photon(:,2,1)
         640           0 :        if (do_heating(2)) prates(:,2)  = p_photon(:,3,1)
         641           0 :        if (do_heating(3)) prates(:,3)  = p_photon(:,4,1)
         642             :
         643           0 :        if (do_heating(5)) prates(:,5)  = p_photon(:,2,2) + p_photon(:,3,2)
         644           0 :        if (do_heating(6)) prates(:,6)  = p_photon(:,2,3) + p_photon(:,3,3)
         645           0 :        if (do_heating(7)) prates(:,7)  = .54_r8*p_photon(:,4,2)
         646           0 :        if (do_heating(8)) prates(:,8)  = .24_r8*p_photon(:,4,2)
         647           0 :        if (do_heating(9)) prates(:,9)  = .22_r8*p_photon(:,4,2)
         648           0 :        if (do_heating(10)) prates(:,10) = .2_r8*p_photon(:,4,3)
         649           0 :        if (do_heating(11)) prates(:,11) = .8_r8*p_photon(:,4,3)
         650           0 :        if (do_heating(12)) prates(:,12) = p_photon(:,5,2)
         651           0 :        if (do_heating(13)) prates(:,13) = p_photon(:,5,3)
         652           0 :        hfactor(:)   = avogad/(cparg(:kbot)*mw(:kbot))
         653           0 :        euv_hrates(kbot+1:nlev) = 0._r8
         654           0 :        euv_hrates(:kbot) = ((prates(kbot:1:-1,1) + prates(kbot:1:-1,2) + prates(kbot:1:-1,3))*o_vmr(:kbot) &
         655             :                      + (prates(kbot:1:-1,5) + prates(kbot:1:-1,7) + prates(kbot:1:-1,8) &
         656             :                         + prates(kbot:1:-1,9) + prates(kbot:1:-1,12))*o2_vmr(:kbot) &
         657             :                      + (prates(kbot:1:-1,6) + prates(kbot:1:-1,10) &
         658           0 :                         + prates(kbot:1:-1,11) + prates(kbot:1:-1,13))*n2_vmr(:kbot))*hfactor(:)
         659             :
         660           0 :        end subroutine heuv
         661             :
         662             :        end module mo_jeuv

Generated by: LCOV version 1.14