LCOV - code coverage report
Current view: | top level - chemistry/mozart - upper_bc.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 71 | 227 | 31.3 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 4 | 6 | 66.7 %

              Line data    Source code
           1             : module upper_bc
           2             :
           3             : !---------------------------------------------------------------------------------
           4             : ! Module to compute the upper boundary conditions for temperature (dry static energy)
           5             : ! and trace gases. Uses the MSIS model, and SNOE and TIME GCM and general prescribed UBC data.
           6             : !---------------------------------------------------------------------------------
           7             :
           8             :   use shr_kind_mod, only: r8 => shr_kind_r8
           9             :   use shr_kind_mod, only: cl => shr_kind_cl
          10             :   use shr_const_mod,only: grav   => shr_const_g,     &   ! gravitational constant (m/s^2)
          11             :                           kboltz => shr_const_boltz, &   ! Boltzmann constant
          12             :                           pi => shr_const_pi,        &   ! pi
          13             :                           rEarth => shr_const_rearth     ! Earth radius
          14             :   use ppgrid,       only: pcols, pver, pverp
          15             :   use constituents, only: pcnst
          16             :   use cam_logfile,  only: iulog
          17             :   use spmd_utils,   only: masterproc
          18             :   use ref_pres,     only: do_molec_diff, ptop_ref
          19             :   use shr_kind_mod, only: cx=>SHR_KIND_CX
          20             :   use cam_abortutils,only: endrun
          21             :   use cam_history,   only: addfld, horiz_only, outfld, fieldname_len
          22             :
          23             :   use upper_bc_file, only: upper_bc_file_readnl, upper_bc_file_specified, upper_bc_file_adv, upper_bc_file_get
          24             :   use infnan,        only: nan, assignment(=)
          25             :
          26             :   implicit none
          27             :   private
          28             :   save
          29             : !
          30             : ! Public interfaces
          31             : !
          32             :   public :: ubc_readnl         ! read namelist options for UBCs
          33             :   public :: ubc_init           ! global initialization
          34             :   public :: ubc_timestep_init  ! time step initialization
          35             :   public :: ubc_get_vals       ! get ubc values for this step
          36             :   public :: ubc_get_flxs       ! get ub fluxes for this step
          37             :   public :: ubc_fixed_conc     ! returns true for constituents that have fixed UBC
          38             :   public :: ubc_fixed_temp     ! true if temperature at upper boundary is fixed
          39             :
          40             :   character(len=64) :: ubc_specifier(pcnst) = 'NOTSET'
          41             :
          42             :   character(len=16) :: ubc_flds(pcnst) = 'NOTSET'
          43             :   character(len=16) :: ubc_file_spfr(pcnst) = ' '
          44             :   character(len=32) :: ubc_source(pcnst) = ' '
          45             :
          46             :   integer :: n_fixed_mmr=0
          47             :   integer :: n_fixed_vmr=0
          48             :
          49             :   integer, allocatable :: fixed_mmr_ndx(:)
          50             :   integer, allocatable :: fixed_vmr_ndx(:)
          51             :   real(r8), allocatable :: fixed_mmr(:)
          52             :   real(r8), allocatable :: fixed_vmr(:)
          53             :
          54             :   integer :: num_infile = 0
          55             :   integer :: num_fixed = 0
          56             :   character(len=2), parameter :: msis_flds(5) = &
          57             :        (/ 'H ','N ','O ','O2','T ' /)
          58             :   character(len=2), parameter :: tgcm_flds(1) = &
          59             :        (/ 'H2' /)
          60             :   character(len=2), parameter :: snoe_flds(1) = &
          61             :        (/ 'NO' /)
          62             :
          63             :   logical, protected :: ubc_fixed_temp =.false.
          64             :   logical :: msis_active =.false.
          65             :   logical :: tgcm_active =.false.
          66             :   logical :: snoe_active =.false.
          67             :
          68             : ! Namelist variables
          69             :   character(len=cl) :: snoe_ubc_file = 'NONE'
          70             :   real(r8)          :: t_pert_ubc  = 0._r8
          71             :   real(r8)          :: no_xfac_ubc = 1._r8
          72             :
          73             :   integer :: h_ndx=-1
          74             :   integer :: h_msis_ndx=-1, n_msis_ndx=-1, o_msis_ndx=-1, o2_msis_ndx=-1
          75             :
          76             :   character(len=cl) :: tgcm_ubc_file = 'NONE'
          77             :   integer           :: tgcm_ubc_cycle_yr = 0
          78             :   integer           :: tgcm_ubc_fixed_ymd = 0
          79             :   integer           :: tgcm_ubc_fixed_tod = 0
          80             :   character(len=32) :: tgcm_ubc_data_type = 'CYCLICAL'
          81             :
          82             :   logical :: apply_upper_bc = .false.
          83             :
          84             :   integer, allocatable :: file_spc_ndx(:)
          85             :   integer, allocatable :: spc_ndx(:)
          86             :   character(len=fieldname_len), allocatable :: hist_names(:)
          87             :
          88             : !================================================================================================
          89             : contains
          90             : !================================================================================================
          91             :
          92             : !-----------------------------------------------------------------------
          93             : !-----------------------------------------------------------------------
          94        1536 :   subroutine ubc_readnl(nlfile)
          95             :     use namelist_utils, only : find_group_name
          96             :     use spmd_utils, only : mpicom, masterprocid, mpi_character, mpi_integer, mpi_real8
          97             :     use string_utils, only : to_lower
          98             :
          99             :     character(len=*), intent(in) :: nlfile
         100             :     integer :: unitn, ierr, m, n, ndx_co, ndx_ar
         101             :
         102             :     character(len=*), parameter :: prefix = 'ubc_readnl: '
         103             :
         104             :     namelist /upper_bc_opts/ tgcm_ubc_file,tgcm_ubc_data_type,tgcm_ubc_cycle_yr,tgcm_ubc_fixed_ymd, &
         105             :                              tgcm_ubc_fixed_tod, snoe_ubc_file, no_xfac_ubc, t_pert_ubc
         106             :     namelist /upper_bc_opts/ ubc_specifier
         107             :
         108        1536 :     if (masterproc) then
         109             :        ! read namelist
         110           2 :        open( newunit=unitn, file=trim(nlfile), status='old' )
         111           2 :        call find_group_name(unitn, 'upper_bc_opts', status=ierr)
         112           2 :        if (ierr == 0) then
         113           0 :           read(unitn, upper_bc_opts, iostat=ierr)
         114           0 :           if (ierr /= 0) then
         115           0 :              call endrun(prefix//'upper_bc_opts: ERROR reading namelist')
         116             :           end if
         117             :        end if
         118           2 :        close(unitn)
         119             :
         120             :        ! log the UBC options
         121           2 :        write(iulog,*) prefix//'tgcm_ubc_file = '//trim(tgcm_ubc_file)
         122           2 :        write(iulog,*) prefix//'tgcm_ubc_data_type = '//trim(tgcm_ubc_data_type)
         123           2 :        write(iulog,*) prefix//'tgcm_ubc_cycle_yr = ', tgcm_ubc_cycle_yr
         124           2 :        write(iulog,*) prefix//'tgcm_ubc_fixed_ymd = ', tgcm_ubc_fixed_ymd
         125           2 :        write(iulog,*) prefix//'tgcm_ubc_fixed_tod = ', tgcm_ubc_fixed_tod
         126           2 :        write(iulog,*) prefix//'snoe_ubc_file = '//trim(snoe_ubc_file)
         127           2 :        write(iulog,*) prefix//'t_pert_ubc = ', t_pert_ubc
         128           2 :        write(iulog,*) prefix//'no_xfac_ubc = ', no_xfac_ubc
         129           2 :        write(iulog,*) prefix//'ubc_specifier : '
         130             :
         131           2 :        n=1
         132           2 :        m=1
         133           2 :        do while(ubc_specifier(n)/='NOTSET')
         134           0 :           write(iulog,'(i4,a)') n,'  '//trim(ubc_specifier(n))
         135             :
         136           0 :           ndx_ar = index(ubc_specifier(n),'->')
         137             :
         138           0 :           if (ndx_ar<1) then
         139           0 :              call endrun(prefix//'ubc_specifier "'//trim(ubc_specifier(n))//'" must include "->"')
         140             :           endif
         141             :
         142           0 :           ubc_source(n) = trim(to_lower(adjustl(ubc_specifier(n)(ndx_ar+2:))))
         143             :
         144           0 :           if (trim(ubc_source(n))=='ubc_file') then
         145           0 :              ubc_file_spfr(m) = trim(ubc_specifier(n)(:ndx_ar-1))
         146           0 :              m=m+1
         147             :           endif
         148           0 :           if (index(ubc_source(n),'mmr')>0) then
         149           0 :              n_fixed_mmr=n_fixed_mmr+1
         150           0 :           else if (index(ubc_source(n),'vmr')>0) then
         151           0 :              n_fixed_vmr=n_fixed_vmr+1
         152             :           end if
         153             :
         154           0 :           ndx_co = index(ubc_specifier(n),':')
         155             :
         156           0 :           if (ndx_co>0) then
         157           0 :              ubc_flds(n) = ubc_specifier(n)(:ndx_co-1)
         158             :           else
         159           0 :              ubc_flds(n) = ubc_specifier(n)(:ndx_ar-1)
         160             :           end if
         161             :
         162           0 :           n=n+1
         163             :        end do
         164           2 :        num_fixed=n-1
         165           2 :        num_infile=m-1
         166             :     end if
         167             :
         168             :
         169             :     ! broadcast to all MPI tasks
         170        1536 :     call mpi_bcast(num_fixed, 1, mpi_integer, masterprocid, mpicom, ierr)
         171        1536 :     if (ierr /= 0) call endrun(prefix//'mpi_bcast error : num_fixed')
         172        1536 :     call mpi_bcast(num_infile, 1, mpi_integer, masterprocid, mpicom, ierr)
         173        1536 :     if (ierr /= 0) call endrun(prefix//'mpi_bcast error : num_infile')
         174        1536 :     call mpi_bcast(n_fixed_mmr, 1, mpi_integer, masterprocid, mpicom, ierr)
         175        1536 :     if (ierr /= 0) call endrun(prefix//'mpi_bcast error : n_fixed_mmr')
         176        1536 :     call mpi_bcast(n_fixed_vmr, 1, mpi_integer, masterprocid, mpicom, ierr)
         177        1536 :     if (ierr /= 0) call endrun(prefix//'mpi_bcast error : n_fixed_vmr')
         178        1536 :     call mpi_bcast(tgcm_ubc_file, len(tgcm_ubc_file), mpi_character, masterprocid, mpicom, ierr)
         179        1536 :     if (ierr /= 0) call endrun(prefix//'mpi_bcast error : tgcm_ubc_file')
         180        1536 :     call mpi_bcast(tgcm_ubc_data_type, len(tgcm_ubc_data_type),mpi_character, masterprocid, mpicom, ierr)
         181        1536 :     if (ierr /= 0) call endrun(prefix//'mpi_bcast error : tgcm_ubc_data_type')
         182        1536 :     call mpi_bcast(tgcm_ubc_cycle_yr, 1, mpi_integer, masterprocid, mpicom, ierr)
         183        1536 :     if (ierr /= 0) call endrun(prefix//'mpi_bcast error : tgcm_ubc_cycle_yr')
         184        1536 :     call mpi_bcast(tgcm_ubc_fixed_ymd, 1, mpi_integer, masterprocid, mpicom, ierr)
         185        1536 :     if (ierr /= 0) call endrun(prefix//'mpi_bcast error : tgcm_ubc_fixed_ymd')
         186        1536 :     call mpi_bcast(tgcm_ubc_fixed_tod, 1, mpi_integer, masterprocid, mpicom, ierr)
         187        1536 :     if (ierr /= 0) call endrun(prefix//'mpi_bcast error : tgcm_ubc_fixed_tod')
         188        1536 :     call mpi_bcast(snoe_ubc_file, len(snoe_ubc_file), mpi_character, masterprocid, mpicom, ierr)
         189        1536 :     if (ierr /= 0) call endrun(prefix//'mpi_bcast error : snoe_ubc_file')
         190        1536 :     call mpi_bcast(t_pert_ubc, 1, mpi_real8, masterprocid, mpicom, ierr)
         191        1536 :     if (ierr /= 0) call endrun(prefix//'mpi_bcast error : t_pert_ubc')
         192        1536 :     call mpi_bcast(no_xfac_ubc,1, mpi_real8, masterprocid, mpicom, ierr)
         193        1536 :     if (ierr /= 0) call endrun(prefix//'mpi_bcast error : no_xfac_ubc')
         194        1536 :     call mpi_bcast(ubc_specifier, pcnst*len(ubc_specifier(1)), mpi_character,masterprocid, mpicom, ierr)
         195        1536 :     if (ierr /= 0) call endrun(prefix//'mpi_bcast error : ubc_specifier')
         196        1536 :     call mpi_bcast(ubc_flds,      pcnst*len(ubc_flds(1)),      mpi_character,masterprocid, mpicom, ierr)
         197        1536 :     if (ierr /= 0) call endrun(prefix//'mpi_bcast error : ubc_flds')
         198        1536 :     call mpi_bcast(ubc_file_spfr, pcnst*len(ubc_file_spfr(1)), mpi_character,masterprocid, mpicom, ierr)
         199        1536 :     if (ierr /= 0) call endrun(prefix//'mpi_bcast error : ubc_file_spfr')
         200        1536 :     call mpi_bcast(ubc_source,    pcnst*len(ubc_source(1)),    mpi_character,masterprocid, mpicom, ierr)
         201        1536 :     if (ierr /= 0) call endrun(prefix//'mpi_bcast error : ubc_source')
         202             :
         203        1536 :     apply_upper_bc = num_fixed>0
         204             :
         205        1536 :     if (apply_upper_bc) then
         206           0 :        call upper_bc_file_readnl(nlfile)
         207             :     end if
         208             :
         209        1536 :   end subroutine ubc_readnl
         210             :
         211             : !===============================================================================
         212             :
         213        1536 :   subroutine ubc_init()
         214             : !-----------------------------------------------------------------------
         215             : ! Initialization of time independent fields for the upper boundary condition
         216             : ! Calls initialization routine for MSIS, TGCM and SNOE
         217             : !-----------------------------------------------------------------------
         218             :     use mo_tgcm_ubc, only: tgcm_ubc_inti
         219             :     use mo_snoe,     only: snoe_inti
         220             :     use mo_msis_ubc, only: msis_ubc_inti
         221             :     use constituents,only: cnst_get_ind
         222             :     use upper_bc_file, only: upper_bc_file_init
         223             :
         224             :     !---------------------------Local workspace-----------------------------
         225             :     logical, parameter :: zonal_avg = .false.
         226             :     integer :: m, mm, ierr
         227             :     integer :: mmrndx, vmrndx, m_mmr, m_vmr
         228             :
         229             :     real(r8) :: val
         230             :     character(len=32) :: str
         231             :     character(len=*), parameter :: prefix = 'ubc_init: '
         232             :     !-----------------------------------------------------------------------
         233             :
         234        1536 :     call cnst_get_ind('H', h_ndx, abort=.false.) ! for H fluxes UBC (WACCMX)
         235             :
         236        1536 :     if (.not.apply_upper_bc) return
         237             :
         238           0 :     if (num_infile>0) then
         239           0 :        call upper_bc_file_init( ubc_file_spfr(:num_infile) )
         240             :     endif
         241             :
         242             :     ! possible MSIS, TGCM, SNOE and ubc_file inputs
         243             :
         244           0 :     mm=1
         245             :
         246           0 :     allocate(hist_names(num_fixed), stat=ierr)
         247           0 :     if (ierr /= 0) call endrun(prefix//'allocate error : hist_names')
         248           0 :     allocate(spc_ndx(num_fixed), stat=ierr)
         249           0 :     if (ierr /= 0) call endrun(prefix//'allocate error : spc_ndx')
         250           0 :     spc_ndx=-1
         251           0 :     if (num_infile>0) then
         252           0 :        allocate(file_spc_ndx(num_infile), stat=ierr)
         253           0 :        if (ierr /= 0) call endrun(prefix//'allocate error : file_spc_ndx')
         254           0 :        file_spc_ndx=-1
         255             :     end if
         256           0 :     if (n_fixed_mmr>0) then
         257           0 :        allocate(fixed_mmr_ndx(n_fixed_mmr), stat=ierr)
         258           0 :        if (ierr /= 0) call endrun(prefix//'allocate error : fixed_mmr_ndx')
         259           0 :        allocate(fixed_mmr(n_fixed_mmr), stat=ierr)
         260           0 :        if (ierr /= 0) call endrun(prefix//'allocate error : fixed_mmr')
         261           0 :        fixed_mmr_ndx=-1
         262           0 :        fixed_mmr = nan
         263             :     end if
         264           0 :     if (n_fixed_vmr>0) then
         265           0 :        allocate(fixed_vmr_ndx(n_fixed_vmr), stat=ierr)
         266           0 :        if (ierr /= 0) call endrun(prefix//'allocate error : fixed_vmr_ndx')
         267           0 :        allocate(fixed_vmr(n_fixed_vmr), stat=ierr)
         268           0 :        if (ierr /= 0) call endrun(prefix//'allocate error : fixed_vmr')
         269           0 :        fixed_vmr_ndx=-1
         270           0 :        fixed_vmr = nan
         271             :     end if
         272             :
         273           0 :     m_mmr = 0
         274           0 :     m_vmr = 0
         275             :
         276           0 :     do m = 1,num_fixed
         277           0 :        hist_names(m) = trim(ubc_flds(m))//'_UBC'
         278           0 :        if (ubc_flds(m)=='T') then
         279           0 :           ubc_fixed_temp=.true.
         280           0 :           spc_ndx(m) = -1
         281           0 :           call addfld(hist_names(m), horiz_only, 'I', 'K', trim(ubc_flds(m))//' at upper boundary' )
         282             :        else
         283           0 :           call cnst_get_ind(ubc_flds(m), spc_ndx(m), abort=.true.)
         284           0 :           call addfld(hist_names(m), horiz_only, 'I', 'kg/kg', trim(ubc_flds(m))//' at upper boundary' )
         285             :        end if
         286             :
         287           0 :        if (trim(ubc_source(m))=='msis') then
         288           0 :           if (do_molec_diff .and. any(msis_flds==ubc_flds(m))) then
         289           0 :              msis_active = .true.
         290           0 :              if (trim(ubc_flds(m))=='H') h_msis_ndx=spc_ndx(m)
         291           0 :              if (trim(ubc_flds(m))=='N') n_msis_ndx=spc_ndx(m)
         292           0 :              if (trim(ubc_flds(m))=='O') o_msis_ndx=spc_ndx(m)
         293           0 :              if (trim(ubc_flds(m))=='O2') o2_msis_ndx=spc_ndx(m)
         294             :           else
         295           0 :              call endrun(prefix//'MSIS is not allowed in this configuration')
         296             :           end if
         297           0 :        else if (trim(ubc_source(m))=='tgcm') then
         298           0 :           if (do_molec_diff .and. any(tgcm_flds==ubc_flds(m))) then
         299           0 :              tgcm_active = .true.
         300             :           else
         301           0 :              call endrun(prefix//'TGCM is not allowed in this configuration')
         302             :           end if
         303           0 :        else if (trim(ubc_source(m))=='snoe') then
         304           0 :           if (do_molec_diff .and. any(snoe_flds==ubc_flds(m))) then
         305           0 :              snoe_active = .true.
         306             :           else
         307           0 :              call endrun(prefix//'SNOE is not allowed in this configuration')
         308             :           end if
         309           0 :        else if (trim(ubc_source(m))=='ubc_file') then
         310           0 :           file_spc_ndx(mm) = spc_ndx(m)
         311           0 :           mm = mm+1
         312             :        else
         313           0 :           mmrndx = index(trim(ubc_source(m)),'mmr')
         314           0 :           vmrndx = index(trim(ubc_source(m)),'vmr')
         315           0 :           if (mmrndx>0 .and. vmrndx>0) then
         316           0 :              call endrun(prefix//'incorrect units in UBC source: '//trim(ubc_source(m)))
         317             :           end if
         318           0 :           if (mmrndx>0) then
         319           0 :              str = ubc_source(m)(:mmrndx-1)
         320           0 :              read(str,*) val
         321           0 :              m_mmr = m_mmr + 1
         322           0 :              fixed_mmr(m_mmr) = val
         323           0 :              fixed_mmr_ndx(m_mmr) = spc_ndx(m)
         324           0 :           else if (vmrndx>0) then
         325           0 :              str = ubc_source(m)(:vmrndx-1)
         326           0 :              read(str,*) val
         327           0 :              m_vmr = m_vmr + 1
         328           0 :              fixed_vmr(m_vmr) = val
         329           0 :              fixed_vmr_ndx(m_vmr) = spc_ndx(m)
         330             :           else
         331           0 :              call endrun(prefix//'unrecognized UBC source: '//trim(ubc_source(m)))
         332             :           end if
         333             :        end if
         334             :     end do
         335             :
         336           0 :     if (tgcm_active) then
         337             :        !-----------------------------------------------------------------------
         338             :        !       ... initialize the tgcm upper boundary module
         339             :        !-----------------------------------------------------------------------
         340             :        call tgcm_ubc_inti( tgcm_ubc_file, tgcm_ubc_data_type, tgcm_ubc_cycle_yr, &
         341           0 :             tgcm_ubc_fixed_ymd, tgcm_ubc_fixed_tod)
         342           0 :        if (masterproc) write(iulog,*) 'ubc_init: after tgcm_ubc_inti'
         343             :     endif
         344             :
         345           0 :     if (snoe_active) then
         346             :        !-----------------------------------------------------------------------
         347             :        !       ... initialize the snoe module
         348             :        !-----------------------------------------------------------------------
         349           0 :        call snoe_inti(snoe_ubc_file)
         350           0 :        if (masterproc) write(iulog,*) 'ubc_init: after snoe_inti'
         351             :     endif
         352             :
         353           0 :     if (msis_active) then
         354             :        !-----------------------------------------------------------------------
         355             :        !       ... initialize the msis module
         356             :        !-----------------------------------------------------------------------
         357           0 :        call msis_ubc_inti( zonal_avg, n_msis_ndx,h_msis_ndx,o_msis_ndx,o2_msis_ndx )
         358           0 :        if (masterproc) write(iulog,*) 'ubc_init: after msis_ubc_inti'
         359             :     endif
         360             :
         361        1536 :   end subroutine ubc_init
         362             :
         363             : !===============================================================================
         364             : !===============================================================================
         365             :
         366           0 :   pure logical function ubc_fixed_conc(name)
         367             :
         368             :     character(len=*), intent(in) :: name
         369             :
         370             :     integer :: m
         371             :
         372           0 :     ubc_fixed_conc = .false.
         373             :
         374           0 :     do m = 1,num_fixed
         375           0 :        if ( trim(ubc_flds(m)) == trim(name) ) then
         376           0 :           ubc_fixed_conc = .true.
         377           0 :           return
         378             :        endif
         379             :     end do
         380             :
         381        1536 :   end function ubc_fixed_conc
         382             :
         383             : !===============================================================================
         384             :
         385      741888 :   subroutine ubc_timestep_init(pbuf2d, state)
         386             : !-----------------------------------------------------------------------
         387             : ! timestep dependent setting
         388             : !-----------------------------------------------------------------------
         389             :
         390             :     use solar_parms_data, only: kp=>solar_parms_kp, ap=>solar_parms_ap, f107=>solar_parms_f107
         391             :     use solar_parms_data, only: f107a=>solar_parms_f107a, f107p=>solar_parms_f107p
         392             :     use mo_msis_ubc,      only: msis_timestep_init
         393             :     use mo_tgcm_ubc,      only: tgcm_timestep_init
         394             :     use mo_snoe,          only: snoe_timestep_init
         395             :     use physics_types,    only: physics_state
         396             :     use ppgrid,           only: begchunk, endchunk
         397             :     use physics_buffer,   only: physics_buffer_desc
         398             :
         399             :     type(physics_state), intent(in) :: state(begchunk:endchunk)
         400             :     type(physics_buffer_desc), pointer :: pbuf2d(:,:)
         401             :
         402      370944 :     if (.not.apply_upper_bc) return
         403             :
         404           0 :     if (num_infile>0) then
         405           0 :        call upper_bc_file_adv( pbuf2d, state )
         406             :     end if
         407           0 :     if (msis_active) then
         408           0 :        call msis_timestep_init( ap, f107p, f107a )
         409             :     end if
         410           0 :     if (tgcm_active) then
         411           0 :        call tgcm_timestep_init( pbuf2d, state )
         412             :     end if
         413           0 :     if (snoe_active) then
         414           0 :        call snoe_timestep_init( kp, f107 )
         415             :     end if
         416             :
         417      370944 :   end subroutine ubc_timestep_init
         418             :
         419             : !===============================================================================
         420             :
         421     1489176 :   subroutine ubc_get_vals (lchnk, ncol, pint, zi, ubc_temp, ubc_mmr)
         422             :
         423             : !-----------------------------------------------------------------------
         424             : ! interface routine for vertical diffusion and pbl scheme
         425             : !-----------------------------------------------------------------------
         426      370944 :     use mo_msis_ubc,      only: get_msis_ubc
         427             :     use mo_snoe,          only: set_no_ubc, ndx_no
         428             :     use mo_tgcm_ubc,      only: set_tgcm_ubc
         429             :     use cam_abortutils,   only: endrun
         430             :     use air_composition,  only: rairv, mbarv ! gas constant, mean mass
         431             :     use constituents,     only: cnst_mw  ! Needed for ubc_flux
         432             :
         433             : !------------------------------Arguments--------------------------------
         434             :     integer,  intent(in)  :: lchnk                 ! chunk identifier
         435             :     integer,  intent(in)  :: ncol                  ! number of atmospheric columns
         436             :     real(r8), intent(in)  :: pint(pcols,pverp)     ! interface pressures
         437             :     real(r8), intent(in)  :: zi(pcols,pverp)       ! interface geoptl height above sfc
         438             :
         439             :     real(r8), intent(out) :: ubc_temp(pcols)      ! upper bndy temperature (K)
         440             :     real(r8), intent(out) :: ubc_mmr(pcols,pcnst)  ! upper bndy mixing ratios (kg/kg)
         441             :
         442             :     !---------------------------Local storage-------------------------------
         443             :     integer :: m                                   ! constituent index
         444             :     real(r8) :: rho_top(pcols)                     ! density at top interface
         445             :     real(r8) :: z_top(pcols)                       ! height of top interface (km)
         446     2978352 :     real(r8) :: vals(pcols,num_infile)
         447             :
         448             :     real(r8), parameter :: m2km = 1.e-3_r8         ! meter to km
         449             :
         450             :     !-----------------------------------------------------------------------
         451             :
         452     1489176 :     ubc_mmr(:,:) = 0._r8
         453     1489176 :     ubc_temp(:) = nan
         454             :
         455     1489176 :     if (.not. apply_upper_bc) return
         456             :
         457             :     ! UBC_FILE
         458           0 :     if (num_infile>0) then
         459           0 :        call upper_bc_file_get(lchnk, ncol, vals)
         460           0 :        do m = 1,num_infile
         461           0 :           if (file_spc_ndx(m)>0) then
         462           0 :              ubc_mmr(:ncol,file_spc_ndx(m)) = vals(:ncol,m)
         463             :           else
         464           0 :              ubc_temp(:ncol) = vals(:ncol,m)
         465             :           end if
         466             :
         467             :        end do
         468             :
         469             :     endif
         470             :
         471             :     ! MSIS
         472           0 :     if (msis_active) then
         473           0 :        call get_msis_ubc( lchnk, ncol, ubc_temp, ubc_mmr )
         474           0 :        if( t_pert_ubc /= 0._r8 ) then
         475           0 :           ubc_temp(:ncol) = ubc_temp(:ncol) + t_pert_ubc
         476           0 :           if( any( ubc_temp(:ncol) < 0._r8 ) ) then
         477           0 :              write(iulog,*) 'ubc_get_vals: msis temp < 0 after applying offset = ',t_pert_ubc
         478           0 :              call endrun('ubc_get_vals: msis temp < 0 after applying t_pert_ubc')
         479             :           end if
         480             :        end if
         481             :     end if
         482             :
         483             :     ! SNOE
         484           0 :     if (snoe_active) then
         485             :
         486           0 :        rho_top(:ncol) = pint(:ncol,1) / (rairv(:ncol,1,lchnk)*ubc_temp(:ncol))
         487           0 :        z_top(:ncol)   = m2km * zi(:ncol,1)
         488             :
         489           0 :        call set_no_ubc  ( lchnk, ncol, z_top, ubc_mmr, rho_top )
         490           0 :        if( ndx_no > 0 .and. no_xfac_ubc /= 1._r8 ) then
         491           0 :           ubc_mmr(:ncol,ndx_no) = no_xfac_ubc * ubc_mmr(:ncol,ndx_no)
         492             :        end if
         493             :
         494             :     endif
         495             :
         496             :     ! TIE-GCM
         497           0 :     if (tgcm_active) then
         498           0 :        call set_tgcm_ubc( lchnk, ncol, ubc_mmr )
         499             :     endif
         500             :
         501             :     ! fixed values
         502           0 :     do m = 1,n_fixed_mmr
         503           0 :        ubc_mmr(:ncol,fixed_mmr_ndx(m)) = fixed_mmr(m)
         504             :     end do
         505           0 :     do m = 1,n_fixed_vmr
         506           0 :        ubc_mmr(:ncol,fixed_vmr_ndx(m)) = cnst_mw(fixed_vmr_ndx(m))*fixed_vmr(m)/mbarv(:ncol,1,lchnk)
         507             :     end do
         508             :
         509             :     ! diagnostic output
         510           0 :     do m = 1,num_fixed
         511           0 :        if (ubc_flds(m)=='T') then
         512           0 :           call outfld(hist_names(m),ubc_temp(:ncol),ncol,lchnk)
         513             :        else
         514           0 :           call outfld(hist_names(m),ubc_mmr(:ncol,spc_ndx(m)),ncol,lchnk)
         515             :        end if
         516             :     end do
         517             :
         518     1489176 :   end subroutine ubc_get_vals
         519             :
         520             : !===============================================================================
         521             :
         522           0 :   subroutine ubc_get_flxs (lchnk, ncol, pint, zi, t, q, phis, ubc_flux)
         523             :
         524     1489176 :     use physconst,       only: avogad, rga
         525             :     use air_composition, only: rairv
         526             :     use constituents,    only: cnst_mw
         527             : !------------------------------Arguments--------------------------------
         528             :     integer,  intent(in)  :: lchnk                 ! chunk identifier
         529             :     integer,  intent(in)  :: ncol                  ! number of atmospheric columns
         530             :     real(r8), intent(in)  :: pint(pcols,pverp)     ! interface pressures
         531             :     real(r8), intent(in)  :: zi(pcols,pverp)       ! interface geoptl height above sfc
         532             :     real(r8), intent(in)  :: t(pcols,pver)         ! midpoint temperature
         533             :     real(r8), intent(in),target :: q(pcols,pver,pcnst)   ! contituent mixing ratios (kg/kg)
         534             :     real(r8), intent(in)  :: phis(pcols)           ! Surface geopotential (m2/s2)
         535             :
         536             :     real(r8), intent(out) :: ubc_flux(pcols,pcnst) ! upper bndy flux (kg/s/m^2)
         537             :
         538             : !---------------------------Local storage-------------------------------
         539             :     integer :: iCol                                ! column loop counter
         540             :
         541             :     real(r8), parameter :: h_escape_flx_factor = 2.03e-13_r8 ! for hydrogen escape flux due to charge exchange
         542             :     ! adopted from TIME-GCM (R. G. Roble, pp. 1-21, AGU Geophys. Monogr. Ser 87, 1995) following
         543             :     ! Liu, S.C., and T. M. Donahue, Mesospheric hydrogen related to exospheric escape mechanisms, J. Atmos. Sci.,
         544             :     ! 31, 1466-1470, 1974. (Equation 4 there). DOI: 10.1175/1520-0469(1974)031<1466:Mhrtee>2.0.Co;2
         545             :     ! https://journals.ametsoc.org/view/journals/atsc/31/5/1520-0469_1974_031_1466_mhrtee_2_0_co_2.xml
         546             :
         547             :     real(r8), parameter :: hfluxlimitfac = 0.72_r8 ! Hydrogen upper boundary flux limiting factor
         548             :
         549             :     real(r8) :: nmbartop                           ! Top level density (rho)
         550             :     real(r8) :: zkt                                ! Factor for H Jean's escape flux calculation
         551             :
         552           0 :     real(r8), pointer :: qh_top(:)                 ! Top level hydrogen mixing ratio (kg/kg)
         553             :
         554           0 :     ubc_flux(:,:) = nan
         555             :
         556           0 :     qh_top => q(:,1,h_ndx)
         557             :
         558           0 :     do iCol = 1, ncol
         559             :        !--------------------------------------------------
         560             :        ! Get total density (rho) at top level
         561             :        !--------------------------------------------------
         562           0 :        nmbartop = 0.5_r8 * (pint(iCol,1) + pint(iCol,2)) / ( rairv(iCol,1,lchnk) * t(iCol,1) )
         563             :
         564             :        !---------------------------------------------------------------------
         565             :        ! Calculate factor for Jean's escape flux once here, used twice below
         566             :        !---------------------------------------------------------------------
         567             :        zkt = (rEarth + ( 0.5_r8 * ( zi(iCol,1) + zi(iCol,2) ) + rga * phis(iCol) ) ) * &
         568           0 :             cnst_mw(h_ndx) / avogad * grav / ( kboltz * t(iCol,1) )
         569             :
         570             :        ubc_flux(iCol,h_ndx) = hfluxlimitfac * SQRT(kboltz/(2.0_r8 * pi * cnst_mw(h_ndx) / avogad)) * &
         571           0 :             qh_top(iCol) * nmbartop * &
         572           0 :             SQRT(t(iCol,1)) * (1._r8 + zkt) * EXP(-zkt)
         573             :
         574             :        ubc_flux(iCol,h_ndx) = ubc_flux(iCol,h_ndx) * &
         575           0 :             (h_escape_flx_factor * qh_top(iCol) * nmbartop / (cnst_mw(h_ndx) / avogad) * t(iCol,1))
         576             :
         577             :     enddo
         578             :
         579           0 :   end subroutine ubc_get_flxs
         580             :
         581             : end module upper_bc

Generated by: LCOV version 1.14