LCOV - code coverage report
Current view: | top level - chemistry/mozart - mo_fstrat.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 34 | 389 | 8.7 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 2 | 4 | 50.0 %

              Line data    Source code
           1             : module mo_fstrat
           2             :   !---------------------------------------------------------------
           3             :   !     ... variables for the upper boundary values
           4             :   !---------------------------------------------------------------
           5             :
           6             :   use shr_kind_mod,     only : r8 => shr_kind_r8
           7             :   use ppgrid,           only : pcols, pver,pverp, begchunk, endchunk
           8             :   use chem_mods,        only : gas_pcnst
           9             :   use cam_abortutils,   only : endrun
          10             :   use cam_pio_utils,    only : cam_pio_openfile
          11             :   use pio
          12             :   use cam_logfile,      only : iulog
          13             :
          14             :   implicit none
          15             :
          16             :   private
          17             :   public  :: fstrat_inti
          18             :   public  :: set_fstrat_vals
          19             :   public  :: set_fstrat_h2o
          20             :   public  :: has_fstrat
          21             :
          22             :   save
          23             :
          24             :   real(r8), parameter :: taurelax = 864000._r8        ! 10 days
          25             :   integer  :: table_nox_ndx = -1
          26             :   integer  :: table_h2o_ndx = -1
          27             :   integer  :: table_ox_ndx  = -1
          28             :   integer  :: no_ndx
          29             :   integer  :: no2_ndx
          30             :   integer  :: h2o_ndx
          31             :   integer  :: ox_ndx
          32             :   integer  :: o3s_ndx
          33             :   integer  :: o3inert_ndx
          34             :   integer  :: o3a_ndx
          35             :   integer  :: xno_ndx
          36             :   integer  :: xno2_ndx
          37             :   integer  :: o3rad_ndx
          38             :   real(r8) :: facrelax
          39             :   real(r8) :: days(12)
          40             :   real(r8), allocatable   :: ub_plevs(:)         ! table midpoint pressure levels (Pa)
          41             :   real(r8), allocatable   :: ub_plevse(:)        ! table edge pressure levels (Pa)
          42             :   integer                 :: ub_nlevs            ! # of levs in ubc file
          43             :   integer                 :: ub_nlat             ! # of lats in ubc file
          44             :   integer                 :: ub_nspecies         ! # of species in ubc file
          45             :   integer                 :: ub_nmonth           ! # of months in ubc file
          46             :   real(r8), allocatable   :: mr_ub(:,:,:,:,:)      ! vmr
          47             :   integer,  allocatable   :: map(:)              ! species indices for ubc species
          48             :   logical :: sim_has_nox
          49             :   integer :: dtime               ! model time step (s)
          50             :   logical :: has_fstrat(gas_pcnst)
          51             :
          52             : contains
          53             :
          54        1536 :   subroutine fstrat_inti( fstrat_file, fstrat_list )
          55             :     !------------------------------------------------------------------
          56             :     !   ... initialize upper boundary values
          57             :     !------------------------------------------------------------------
          58             :
          59             :     use mo_constants,  only : d2r
          60             :     use phys_grid,     only : get_ncols_p, get_rlat_all_p
          61             :
          62             :     use time_manager,  only : get_step_size
          63             :     use time_manager,  only : get_calday
          64             :     use ioFileMod,     only : getfil
          65             :     use spmd_utils,    only : masterproc
          66             :     use mo_tracname,   only : solsym
          67             :     use chem_mods,     only : gas_pcnst
          68             :     use mo_chem_utls,  only : get_spc_ndx, get_inv_ndx
          69             :     use constituents,  only : pcnst
          70             :     use interpolate_data, only : interp_type, lininterp_init, lininterp
          71             :
          72             :     character(len=*), intent(in) :: fstrat_file
          73             :     character(len=*), intent(in) :: fstrat_list(:)
          74             :
          75             :     !------------------------------------------------------------------
          76             :     !   ... local variables
          77             :     !------------------------------------------------------------------
          78             :     real(r8), parameter :: mb2pa = 100._r8
          79             :
          80             :     integer :: i, j, nchar
          81             :     integer :: spcno, lev, month, ierr
          82             :     integer :: vid, ndims
          83             :     type(file_desc_t) :: ncid
          84             :     integer :: dimid_lat, dimid_lev, dimid_species, dimid_month
          85             :     integer :: dimid(4)
          86             :     integer :: start(4)
          87             :     integer :: count(4)
          88             :     integer, parameter :: dates(12) = (/ 116, 214, 316, 415,  516,  615, &
          89             :          716, 816, 915, 1016, 1115, 1216 /)
          90             :     integer :: ncols, c
          91        1536 :     real(r8), allocatable :: mr_ub_in(:,:,:,:)
          92        1536 :     real(r8), allocatable :: lat(:)
          93             :     real(r8) :: to_lats(pcols)
          94             :     character(len=80) :: attribute
          95             :     character(len=8)  :: wrk_name
          96        1536 :     character(len=25), allocatable :: ub_species_names(:)
          97             :     character(len=256) :: locfn
          98             :     type(interp_type) :: lat_wgts
          99             :
         100             :
         101             :     !-----------------------------------------------------------------------
         102             :     !       ... get species indicies
         103             :     !-----------------------------------------------------------------------
         104        1536 :     no_ndx      = get_spc_ndx( 'NO' )
         105        1536 :     no2_ndx     = get_spc_ndx( 'NO2' )
         106        1536 :     sim_has_nox = no_ndx > 0 .or. no2_ndx > 0
         107        1536 :     ox_ndx      = get_spc_ndx( 'OX' )
         108        1536 :     if( ox_ndx < 1 ) then
         109        1536 :        ox_ndx = get_spc_ndx( 'O3' )
         110             :     end if
         111        1536 :     o3s_ndx     = get_spc_ndx( 'O3S' )
         112        1536 :     o3inert_ndx = get_spc_ndx( 'O3INERT' )
         113        1536 :     o3rad_ndx   = get_spc_ndx( 'O3RAD' )
         114        1536 :     o3a_ndx     = get_spc_ndx( 'O3A' )
         115        1536 :     xno_ndx      = get_spc_ndx( 'XNO' )
         116        1536 :     xno2_ndx      = get_spc_ndx( 'XNO2' )
         117             :
         118        1536 :     if(masterproc) then
         119           2 :        if (ox_ndx > 0) then
         120           0 :           if ( .not. any(fstrat_list == 'O3') ) then
         121           0 :              write(iulog,*) 'fstrat_inti: ***WARNING*** O3 is not include in fstrat_list namelist variable'
         122             :           endif
         123             :        endif
         124             :     end if
         125        1536 :     h2o_ndx     = get_spc_ndx( 'H2O' )
         126        1536 :     if( h2o_ndx < 0 ) then
         127           0 :        h2o_ndx  = get_inv_ndx( 'H2O' )
         128             :     end if
         129             :
         130        1536 :     has_fstrat(:) = .false.
         131             :
         132        1536 :     do i = 1,pcnst
         133             :
         134        1536 :        if ( len_trim(fstrat_list(i))==0 ) exit
         135             :
         136           0 :        j = get_spc_ndx(fstrat_list(i))
         137             :
         138        1536 :        if ( j > 0 ) then
         139           0 :           has_fstrat(j) = .true.
         140             :        else
         141           0 :           write(iulog,*) 'fstrat_inti: '//trim(fstrat_list(i))//' is not included in species set'
         142           0 :           call endrun('fstrat_inti: invalid fixed stratosphere species')
         143             :        endif
         144             :
         145             :     enddo
         146             :
         147       49152 :     if (.not. any(has_fstrat(:))) return
         148             :
         149             :     !-----------------------------------------------------------------------
         150             :     !       ... open netcdf file
         151             :     !-----------------------------------------------------------------------
         152           0 :     call getfil (fstrat_file, locfn, 0)
         153           0 :     call cam_pio_openfile (ncid, trim(locfn), PIO_NOWRITE)
         154             :     !-----------------------------------------------------------------------
         155             :     !       ... get latitude
         156             :     !-----------------------------------------------------------------------
         157           0 :     ierr = pio_inq_dimid( ncid, 'lat', dimid_lat )
         158           0 :     ierr = pio_inq_dimlen( ncid, dimid_lat, ub_nlat )
         159           0 :     allocate( lat(ub_nlat), stat=ierr )
         160           0 :     if( ierr /= 0 ) then
         161           0 :        write(iulog,*) 'fstrat_inti: lat allocation error = ',ierr
         162           0 :        call endrun
         163             :     end if
         164           0 :     ierr = pio_inq_varid( ncid, 'lat', vid )
         165           0 :     ierr = pio_get_var( ncid, vid, lat )
         166           0 :     lat(:ub_nlat) = lat(:ub_nlat) * d2r
         167             :
         168           0 :     if( ierr /= 0 ) then
         169           0 :        write(iulog,*) 'fstrat_inti: failed to deallocate lat; ierr = ',ierr
         170           0 :        call endrun
         171             :     end if
         172             :
         173             :     !-----------------------------------------------------------------------
         174             :     !       ... get vertical coordinate (if necessary, convert units to pa)
         175             :     !-----------------------------------------------------------------------
         176           0 :     ierr = pio_inq_dimid( ncid, 'lev', dimid_lev )
         177           0 :     ierr = pio_inq_dimlen( ncid, dimid_lev, ub_nlevs )
         178           0 :     allocate( ub_plevs(ub_nlevs), stat=ierr )
         179           0 :     if( ierr /= 0 ) then
         180           0 :        write(iulog,*) 'fstrat_inti: ub_plevs allocation error = ',ierr
         181           0 :        call endrun
         182             :     end if
         183           0 :     ierr = pio_inq_varid( ncid, 'lev', vid )
         184           0 :     ierr = pio_get_var( ncid, vid, ub_plevs )
         185           0 :     attribute(:) = ' '
         186           0 :     call pio_seterrorhandling(ncid, pio_BCAST_error)
         187           0 :     ierr = pio_get_att( ncid, vid, 'units', attribute )
         188           0 :     call pio_seterrorhandling(ncid, pio_INTERNAL_error)
         189           0 :     if( ierr == PIO_noerr )then
         190           0 :        if( trim(attribute) == 'mb' .or. trim(attribute) == 'hpa' )then
         191           0 :           if (masterproc) write(iulog,*) 'fstrat_inti: units for lev = ',trim(attribute),'... converting to pa'
         192           0 :           ub_plevs(:) = mb2pa * ub_plevs(:)
         193           0 :        else if( trim(attribute) /= 'pa' .and. trim(attribute) /= 'pa' )then
         194           0 :           write(iulog,*) 'fstrat_inti: unknown units for lev, units=*',trim(attribute),'*'
         195           0 :           write(iulog,*) 'fstrat_inti: ',attribute=='mb',trim(attribute)=='mb',attribute(1:2)=='mb'
         196           0 :           call endrun
         197             :        end if
         198             :     else
         199           0 :        write(iulog,*) 'fstrat_inti: warning! units attribute for lev missing, assuming mb'
         200           0 :        ub_plevs(:) = mb2pa * ub_plevs(:)
         201             :     end if
         202             :     !-----------------------------------------------------------------------
         203             :     !       ... get time and species dimensions
         204             :     !-----------------------------------------------------------------------
         205           0 :     ierr = pio_inq_dimid( ncid, 'month', dimid_month )
         206           0 :     ierr = pio_inq_dimlen( ncid, dimid_month, ub_nmonth )
         207           0 :     if( ub_nmonth /= 12 )then
         208           0 :        write(iulog,*) 'fstrat_inti: error! number of months = ',ub_nmonth,', expecting 12'
         209           0 :        call endrun
         210             :     end if
         211           0 :     ierr = pio_inq_dimid( ncid, 'species', dimid_species )
         212           0 :     ierr = pio_inq_dimlen( ncid, dimid_species, ub_nspecies )
         213             :
         214             :     !------------------------------------------------------------------
         215             :     !   ... allocate arrays
         216             :     !------------------------------------------------------------------
         217           0 :     allocate( mr_ub_in(ub_nlat,ub_nspecies,ub_nmonth,ub_nlevs), stat=ierr )
         218           0 :     if( ierr /= 0 ) then
         219           0 :        write(iulog,*) 'fstrat_inti: mr_ub_in allocation error = ',ierr
         220           0 :        call endrun
         221             :     end if
         222           0 :     allocate( mr_ub(pcols,ub_nspecies,ub_nmonth,ub_nlevs,begchunk:endchunk), stat=ierr )
         223           0 :     if( ierr /= 0 ) then
         224           0 :        write(iulog,*) 'fstrat_inti: mr_ub allocation error = ',ierr
         225           0 :        call endrun
         226             :     end if
         227             :
         228             :     !------------------------------------------------------------------
         229             :     !   ... read in the species names
         230             :     !------------------------------------------------------------------
         231             :
         232           0 :     ierr = pio_inq_varid( ncid, 'specname', vid )
         233           0 :     ierr = pio_inq_vardimid( ncid, vid, dimid )
         234           0 :     ierr = pio_inq_dimlen( ncid, dimid(1), nchar )
         235           0 :     allocate( ub_species_names(ub_nspecies), stat=ierr )
         236           0 :     if( ierr /= 0 ) then
         237           0 :        write(iulog,*) 'fstrat_inti: ub_species_names allocation error = ',ierr
         238           0 :        call endrun
         239             :     end if
         240           0 :     allocate( map(ub_nspecies), stat=ierr )
         241           0 :     if( ierr /= 0 ) then
         242           0 :        write(iulog,*) 'fstrat_inti: map allocation error = ',ierr
         243           0 :        call endrun
         244             :     end if
         245             :
         246           0 :     table_loop :  do i = 1,ub_nspecies
         247           0 :        start(:2) = (/ 1, i /)
         248           0 :        count(:2) = (/ nchar, 1 /)
         249           0 :        ub_species_names(i)(:) = ' '
         250           0 :        ierr = pio_get_var( ncid, vid, start(:2), count(:2), ub_species_names(i:i))
         251           0 :        if( trim(ub_species_names(i)) == 'NOX' ) then
         252           0 :           table_nox_ndx = i
         253           0 :        else if( trim(ub_species_names(i)) == 'H2O' ) then
         254           0 :           table_h2o_ndx = i
         255           0 :        else if( trim(ub_species_names(i)) == 'OX' ) then
         256           0 :           table_ox_ndx = i
         257             :        end if
         258           0 :        map(i) = 0
         259           0 :        do j = 1,gas_pcnst
         260           0 :           if( trim(ub_species_names(i)) == trim(solsym(j)) ) then
         261           0 :              if( has_fstrat(j) ) then
         262           0 :                 map(i) = j
         263           0 :                 if( masterproc ) write(iulog,*) 'fstrat_inti: '//trim(solsym(j))//' is fixed in stratosphere'
         264           0 :                 exit
         265             :              end if
         266             :           endif
         267             :        end do
         268           0 :        if( map(i) == 0 ) then
         269           0 :           if( trim(ub_species_names(i)) == 'OX' ) then
         270           0 :              if( o3rad_ndx > 0 ) then
         271           0 :                 wrk_name = 'O3RAD'
         272             :              else
         273           0 :                 wrk_name = 'O3'
         274             :              end if
         275           0 :              do j = 1,gas_pcnst
         276           0 :                 if( trim(wrk_name) == trim(solsym(j)) ) then
         277           0 :                    if( has_fstrat(j) ) then
         278           0 :                       if( masterproc ) write(iulog,*) 'fstrat_inti: '//trim(solsym(j))//' is fixed in stratosphere'
         279           0 :                       map(i) = j
         280           0 :                       exit
         281             :                    end if
         282             :                 end if
         283             :              end do
         284           0 :              if( map(i) == 0 ) then
         285           0 :                 write(iulog,*) 'fstrat_inti: ubc table species ',trim(ub_species_names(i)), ' not used'
         286             :              end if
         287           0 :           else if( (trim(ub_species_names(i)) /= 'NOX') .and. (trim(ub_species_names(i)) /= 'H2O') ) then
         288           0 :              write(iulog,*) 'fstrat_inti: ubc table species ',trim(ub_species_names(i)), ' not used'
         289             :           end if
         290             :        end if
         291             :     end do table_loop
         292             :
         293           0 :     if( table_nox_ndx > 0 ) then
         294           0 :        if ( any(fstrat_list(:) == 'NO') .or. any(fstrat_list(:) == 'NO2') ) then
         295           0 :           map(table_nox_ndx) = gas_pcnst + 1
         296             :        else
         297           0 :           write(iulog,*) 'fstrat_inti: ubc table species ',trim(ub_species_names(table_nox_ndx)), ' not used'
         298             :        endif
         299             :     end if
         300           0 :     if( table_h2o_ndx > 0 ) then
         301           0 :        if ( h2o_ndx > 0 ) then
         302           0 :           map(table_h2o_ndx) = gas_pcnst + 2
         303             :        else
         304           0 :           write(iulog,*) 'fstrat_inti: ubc table species ',trim(ub_species_names(table_h2o_ndx)), ' not used'
         305             :        endif
         306             :     end if
         307             :
         308           0 :     if (masterproc) write(iulog,*) 'fstrat_inti: h2o_ndx, table_h2o_ndx  = ', h2o_ndx, table_h2o_ndx
         309             :
         310             :     !------------------------------------------------------------------
         311             :     !   ... check dimensions for vmr variable
         312             :     !------------------------------------------------------------------
         313           0 :     ierr = pio_inq_varid( ncid, 'vmr', vid )
         314           0 :     ierr = pio_inq_varndims( ncid, vid, ndims )
         315           0 :     if( ndims /= 4 ) then
         316           0 :        write(iulog,*) 'fstrat_inti: error! variable vmr has ndims = ',ndims,', expecting 4'
         317           0 :        call endrun
         318             :     end if
         319           0 :     ierr = pio_inq_vardimid( ncid, vid, dimid )
         320             :     if( dimid(1) /= dimid_lat .or. dimid(2) /= dimid_species .or. &
         321           0 :          dimid(3) /= dimid_month .or. dimid(4) /= dimid_lev )then
         322             :        write(iulog,*) 'fstrat_inti: error! dimensions in wrong order for variable vmr,'// &
         323           0 :             'expecting (lat,species,month,lev)'
         324           0 :        call endrun
         325             :     end if
         326             :
         327             :     !------------------------------------------------------------------
         328             :     !   ... read in the ub mixing ratio values
         329             :     !------------------------------------------------------------------
         330           0 :     start = (/ 1, 1, 1, 1 /)
         331           0 :     count = (/ ub_nlat, ub_nspecies, ub_nmonth, ub_nlevs /)
         332             :
         333           0 :     ierr = pio_get_var( ncid, vid, start, count, mr_ub_in )
         334           0 :     call pio_closefile (ncid)
         335             :     !--------------------------------------------------------------------
         336             :     !   ... regrid
         337             :     !--------------------------------------------------------------------
         338           0 :     do c=begchunk,endchunk
         339           0 :        ncols = get_ncols_p(c)
         340           0 :        call get_rlat_all_p(c, pcols, to_lats)
         341           0 :        call lininterp_init(lat,ub_nlat, to_lats, ncols, 1, lat_wgts)
         342             :
         343           0 :        do lev = 1,ub_nlevs
         344           0 :           do month = 1,ub_nmonth
         345           0 :              do spcno = 1,ub_nspecies
         346           0 :                 if( map(spcno) > 0 ) then
         347             :
         348             :                    call lininterp(mr_ub_in(:,spcno, month, lev), ub_nlat, mr_ub(:,spcno,month,lev,c), &
         349           0 :                         ncols, lat_wgts)
         350             : #ifdef DEBUG
         351           0 :                    if( lev == 25 .and. month == 1 .and. spcno == 1 ) then
         352           0 :                       write(iulog,*) 'mr_ub_in='
         353           0 :                       write(iulog,'(10f7.1)') mr_ub_in(:,spcno,month,lev)*1.e9_r8
         354           0 :                       write(iulog,*) 'mr_ub='
         355           0 :                       write(iulog,'(10f7.1)') mr_ub(:,spcno,month,lev,c)*1.e9_r8
         356             :                    end if
         357             : #endif
         358             :                 end if
         359             :
         360             :              end do
         361             :           end do
         362             :        end do
         363             :     end do
         364             :     !--------------------------------------------------------
         365             :     !   ... initialize the monthly day of year times
         366             :     !--------------------------------------------------------
         367           0 :     do month = 1,12
         368           0 :        days(month) = get_calday( dates(month), 0 )
         369             :     end do
         370             :
         371             :     !--------------------------------------------------------
         372             :     !           ... set up the relaxation for lower stratosphere
         373             :     !--------------------------------------------------------
         374             :     !   ... taurelax = relaxation timescale (in sec)
         375             :     !           facrelax = fractional relaxation towards ubc
         376             :     !            1 => use ubc
         377             :     !            0 => ignore ubc, use model concentrations
         378             :     !--------------------------------------------------------
         379           0 :     dtime = get_step_size()
         380           0 :     facrelax = 1._r8 - exp( -real(dtime)/taurelax )
         381             :
         382             :     !--------------------------------------------------------
         383             :     !   ... setup conserving interp for OX
         384             :     !--------------------------------------------------------
         385           0 :     if( table_ox_ndx > 0 ) then
         386           0 :        allocate( ub_plevse(ub_nlevs-1), stat=ierr )
         387           0 :        if( ierr /= 0 ) then
         388           0 :           write(iulog,*) 'fstrat_inti: ub_plevse allocation error = ',ierr
         389           0 :           call endrun
         390             :        end if
         391           0 :        ub_plevse(1:ub_nlevs-1) = .5_r8*(ub_plevs(1:ub_nlevs-1) + ub_plevs(2:ub_nlevs))
         392             :     end if
         393             :
         394        3072 :   end subroutine fstrat_inti
         395             :
         396     1489176 :   subroutine set_fstrat_vals( vmr, pmid, pint, ltrop, calday, ncol,lchnk )
         397             :
         398             :     !--------------------------------------------------------------------
         399             :     !   ... set the upper boundary values for :
         400             :     !           ox, nox, hno3, ch4, co, n2o, n2o5 & stratospheric o3
         401             :     !--------------------------------------------------------------------
         402             :
         403             :     !--------------------------------------------------------------------
         404             :     !   ... dummy args
         405             :     !--------------------------------------------------------------------
         406             :     integer,  intent(in)    :: lchnk             ! chunk number
         407             :     integer,  intent(in)    :: ncol              ! columns in chunk
         408             :     integer,  intent(in)    :: ltrop(pcols)      ! tropopause vertical index
         409             :     real(r8), intent(in)    :: calday            ! day of year including fraction
         410             :     real(r8), intent(in)    :: pmid(pcols,pver)  ! midpoint pressure (Pa)
         411             :     real(r8), intent(in)    :: pint(pcols,pverp) ! interface pressure (Pa)
         412             :     real(r8), intent(inout) :: vmr(ncol,pver,gas_pcnst) ! species concentrations (mol/mol)
         413             :
         414             :     !--------------------------------------------------------------------
         415             :     !   ... local variables
         416             :     !--------------------------------------------------------------------
         417             :     integer, parameter :: zlower = pver
         418             :
         419             :     integer  ::  m, last, next, i, k, k1, km
         420             :     integer  ::  astat
         421     2978352 :     integer  ::  kmax(ncol)
         422             :     integer  ::  levrelax
         423     2978352 :     integer  ::  kl(ncol,zlower)
         424     2978352 :     integer  ::  ku(ncol,zlower)
         425             :     real(r8) ::  vmrrelax
         426             :     real(r8) ::  pinterp
         427             :     real(r8) ::  nox_ubc, xno, xno2, rno
         428             :     real(r8) ::  dels
         429     1489176 :     real(r8) ::  delp(ncol,zlower)
         430             :     real(r8) ::  pint_vals(2)
         431     1489176 :     real(r8), allocatable :: table_ox(:)
         432             :
         433    47653632 :     if (.not. any(has_fstrat(:))) return
         434             :
         435             :     !--------------------------------------------------------
         436             :     !   ... setup the time interpolation
         437             :     !--------------------------------------------------------
         438           0 :     if( calday < days(1) ) then
         439           0 :        next = 1
         440           0 :        last = 12
         441           0 :        dels = (365._r8 + calday - days(12)) / (365._r8 + days(1) - days(12))
         442           0 :     else if( calday >= days(12) ) then
         443           0 :        next = 1
         444           0 :        last = 12
         445           0 :        dels = (calday - days(12)) / (365._r8 + days(1) - days(12))
         446             :     else
         447           0 :        do m = 11,1,-1
         448           0 :           if( calday >= days(m) ) then
         449             :              exit
         450             :           end if
         451             :        end do
         452           0 :        last = m
         453           0 :        next = m + 1
         454           0 :        dels = (calday - days(m)) / (days(m+1) - days(m))
         455             :     end if
         456           0 :     dels = max( min( 1._r8,dels ),0._r8 )
         457             :
         458             :     !--------------------------------------------------------
         459             :     !   ... setup the pressure interpolation
         460             :     !--------------------------------------------------------
         461           0 :     do k = 1,zlower
         462           0 :        do i = 1,ncol
         463           0 :           if( pmid(i,k) <= ub_plevs(1) ) then
         464           0 :              kl(i,k) = 1
         465           0 :              ku(i,k) = 1
         466           0 :              delp(i,k) = 0._r8
         467           0 :           else if( pmid(i,k) >= ub_plevs(ub_nlevs) ) then
         468           0 :              kl(i,k) = ub_nlevs
         469           0 :              ku(i,k) = ub_nlevs
         470           0 :              delp(i,k) = 0._r8
         471             :           else
         472             :              pinterp = pmid(i,k)
         473           0 :              do k1 = 2,ub_nlevs
         474           0 :                 if( pinterp <= ub_plevs(k1) ) then
         475           0 :                    ku(i,k) = k1
         476           0 :                    kl(i,k) = k1 - 1
         477           0 :                    delp(i,k) = log( pinterp/ub_plevs(kl(i,k)) ) &
         478           0 :                         / log( ub_plevs(ku(i,k))/ub_plevs(kl(i,k)) )
         479           0 :                    exit
         480             :                 end if
         481             :              end do
         482             :           end if
         483             :        end do
         484             :     end do
         485             :
         486             :     !--------------------------------------------------------
         487             :     !   ... find max level less than 50 mb
         488             :     !           fix UB vals from top of model to this level
         489             :     !--------------------------------------------------------
         490           0 :     do i = 1,ncol
         491           0 :        do k = 2,pver
         492           0 :           if( pmid(i,k) > 50.e2_r8 ) then
         493           0 :              kmax(i) = k
         494           0 :              exit
         495             :           end if
         496             :        end do
         497             :     end do
         498             :
         499             :     !--------------------------------------------------------
         500             :     !   ... setup for ox conserving interp
         501             :     !--------------------------------------------------------
         502           0 :     if( table_ox_ndx > 0 ) then
         503           0 :        if( map(table_ox_ndx) > 0 ) then
         504           0 :           allocate( table_ox(ub_nlevs-2),stat=astat )
         505           0 :           if( astat /= 0 ) then
         506           0 :              write(iulog,*) 'set_fstrat_vals: table_ox allocation error = ',astat
         507           0 :              call endrun
         508             :           end if
         509             : #ifdef UB_DEBUG
         510             :           write(iulog,*) ' '
         511             :           write(iulog,*) 'set_fstrat_vals: ub_nlevs = ',ub_nlevs
         512             :           write(iulog,*) 'set_fstrat_vals: ub_plevse'
         513             :           write(iulog,'(1p5g15.7)') ub_plevse(:)
         514             :           write(iulog,*) ' '
         515             : #endif
         516             :        end if
         517             :     end if
         518             :
         519             :     !--------------------------------------------------------
         520             :     !   ... set the mixing ratios at upper boundary
         521             :     !--------------------------------------------------------
         522           0 :     species_loop : do m = 1,ub_nspecies
         523           0 :        ub_overwrite : if( map(m) > 0 ) then
         524           0 :           if( m == table_ox_ndx ) then
         525           0 :              do i = 1,ncol
         526             :
         527           0 :                 table_ox(1:ub_nlevs-2) = mr_ub(i,m,last,2:ub_nlevs-1,lchnk) &
         528           0 :                      + dels*(mr_ub(i,m,next,2:ub_nlevs-1,lchnk) &
         529           0 :                      - mr_ub(i,m,last,2:ub_nlevs-1,lchnk))
         530             : #ifdef UB_DEBUG
         531             :                 write(iulog,*) 'set_fstrat_vals: table_ox @ i = ',i
         532             :                 write(iulog,'(1p5g15.7)') table_ox(:)
         533             :                 write(iulog,*) ' '
         534             : #endif
         535             :
         536           0 :                 km = kmax(i)
         537             : #ifdef UB_DEBUG
         538             :                 write(iulog,*) 'set_fstrat_vals: pint with km = ',km
         539             :                 write(iulog,'(1p5g15.7)') pint(i,:km+1)
         540             :                 write(iulog,*) ' '
         541             :                 write(iulog,*) 'set_fstrat_vals: pmid with km = ',km
         542             :                 write(iulog,'(1p5g15.7)') pmid(i,:km)
         543             :                 write(iulog,*) ' '
         544             : #endif
         545           0 :                 call rebin( ub_nlevs-2, km, ub_plevse, pint(i,:km+1), table_ox, vmr(i,:km,map(m)) )
         546             : #ifdef UB_DEBUG
         547             :                 write(iulog,*) 'set_fstrat_vals: ub o3 @ i = ',i
         548             :                 write(iulog,'(1p5g15.7)') vmr(i,:km,map(m))
         549             : #endif
         550             :              end do
         551           0 :              deallocate( table_ox )
         552           0 :              cycle species_loop
         553             :           end if
         554           0 :           do i = 1,ncol
         555           0 :              do k = 1,kmax(i)
         556           0 :                 pint_vals(1) = mr_ub(i,m,last,kl(i,k),lchnk) &
         557             :                      + delp(i,k) &
         558           0 :                      * (mr_ub(i,m,last,ku(i,k),lchnk) &
         559           0 :                      - mr_ub(i,m,last,kl(i,k),lchnk))
         560           0 :                 pint_vals(2) = mr_ub(i,m,next,kl(i,k),lchnk) &
         561             :                      + delp(i,k) &
         562             :                      * (mr_ub(i,m,next,ku(i,k),lchnk) &
         563           0 :                      - mr_ub(i,m,next,kl(i,k),lchnk))
         564           0 :                 if( m /= table_nox_ndx .and. m /= table_h2o_ndx ) then
         565           0 :                    vmr(i,k,map(m)) = pint_vals(1) &
         566           0 :                         + dels * (pint_vals(2) - pint_vals(1))
         567           0 :                 else if( m == table_nox_ndx .and. sim_has_nox ) then
         568             :
         569           0 :                    nox_ubc = pint_vals(1) + dels * (pint_vals(2) - pint_vals(1))
         570           0 :                    if( no_ndx > 0 ) then
         571           0 :                       xno  = vmr(i,k,no_ndx)
         572             :                    else
         573             :                       xno  = 0._r8
         574             :                    end if
         575           0 :                    if( no2_ndx > 0 ) then
         576           0 :                       xno2 = vmr(i,k,no2_ndx)
         577             :                    else
         578             :                       xno2 = 0._r8
         579             :                    end if
         580           0 :                    rno  = xno / (xno + xno2)
         581           0 :                    if( no_ndx > 0 ) then
         582           0 :                       vmr(i,k,no_ndx)  = rno * nox_ubc
         583             :                    end if
         584           0 :                    if( no2_ndx > 0 ) then
         585           0 :                       vmr(i,k,no2_ndx) = (1._r8 - rno) * nox_ubc
         586             :                    end if
         587             :                 end if
         588             :              end do
         589             :           end do
         590             :        end if ub_overwrite
         591             :     end do species_loop
         592             :
         593           0 :     col_loop2 : do i = 1,ncol
         594             :        !--------------------------------------------------------
         595             :        !        ... relax lower stratosphere to extended ubc
         596             :        !           check to make sure ubc is not being imposed too low
         597             :        !           levrelax = lowest model level (highest pressure)
         598             :        !                      in which to relax to ubc
         599             :        !--------------------------------------------------------
         600           0 :        levrelax = ltrop(i)
         601           0 :        do while( pmid(i,levrelax) > ub_plevs(ub_nlevs) )
         602           0 :           levrelax = levrelax - 1
         603             :        end do
         604             : #ifdef DEBUG
         605           0 :        if( levrelax /= ltrop(i) ) then
         606           0 :           write(iulog,*) 'warning -- raised ubc: ',i,        &
         607           0 :              ltrop(i)-1,nint(pmid(i,ltrop(i)-1)/100._r8),'mb -->', &
         608           0 :              levrelax,nint(pmid(i,levrelax)/100._r8),'mb'
         609             :        end if
         610             : #endif
         611           0 :        level_loop2 : do k = kmax(i)+1,levrelax
         612           0 :           if( sim_has_nox ) then
         613           0 :              if( no_ndx > 0 ) then
         614           0 :                 xno  = vmr(i,k,no_ndx)
         615             :              else
         616             :                 xno  = 0._r8
         617             :              end if
         618           0 :              if( no2_ndx > 0 ) then
         619           0 :                 xno2 = vmr(i,k,no2_ndx)
         620             :              else
         621             :                 xno2 = 0._r8
         622             :              end if
         623           0 :              rno     = xno / (xno + xno2)
         624           0 :              nox_ubc = xno + xno2
         625             :           end if
         626           0 :           do m = 1,ub_nspecies
         627           0 :              if( map(m) < 1 ) then
         628             :                 cycle
         629             :              end if
         630           0 :              pint_vals(1) = mr_ub(i,m,last,kl(i,k),lchnk) &
         631             :                   + delp(i,k) &
         632           0 :                   * (mr_ub(i,m,last,ku(i,k),lchnk) &
         633           0 :                   - mr_ub(i,m,last,kl(i,k),lchnk))
         634           0 :              pint_vals(2) = mr_ub(i,m,next,kl(i,k),lchnk) &
         635             :                   + delp(i,k) &
         636             :                   * (mr_ub(i,m,next,ku(i,k),lchnk) &
         637           0 :                   - mr_ub(i,m,next,kl(i,k),lchnk))
         638             :              vmrrelax = pint_vals(1) &
         639           0 :                   + dels * (pint_vals(2) - pint_vals(1))
         640           0 :              if( m /= table_nox_ndx .and. m /= table_h2o_ndx ) then
         641           0 :                 vmr(i,k,map(m)) = vmr(i,k,map(m)) &
         642           0 :                      + (vmrrelax - vmr(i,k,map(m))) * facrelax
         643           0 :              else if( m == table_nox_ndx .and. sim_has_nox) then
         644             :
         645           0 :                 nox_ubc = nox_ubc + (vmrrelax - nox_ubc) * facrelax
         646             :              end if
         647             :           end do
         648           0 :           if( sim_has_nox ) then
         649           0 :              if( no_ndx > 0 ) then
         650           0 :                 vmr(i,k,no_ndx)  = rno * nox_ubc
         651             :              end if
         652           0 :              if( no2_ndx > 0 ) then
         653           0 :                 vmr(i,k,no2_ndx) = (1._r8 - rno) * nox_ubc
         654             :              end if
         655             :           end if
         656             :        end do level_loop2
         657             :
         658             :        !--------------------------------------------------------
         659             :        !       ... set O3S and O3INERT to OX when no synoz
         660             :        !--------------------------------------------------------
         661           0 :        if( ox_ndx > 0 ) then
         662           0 :           if( o3s_ndx > 0 ) then
         663           0 :              vmr(i,:ltrop(i),o3s_ndx)     = vmr(i,:ltrop(i),ox_ndx)
         664             :           end if
         665           0 :           if( o3inert_ndx > 0 ) then
         666           0 :              vmr(i,:ltrop(i),o3inert_ndx) = vmr(i,:ltrop(i),ox_ndx)
         667             :           end if
         668             :        end if
         669             :
         670             :
         671           0 :        if ( o3a_ndx > 0 ) then
         672           0 :           vmr(i,:ltrop(i),o3a_ndx) = (1._r8 - facrelax ) * vmr(i,:ltrop(i),o3a_ndx)
         673             :        endif
         674           0 :        if ( xno_ndx > 0 ) then
         675           0 :           vmr(i,:ltrop(i),xno_ndx) = (1._r8 - facrelax ) * vmr(i,:ltrop(i),xno_ndx)
         676             :        endif
         677           0 :        if ( xno2_ndx > 0 ) then
         678           0 :           vmr(i,:ltrop(i),xno2_ndx) = (1._r8 - facrelax ) * vmr(i,:ltrop(i),xno2_ndx)
         679             :        endif
         680             :
         681             :     end do col_loop2
         682             :
         683        1536 :   end subroutine set_fstrat_vals
         684             :
         685           0 :   subroutine set_fstrat_h2o( h2o, pmid, ltrop, calday, ncol, lchnk )
         686             :     !--------------------------------------------------------------------
         687             :     !   ... set the h2o upper boundary values
         688             :     !--------------------------------------------------------------------
         689             :
         690             :     !--------------------------------------------------------------------
         691             :     !   ... dummy args
         692             :     !--------------------------------------------------------------------
         693             :     integer,  intent(in)    :: ncol              ! columns in chunk
         694             :     integer,  intent(in)    :: ltrop(pcols)      ! tropopause vertical index
         695             :     integer,  intent(in)    :: lchnk
         696             :     real(r8), intent(in)    :: calday            ! day of year including fraction
         697             :     real(r8), intent(in)    :: pmid(pcols,pver)  ! midpoint pressure (Pa)
         698             :     real(r8), intent(inout) :: h2o(ncol,pver)    ! h2o concentration (mol/mol)
         699             :
         700             :     !--------------------------------------------------------------------
         701             :     !   ... local variables
         702             :     !--------------------------------------------------------------------
         703             :     integer, parameter :: zlower = pver
         704             :
         705             :     integer  ::  m, last, next, i, k, k1
         706           0 :     integer  ::  kmax(ncol)
         707             :     integer  ::  levrelax
         708           0 :     integer  ::  kl(ncol,zlower)
         709           0 :     integer  ::  ku(ncol,zlower)
         710             :     real(r8) ::  vmrrelax
         711             :     real(r8) ::  pinterp
         712             :     real(r8) ::  dels
         713           0 :     real(r8) ::  delp(ncol,zlower)
         714             :     real(r8) ::  pint_vals(2)
         715             :
         716           0 :     h2o_overwrite : if( h2o_ndx > 0 .and. table_h2o_ndx > 0 ) then
         717             :        !--------------------------------------------------------
         718             :        !        ... setup the time interpolation
         719             :        !--------------------------------------------------------
         720           0 :        if( calday < days(1) ) then
         721           0 :           next = 1
         722           0 :           last = 12
         723           0 :           dels = (365._r8 + calday - days(12)) / (365._r8 + days(1) - days(12))
         724           0 :        else if( calday >= days(12) ) then
         725           0 :           next = 1
         726           0 :           last = 12
         727           0 :           dels = (calday - days(12)) / (365._r8 + days(1) - days(12))
         728             :        else
         729           0 :           do m = 11,1,-1
         730           0 :              if( calday >= days(m) ) then
         731             :                 exit
         732             :              end if
         733             :           end do
         734           0 :           last = m
         735           0 :           next = m + 1
         736           0 :           dels = (calday - days(m)) / (days(m+1) - days(m))
         737             :        end if
         738           0 :        dels = max( min( 1._r8,dels ),0._r8 )
         739             :
         740             :        !--------------------------------------------------------
         741             :        !        ... setup the pressure interpolation
         742             :        !--------------------------------------------------------
         743           0 :        do k = 1,zlower
         744           0 :           do i = 1,ncol
         745           0 :              if( pmid(i,k) <= ub_plevs(1) ) then
         746           0 :                 kl(i,k) = 1
         747           0 :                 ku(i,k) = 1
         748           0 :                 delp(i,k) = 0._r8
         749           0 :              else if( pmid(i,k) >= ub_plevs(ub_nlevs) ) then
         750           0 :                 kl(i,k) = ub_nlevs
         751           0 :                 ku(i,k) = ub_nlevs
         752           0 :                 delp(i,k) = 0._r8
         753             :              else
         754             :                 pinterp = pmid(i,k)
         755           0 :                 do k1 = 2,ub_nlevs
         756           0 :                    if( pinterp <= ub_plevs(k1) ) then
         757           0 :                       ku(i,k) = k1
         758           0 :                       kl(i,k) = k1 - 1
         759           0 :                       delp(i,k) = log( pinterp/ub_plevs(kl(i,k)) ) &
         760           0 :                            / log( ub_plevs(ku(i,k))/ub_plevs(kl(i,k)) )
         761           0 :                       exit
         762             :                    end if
         763             :                 end do
         764             :              end if
         765             :           end do
         766             :        end do
         767             :
         768             :        !--------------------------------------------------------
         769             :        !        ... Find max level less than 50 mb
         770             :        !           fix UB vals from top of model to this level
         771             :        !--------------------------------------------------------
         772           0 :        do i = 1,ncol
         773           0 :           do k = 2,pver
         774           0 :              if( pmid(i,k) > 50.e2_r8 ) then
         775           0 :                 kmax(i) = k
         776           0 :                 exit
         777             :              end if
         778             :           end do
         779             :        end do
         780             :        !--------------------------------------------------------
         781             :        !        ... set the mixing ratio at upper boundary
         782             :        !--------------------------------------------------------
         783           0 :        m = table_h2o_ndx
         784           0 :        do i = 1,ncol
         785           0 :           do k = 1,kmax(i)
         786           0 :              pint_vals(1) = mr_ub(i,m,last,kl(i,k),lchnk) &
         787             :                   + delp(i,k) &
         788           0 :                   * (mr_ub(i,m,last,ku(i,k),lchnk) &
         789           0 :                   - mr_ub(i,m,last,kl(i,k),lchnk))
         790           0 :              pint_vals(2) = mr_ub(i,m,next,kl(i,k),lchnk) &
         791             :                   + delp(i,k) &
         792             :                   * (mr_ub(i,m,next,ku(i,k),lchnk) &
         793           0 :                   - mr_ub(i,m,next,kl(i,k),lchnk))
         794             :              h2o(i,k) = pint_vals(1) &
         795           0 :                   + dels * (pint_vals(2) - pint_vals(1))
         796             :           end do
         797             :        end do
         798             :
         799           0 :        col_loop2 : do i = 1,ncol
         800             :           !--------------------------------------------------------
         801             :           !     ... relax lower stratosphere to extended ubc
         802             :           !           check to make sure ubc is not being imposed too low
         803             :           !           levrelax = lowest model level (highest pressure)
         804             :           !                      in which to relax to ubc
         805             :           !--------------------------------------------------------
         806           0 :           levrelax = ltrop(i)
         807           0 :           do while( pmid(i,levrelax) > ub_plevs(ub_nlevs) )
         808           0 :              levrelax = levrelax - 1
         809             :           end do
         810             : #ifdef DEBUG
         811           0 :           if( levrelax /= ltrop(i) ) then
         812           0 :              write(iulog,*) 'warning -- raised ubc: ',i,          &
         813           0 :                 ltrop(i)-1,nint(pmid(i,ltrop(i)-1)/100._r8),'mb -->', &
         814           0 :                 levrelax,nint(pmid(i,levrelax)/100._r8),'mb'
         815             :           end if
         816             : #endif
         817           0 :           do k = kmax(i)+1,levrelax
         818           0 :              pint_vals(1) = mr_ub(i,m,last,kl(i,k),lchnk) &
         819             :                   + delp(i,k) &
         820           0 :                   * (mr_ub(i,m,last,ku(i,k),lchnk) &
         821           0 :                   - mr_ub(i,m,last,kl(i,k),lchnk))
         822           0 :              pint_vals(2) = mr_ub(i,m,next,kl(i,k),lchnk) &
         823             :                   + delp(i,k) &
         824             :                   * (mr_ub(i,m,next,ku(i,k),lchnk) &
         825           0 :                   - mr_ub(i,m,next,kl(i,k),lchnk))
         826             :              vmrrelax = pint_vals(1) &
         827           0 :                   + dels * (pint_vals(2) - pint_vals(1))
         828           0 :              h2o(i,k) = h2o(i,k) + (vmrrelax - h2o(i,k)) * facrelax
         829             :           end do
         830             :        end do col_loop2
         831             :     end if h2o_overwrite
         832             :
         833           0 :   end subroutine set_fstrat_h2o
         834             :
         835           0 :   subroutine rebin( nsrc, ntrg, src_x, trg_x, src, trg )
         836             :     !---------------------------------------------------------------
         837             :     !   ... rebin src to trg
         838             :     !---------------------------------------------------------------
         839             :
         840             :     !---------------------------------------------------------------
         841             :     !   ... dummy arguments
         842             :     !---------------------------------------------------------------
         843             :     integer, intent(in)   :: nsrc                  ! dimension source array
         844             :     integer, intent(in)   :: ntrg                  ! dimension target array
         845             :     real(r8), intent(in)  :: src_x(nsrc+1)         ! source coordinates
         846             :     real(r8), intent(in)  :: trg_x(ntrg+1)         ! target coordinates
         847             :     real(r8), intent(in)  :: src(nsrc)             ! source array
         848             :     real(r8), intent(out) :: trg(ntrg)             ! target array
         849             :
         850             :     !---------------------------------------------------------------
         851             :     !   ... local variables
         852             :     !---------------------------------------------------------------
         853             :     integer  :: i
         854             :     integer  :: si, si1
         855             :     integer  :: sil, siu
         856             :     real(r8) :: y
         857             :     real(r8) :: sl, su
         858             :     real(r8) :: tl, tu
         859             :
         860             :     !---------------------------------------------------------------
         861             :     !   ... check interval overlap
         862             :     !---------------------------------------------------------------
         863             :     !     if( trg_x(1) < src_x(1) .or. trg_x(ntrg+1) > src_x(nsrc+1) ) then
         864             :     !        write(iulog,*) 'rebin: target grid is outside source grid'
         865             :     !        write(iulog,*) '       target grid from ',trg_x(1),' to ',trg_x(ntrg+1)
         866             :     !        write(iulog,*) '       source grid from ',src_x(1),' to ',src_x(nsrc+1)
         867             :     !        call endrun
         868             :     !     end if
         869             :
         870           0 :     do i = 1,ntrg
         871           0 :        tl = trg_x(i)
         872           0 :        if( tl < src_x(nsrc+1) ) then
         873           0 :           do sil = 1,nsrc+1
         874           0 :              if( tl <= src_x(sil) ) then
         875             :                 exit
         876             :              end if
         877             :           end do
         878           0 :           tu = trg_x(i+1)
         879           0 :           do siu = 1,nsrc+1
         880           0 :              if( tu <= src_x(siu) ) then
         881             :                 exit
         882             :              end if
         883             :           end do
         884           0 :           y   = 0._r8
         885           0 :           sil = max( sil,2 )
         886           0 :           siu = min( siu,nsrc+1 )
         887           0 :           do si = sil,siu
         888           0 :              si1 = si - 1
         889           0 :              sl  = max( tl,src_x(si1) )
         890           0 :              su  = min( tu,src_x(si) )
         891           0 :              y   = y + (su - sl)*src(si1)
         892             :           end do
         893           0 :           trg(i) = y/(trg_x(i+1) - trg_x(i))
         894             :        else
         895           0 :           trg(i) = 0._r8
         896             :        end if
         897             :     end do
         898             :
         899           0 :   end subroutine rebin
         900             :
         901             : end module mo_fstrat

Generated by: LCOV version 1.14