LCOV - code coverage report
Current view: | top level - physics/cam - chem_surfvals.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 127 | 273 | 46.5 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 5 | 8 | 62.5 %

              Line data    Source code
           1             :
           2             : module chem_surfvals
           3             :
           4             : !-----------------------------------------------------------------------------------
           5             : ! Purpose: Provides greenhouse gas (ghg) values at the Earth's surface.
           6             : !          These values may be time dependent.
           7             : !
           8             : ! Author: Brian Eaton (assembled module from existing scattered code pieces)
           9             : !-----------------------------------------------------------------------------------
          10             :
          11             :    use shr_kind_mod,   only: r8=>shr_kind_r8
          12             :    use spmd_utils,     only: masterproc
          13             :    use time_manager,   only: get_curr_date, get_start_date, is_end_curr_day, &
          14             :                              timemgr_datediff, get_curr_calday
          15             :    use cam_abortutils, only: endrun
          16             :    use netcdf
          17             :    use error_messages, only: handle_ncerr
          18             :    use cam_logfile,    only: iulog
          19             :    use m_types,        only: time_ramp
          20             :    use constituents,   only: pcnst
          21             :
          22             : !-----------------------------------------------------------------------
          23             : !- module boilerplate --------------------------------------------------
          24             : !-----------------------------------------------------------------------
          25             :    implicit none
          26             :    private                   ! Make default access private
          27             :    save
          28             :
          29             :    public ::&
          30             :       chem_surfvals_readnl,  &! read namelist input
          31             :       chem_surfvals_init,    &! initialize options that depend on namelist input
          32             :       chem_surfvals_set,     &! set ghg surface values when scenario_ghg is 'RAMPED' or 'CHEM_LBC_FILE'
          33             :       chem_surfvals_get,     &! return surface values for: CO2VMR, CO2MMR, CH4VMR
          34             :                               ! N2OVMR, F11VMR, and F12VMR
          35             :       chem_surfvals_co2_rad   ! return co2 for radiation
          36             :
          37             :    public :: flbc_list
          38             :
          39             : ! Private module data
          40             :
          41             :    ! Default values for namelist variables -- now set by build-namelist
          42             :    real(r8) :: o2mmr = .23143_r8               ! o2 mass mixing ratio
          43             :    real(r8) :: co2vmr_rad = -1.0_r8            ! co2 vmr override for radiation
          44             :    real(r8) :: co2vmr = -1.0_r8                ! co2   volume mixing ratio
          45             :    real(r8) :: n2ovmr = -1.0_r8                ! n2o   volume mixing ratio
          46             :    real(r8) :: ch4vmr = -1.0_r8                ! ch4   volume mixing ratio
          47             :    real(r8) :: f11vmr = -1.0_r8                ! cfc11 volume mixing ratio
          48             :    real(r8) :: f12vmr = -1.0_r8                ! cfc12 volume mixing ratio
          49             :    character(len=16) :: scenario_ghg = 'FIXED' ! 'FIXED','RAMPED', 'RAMP_CO2_ONLY', 'CHEM_LBC_FILE'
          50             :    integer  :: rampYear_ghg = 0                ! ramped gases fixed at this year (if > 0)
          51             :    character(len=256) :: bndtvghg = 'NONE'     ! filename for ramped data
          52             :    integer  :: ramp_co2_start_ymd = 0          ! start date for co2 ramping (yyyymmdd)
          53             :    real(r8) :: ramp_co2_annual_rate = 1.0_r8      ! % amount of co2 ramping per yr; default is 1%
          54             :    real(r8) :: ramp_co2_cap = -9999.0_r8          ! co2 ramp cap if rate>0, floor otherwise
          55             :                                                ! as multiple or fraction of inital value
          56             :                                                ! ex. 4.0 => cap at 4x initial co2 setting
          57             :    integer  :: ghg_yearStart_model = 0         ! model start year
          58             :    integer  :: ghg_yearStart_data  = 0         ! data  start year
          59             :
          60             :    logical  :: ghg_use_calendar                ! true => data year = model year
          61             :    logical :: doRamp_ghg    ! true => turn on ramping for ghg
          62             :    logical :: ramp_just_co2 ! true => ramping to be done just for co2 and not other ghg's
          63             :    integer :: fixYear_ghg   ! year at which Ramped gases are fixed
          64             :    integer :: co2_start     ! date at which co2 begins ramping
          65             :    real(r8) :: co2_daily_factor    ! daily multiplier to achieve annual rate of co2 ramp
          66             :    real(r8) :: co2_limit    ! value of co2vmr where ramping ends
          67             :    real(r8) :: co2_base     ! initial co2 volume mixing ratio, before any ramping
          68             :    integer :: ntim = -1               ! number of yearly data values
          69             :    integer,  allocatable :: yrdata(:) ! yearly data values
          70             :    real(r8), allocatable :: co2(:)    ! co2 mixing ratios in ppmv
          71             :    real(r8), allocatable :: ch4(:)    ! ppbv
          72             :    real(r8), allocatable :: n2o(:)    ! ppbv
          73             :    real(r8), allocatable :: f11(:)    ! pptv
          74             :    real(r8), allocatable :: f12(:)    ! pptv
          75             :    real(r8), allocatable :: adj(:)    ! unitless adjustment factor for f11 & f12
          76             :
          77             :    ! fixed lower boundary
          78             :
          79             :    character(len=256) :: flbc_file = 'NONE'
          80             :    character(len=16)  :: flbc_list(pcnst) = ''
          81             :    type(time_ramp)    :: flbc_timing     != time_ramp( "CYCLICAL",  19970101, 0 )
          82             :
          83             : !=========================================================================================
          84             : contains
          85             : !=========================================================================================
          86             :
          87        1536 : subroutine chem_surfvals_readnl(nlfile)
          88             :
          89             :    ! Read chem_surfvals_nl namelist group.
          90             :
          91             :    use namelist_utils,  only: find_group_name
          92             :    use units,           only: getunit, freeunit
          93             :    use spmd_utils,      only: mpicom, mstrid=>masterprocid, mpi_logical, mpi_integer, &
          94             :                               mpi_real8, mpi_character
          95             :
          96             :    character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
          97             :
          98             :    ! Local variables
          99             :    integer :: unitn, ierr, i
         100             :    character(len=*), parameter :: sub = 'chem_surfvals_readnl'
         101             :
         102             :    character(len=8)   :: flbc_type = 'CYCLICAL'     ! 'CYCLICAL' | 'SERIAL' | 'FIXED'
         103             :    integer            :: flbc_cycle_yr = 0
         104             :    integer            :: flbc_fixed_ymd = 0
         105             :    integer            :: flbc_fixed_tod = 0
         106             :
         107             :    namelist /chem_surfvals_nl/ co2vmr, n2ovmr, ch4vmr, f11vmr, f12vmr, &
         108             :                                co2vmr_rad, scenario_ghg, rampyear_ghg, bndtvghg, &
         109             :                                ramp_co2_start_ymd, ramp_co2_annual_rate, ramp_co2_cap, &
         110             :                                ghg_yearStart_model, ghg_yearStart_data
         111             :
         112             :    namelist /chem_surfvals_nl/ flbc_type, flbc_cycle_yr, flbc_fixed_ymd, flbc_fixed_tod, flbc_list, flbc_file
         113             :
         114             :    !-----------------------------------------------------------------------------
         115             :
         116        1536 :    if (masterproc) then
         117           2 :       unitn = getunit()
         118           2 :       open( unitn, file=trim(nlfile), status='old' )
         119           2 :       call find_group_name(unitn, 'chem_surfvals_nl', status=ierr)
         120           2 :       if (ierr == 0) then
         121           2 :          read(unitn, chem_surfvals_nl, iostat=ierr)
         122           2 :          if (ierr /= 0) then
         123           0 :             call endrun(sub // ':: ERROR reading namelist')
         124             :          end if
         125             :       end if
         126           2 :       close(unitn)
         127           2 :       call freeunit(unitn)
         128             :    end if
         129             :
         130        1536 :    call mpi_bcast(co2vmr, 1, mpi_real8, mstrid, mpicom, ierr)
         131        1536 :    if (ierr /= 0) call endrun(sub//" mpi_bcast: co2vmr")
         132        1536 :    call mpi_bcast(n2ovmr, 1, mpi_real8, mstrid, mpicom, ierr)
         133        1536 :    if (ierr /= 0) call endrun(sub//" mpi_bcast: n2ovmr")
         134        1536 :    call mpi_bcast(ch4vmr, 1, mpi_real8, mstrid, mpicom, ierr)
         135        1536 :    if (ierr /= 0) call endrun(sub//" mpi_bcast: ch4vmr")
         136        1536 :    call mpi_bcast(f11vmr, 1, mpi_real8, mstrid, mpicom, ierr)
         137        1536 :    if (ierr /= 0) call endrun(sub//" mpi_bcast: f11vmr")
         138        1536 :    call mpi_bcast(f12vmr, 1, mpi_real8, mstrid, mpicom, ierr)
         139        1536 :    if (ierr /= 0) call endrun(sub//" mpi_bcast: f12vmr")
         140        1536 :    call mpi_bcast(co2vmr_rad, 1, mpi_real8, mstrid, mpicom, ierr)
         141        1536 :    if (ierr /= 0) call endrun(sub//" mpi_bcast: co2vmr_rad")
         142        1536 :    call mpi_bcast(scenario_ghg, len(scenario_ghg), mpi_character, mstrid, mpicom, ierr)
         143        1536 :    if (ierr /= 0) call endrun(sub//" mpi_bcast: scenario_ghg")
         144        1536 :    call mpi_bcast(rampyear_ghg, 1, mpi_integer, mstrid, mpicom, ierr)
         145        1536 :    if (ierr /= 0) call endrun(sub//" mpi_bcast: rampyear_ghg")
         146        1536 :    call mpi_bcast(bndtvghg, len(bndtvghg), mpi_character, mstrid, mpicom, ierr)
         147        1536 :    if (ierr /= 0) call endrun(sub//" mpi_bcast: bndtvghg")
         148        1536 :    call mpi_bcast(ramp_co2_start_ymd, 1, mpi_integer, mstrid, mpicom, ierr)
         149        1536 :    if (ierr /= 0) call endrun(sub//" mpi_bcast: ramp_co2_start_ymd")
         150        1536 :    call mpi_bcast(ramp_co2_annual_rate, 1, mpi_real8, mstrid, mpicom, ierr)
         151        1536 :    if (ierr /= 0) call endrun(sub//" mpi_bcast: ramp_co2_annual_rate")
         152        1536 :    call mpi_bcast(ramp_co2_cap, 1, mpi_real8, mstrid, mpicom, ierr)
         153        1536 :    if (ierr /= 0) call endrun(sub//" mpi_bcast: ramp_co2_cap")
         154        1536 :    call mpi_bcast(ghg_yearstart_model, 1, mpi_integer, mstrid, mpicom, ierr)
         155        1536 :    if (ierr /= 0) call endrun(sub//" mpi_bcast: ghg_yearstart_model")
         156        1536 :    call mpi_bcast(ghg_yearstart_data, 1, mpi_integer, mstrid, mpicom, ierr)
         157        1536 :    if (ierr /= 0) call endrun(sub//" mpi_bcast: ghg_yearstart_data")
         158        1536 :    call mpi_bcast(flbc_type, len(flbc_type), mpi_character, mstrid, mpicom, ierr)
         159        1536 :    if (ierr /= 0) call endrun(sub//" mpi_bcast: flbc_type")
         160        1536 :    call mpi_bcast(flbc_cycle_yr, 1, mpi_integer, mstrid, mpicom, ierr)
         161        1536 :    if (ierr /= 0) call endrun(sub//" mpi_bcast: flbc_cycle_yr")
         162        1536 :    call mpi_bcast(flbc_fixed_ymd, 1, mpi_integer, mstrid, mpicom, ierr)
         163        1536 :    if (ierr /= 0) call endrun(sub//" mpi_bcast: flbc_fixed_ymd")
         164        1536 :    call mpi_bcast(flbc_fixed_tod, 1, mpi_integer, mstrid, mpicom, ierr)
         165        1536 :    if (ierr /= 0) call endrun(sub//" mpi_bcast: flbc_fixed_tod")
         166        1536 :    call mpi_bcast(flbc_list, len(flbc_list(1))*pcnst, mpi_character, mstrid, mpicom, ierr)
         167        1536 :    if (ierr /= 0) call endrun(sub//" mpi_bcast: flbc_list")
         168        1536 :    call mpi_bcast(flbc_file, len(flbc_file), mpi_character, mstrid, mpicom, ierr)
         169        1536 :    if (ierr /= 0) call endrun(sub//" mpi_bcast: flbc_file")
         170             :
         171        1536 :    flbc_timing%type      = flbc_type
         172        1536 :    flbc_timing%cycle_yr  = flbc_cycle_yr
         173        1536 :    flbc_timing%fixed_ymd = flbc_fixed_ymd
         174        1536 :    flbc_timing%fixed_tod = flbc_fixed_tod
         175             :
         176        1536 :    if ( (bndtvghg.ne.'NONE') .and. (flbc_file.ne.'NONE') ) then
         177           0 :       call endrun(sub//': Cannot specify both bndtvghg and flbc_file ')
         178             :    endif
         179             :
         180        1536 :    if (masterproc) then
         181           2 :       write(iulog,*) ' '
         182           2 :       write(iulog,*) sub//': Settings for control of GHG surface values '
         183           2 :       write(iulog,*) '  scenario_ghg = '//trim(scenario_ghg)
         184             :
         185           2 :       if (scenario_ghg == 'FIXED' .or. scenario_ghg == 'RAMP_CO2_ONLY') then
         186             :
         187           0 :          if (scenario_ghg == 'RAMP_CO2_ONLY') then
         188           0 :             write(iulog,*) '  CO2 will be ramped as follows:'
         189           0 :             write(iulog,*) '    Initial co2vmr       = ', co2vmr
         190           0 :             write(iulog,*) '    ramp_co2_annual_rate = ', ramp_co2_annual_rate
         191           0 :             write(iulog,*) '    ramp_co2_cap         = ', ramp_co2_cap
         192           0 :             if (ramp_co2_start_ymd == 0) then
         193           0 :                write(iulog,*) '    ramp starts at initial run time '
         194             :             else
         195           0 :                write(iulog,*) '    ramp_co2_start_ymd   = ', ramp_co2_start_ymd
         196             :             end if
         197             :
         198             :          else
         199           0 :             write(iulog,*) '  CO2 will be fixed:'
         200           0 :             write(iulog,*) '    co2vmr = ', co2vmr
         201             :          end if
         202             :
         203           0 :          write(iulog,*) '  Other GHG values will be fixed as follows:'
         204           0 :          write(iulog,*) '    n2ovmr = ', n2ovmr
         205           0 :          write(iulog,*) '    ch4vmr = ', ch4vmr
         206           0 :          write(iulog,*) '    f11vmr = ', f11vmr
         207           0 :          write(iulog,*) '    f12vmr = ', f12vmr
         208             :
         209           2 :       else if (scenario_ghg == 'RAMPED') then
         210           0 :          write(iulog,*) '  GHG values time interpolated from global annual averages:'
         211           0 :          write(iulog,*) '    bndtvghg = ', trim(bndtvghg)
         212           0 :          if (rampYear_ghg > 0) then
         213           0 :             write(iulog,*) '  Data from bndtvghg FIXED at year ', rampYear_ghg
         214             :          end if
         215             :
         216           2 :       else if (scenario_ghg == 'CHEM_LBC_FILE') then
         217           2 :          write(iulog,*) '  GHG values time interpolated from LBC file:'
         218           2 :          write(iulog,*) '    flbc_file = ', trim(flbc_file)
         219           2 :          write(iulog,*) '    flbc_type = ', trim(flbc_type)
         220           2 :          if (flbc_type == 'CYCLICAL') then
         221           0 :             write(iulog,*) '    flbc_cycle_yr = ', flbc_cycle_yr
         222           2 :          else if (flbc_type == 'FIXED') then
         223           0 :             write(iulog,*) '    flbc_fixed_ymd = ', flbc_fixed_ymd
         224           0 :             write(iulog,*) '    flbc_fixed_tod = ', flbc_fixed_tod
         225             :          end if
         226           2 :          write(iulog,*) '  Species from LBC file:'
         227          14 :          do i = 1, pcnst
         228          14 :             if (flbc_list(i) == ' ') exit
         229          14 :             write(iulog,*) '    '//trim(flbc_list(i))
         230             :          end do
         231             :
         232             :       else
         233             :          call endrun (sub//': scenario_ghg must be set to either FIXED, RAMPED, RAMP_CO2_ONLY, &
         234           0 :                    & or CHEM_LBC_FILE')
         235             :
         236             :       end if
         237             :
         238           2 :       if (co2vmr_rad > 0._r8) then
         239           0 :          write(iulog,*) '  *** The CO2 prescribed by scenario_ghg has been OVERRIDDEN by ***'
         240           0 :          write(iulog,*) '      co2vmr_rad = ', co2vmr_rad
         241             :       end if
         242             :
         243           2 :       write(iulog,*) ' '
         244             :    end if
         245             :
         246        1536 : end subroutine chem_surfvals_readnl
         247             :
         248             : !================================================================================================
         249             :
         250        1536 : subroutine chem_surfvals_init()
         251             :
         252             : !-----------------------------------------------------------------------
         253             : !
         254             : ! Purpose:
         255             : ! Initialize the ramp options that are controlled by namelist input.
         256             : ! Set surface values at initial time.
         257             : ! N.B. This routine must be called after the time manager has been initialized
         258             : !      since chem_surfvals_set calls time manager methods.
         259             : !
         260             : ! Author: B. Eaton - merged code from parse_namelist and rampnl_ghg.
         261             : !
         262             : !-----------------------------------------------------------------------
         263             :
         264             :    use infnan,       only: posinf, assignment(=)
         265             :    use mo_flbc,      only: flbc_inti
         266             :    use phys_control, only: use_simple_phys
         267             :
         268             :    !---------------------------Local variables-----------------------------
         269             :    integer :: yr, mon, day, ncsec
         270             :    character(len=*), parameter :: sub = 'chem_surfvals_init'
         271             :    !-----------------------------------------------------------------------
         272             :
         273        1536 :    if (use_simple_phys) return
         274             :
         275        1536 :    if (scenario_ghg == 'FIXED') then
         276           0 :       doRamp_ghg = .false.
         277           0 :       ramp_just_co2 = .false.
         278             :
         279        1536 :    else if (scenario_ghg == 'RAMPED') then
         280           0 :       doRamp_ghg = .true.
         281           0 :       ramp_just_co2 = .false.
         282           0 :       call ghg_ramp_read
         283             :
         284           0 :       fixYear_ghg = rampYear_ghg
         285             :
         286           0 :       call chem_surfvals_set()
         287             :
         288        1536 :    else if (scenario_ghg == 'RAMP_CO2_ONLY') then
         289             :
         290           0 :       if(ramp_co2_start_ymd == 0) then
         291             :          ! by default start the ramp at the initial run time
         292           0 :          call get_start_date(yr, mon, day, ncsec)
         293           0 :          ramp_co2_start_ymd = yr*10000 + mon*100 + day
         294             :       end if
         295           0 :       co2_start = ramp_co2_start_ymd
         296             :
         297           0 :       if(ramp_co2_annual_rate <= -100.0_r8) then
         298           0 :          write(iulog,*) 'RAMP_CO2:  invalid ramp_co2_annual_rate= ',ramp_co2_annual_rate
         299           0 :          call endrun (sub//': RAMP_CO2_ANNUAL_RATE must be greater than -100.0')
         300             :       end if
         301             :
         302           0 :       doRamp_ghg = .true.
         303           0 :       ramp_just_co2 = .true.
         304           0 :       co2_base = co2vmr        ! save initial setting
         305             :
         306           0 :       co2_daily_factor = (ramp_co2_annual_rate*0.01_r8+1.0_r8)**(1.0_r8/365.0_r8)
         307             :
         308           0 :       if (ramp_co2_cap > 0.0_r8) then
         309           0 :          co2_limit = ramp_co2_cap * co2_base
         310             :       else                                  ! if no cap/floor specified, provide default
         311           0 :          if (ramp_co2_annual_rate < 0.0_r8) then
         312           0 :             co2_limit = 0.0_r8
         313             :          else
         314           0 :             co2_limit = posinf
         315             :          end if
         316             :       end if
         317             :
         318           0 :       if ((ramp_co2_annual_rate<0.0_r8 .and. co2_limit>co2_base) .or. &
         319             :           (ramp_co2_annual_rate>0.0_r8 .and. co2_limit<co2_base)) then
         320           0 :          write(iulog,*) sub//': ramp_co2_cap is unreachable'
         321           0 :          write(iulog,*) sub//': ramp_co2_annual_rate= ',ramp_co2_annual_rate,' ramp_co2_cap= ',ramp_co2_cap
         322           0 :          call endrun(sub//':  ramp_co2_annual_rate and ramp_co2_cap incompatible')
         323             :       end if
         324             :
         325           0 :       call chem_surfvals_set()
         326             :
         327        1536 :    else if (scenario_ghg == 'CHEM_LBC_FILE') then
         328             :       ! set by lower boundary conditions file
         329        1536 :       call flbc_inti( flbc_file, flbc_list, flbc_timing, co2vmr, ch4vmr, n2ovmr, f11vmr, f12vmr )
         330        1536 :       call chem_surfvals_set()
         331             :
         332             :    endif
         333             :
         334        1536 :    if (masterproc) then
         335           2 :       write(iulog,*) ' '
         336           2 :       write(iulog,*) 'chem_surfvals_init: Initial ghg surface values:'
         337           2 :       write(iulog,*) '  co2 volume mixing ratio = ', chem_surfvals_co2_rad(vmr_in=.true.)
         338           2 :       write(iulog,*) '  ch4 volume mixing ratio = ', ch4vmr
         339           2 :       write(iulog,*) '  n2o volume mixing ratio = ', n2ovmr
         340           2 :       write(iulog,*) '  f11 volume mixing ratio = ', f11vmr
         341           2 :       write(iulog,*) '  f12 volume mixing ratio = ', f12vmr
         342           2 :       write(iulog,*) ' '
         343             :    end if
         344             :
         345        1536 : end subroutine chem_surfvals_init
         346             :
         347             : !=========================================================================================
         348             :
         349           0 : subroutine ghg_ramp_read()
         350             :
         351             : !-----------------------------------------------------------------------
         352             : !
         353             : ! Purpose:
         354             : ! Read ramped greenhouse gas surface data.
         355             : !
         356             : ! Author: T. Henderson
         357             : !
         358             : !-----------------------------------------------------------------------
         359             :
         360        1536 :    use ioFileMod, only: getfil
         361             : #if ( defined SPMD )
         362             :    use mpishorthand, only: mpicom, mpiint, mpir8
         363             : #endif
         364             :    character(len=*), parameter :: subname = 'ghg_ramp_read'
         365             :
         366             : !---------------------------Local variables-----------------------------
         367             :    integer :: ncid
         368             :    integer :: co2_id
         369             :    integer :: ch4_id
         370             :    integer :: n2o_id
         371             :    integer :: f11_id
         372             :    integer :: f12_id
         373             :    integer :: adj_id
         374             :    integer :: date_id
         375             :    integer :: time_id
         376             :    integer :: ierror
         377             :    character(len=256) :: locfn          ! netcdf local filename to open
         378             :
         379           0 :    if (masterproc) then
         380           0 :      call getfil (bndtvghg, locfn, 0)
         381           0 :      call handle_ncerr( nf90_open (trim(locfn), NF90_NOWRITE, ncid),subname,__LINE__)
         382             :
         383           0 :      write(iulog,*)'GHG_RAMP_READ:  reading ramped greenhouse gas surface data from file ',trim(locfn)
         384             :
         385           0 :      call handle_ncerr( nf90_inq_varid( ncid, 'date', date_id ),subname,__LINE__)
         386           0 :      call handle_ncerr( nf90_inq_varid( ncid, 'CO2', co2_id ),subname,__LINE__)
         387           0 :      call handle_ncerr( nf90_inq_varid( ncid, 'CH4', ch4_id ),subname,__LINE__)
         388           0 :      call handle_ncerr( nf90_inq_varid( ncid, 'N2O', n2o_id ),subname,__LINE__)
         389           0 :      call handle_ncerr( nf90_inq_varid( ncid, 'f11', f11_id ),subname,__LINE__)
         390           0 :      call handle_ncerr( nf90_inq_varid( ncid, 'f12', f12_id ),subname,__LINE__)
         391           0 :      call handle_ncerr( nf90_inq_varid( ncid, 'adj', adj_id ),subname,__LINE__)
         392           0 :      call handle_ncerr( nf90_inq_dimid( ncid, 'time', time_id ),subname,__LINE__)
         393           0 :      call handle_ncerr( nf90_inquire_dimension( ncid, time_id, len=ntim ),subname,__LINE__)
         394             :
         395             :    endif
         396             : #if (defined SPMD )
         397           0 :    call mpibcast (ntim, 1, mpiint, 0, mpicom)
         398             : #endif
         399             :    ! these arrays are never deallocated
         400             :    allocate ( yrdata(ntim), co2(ntim), ch4(ntim), n2o(ntim),    &
         401           0 :                  f11(ntim), f12(ntim), adj(ntim), stat=ierror )
         402           0 :    if (ierror /= 0) then
         403           0 :      write(iulog,*)'GHG_RAMP_READ:  ERROR, allocate() failed!'
         404           0 :      call endrun
         405             :    endif
         406           0 :    if (masterproc) then
         407           0 :      call handle_ncerr( nf90_get_var (ncid, date_id, yrdata ),subname,__LINE__)
         408           0 :      yrdata = yrdata / 10000
         409           0 :      call handle_ncerr( nf90_get_var (ncid, co2_id, co2 ),subname,__LINE__)
         410           0 :      call handle_ncerr( nf90_get_var (ncid, ch4_id, ch4 ),subname,__LINE__)
         411           0 :      call handle_ncerr( nf90_get_var (ncid, n2o_id, n2o ),subname,__LINE__)
         412           0 :      call handle_ncerr( nf90_get_var (ncid, f11_id, f11 ),subname,__LINE__)
         413           0 :      call handle_ncerr( nf90_get_var (ncid, f12_id, f12 ),subname,__LINE__)
         414           0 :      call handle_ncerr( nf90_get_var (ncid, adj_id, adj ),subname,__LINE__)
         415           0 :      call handle_ncerr( nf90_close (ncid),subname,__LINE__)
         416           0 :      write(iulog,*)'GHG_RAMP_READ:  successfully read ramped greenhouse gas surface data from years ',&
         417           0 :         yrdata(1),' through ',yrdata(ntim)
         418             :    endif
         419             : #if (defined SPMD )
         420           0 :    call mpibcast (co2, ntim, mpir8, 0, mpicom)
         421           0 :    call mpibcast (ch4, ntim, mpir8, 0, mpicom)
         422           0 :    call mpibcast (n2o, ntim, mpir8, 0, mpicom)
         423           0 :    call mpibcast (f11, ntim, mpir8, 0, mpicom)
         424           0 :    call mpibcast (f12, ntim, mpir8, 0, mpicom)
         425           0 :    call mpibcast (adj, ntim, mpir8, 0, mpicom)
         426           0 :    call mpibcast (yrdata, ntim, mpiint, 0, mpicom)
         427             : #endif
         428             :
         429           0 :    return
         430             :
         431             : end subroutine ghg_ramp_read
         432             :
         433             : !=========================================================================================
         434             :
         435    15420432 : function chem_surfvals_get(name)
         436             :   use physconst,    only: mwdry, mwco2
         437             :
         438             :   character(len=*), intent(in) :: name
         439             :
         440             :   real(r8) :: rmwco2
         441             :   real(r8) :: chem_surfvals_get
         442             :
         443    15420432 :   rmwco2 = mwco2/mwdry    ! ratio of molecular weights of co2 to dry air
         444     1495368 :   select case (name)
         445             :   case ('CO2VMR')
         446     1495368 :      chem_surfvals_get = co2vmr
         447             :   case ('CO2MMR')
         448           0 :      chem_surfvals_get = rmwco2 * co2vmr
         449             :   case ('N2OVMR')
         450     3107424 :      chem_surfvals_get = n2ovmr
         451             :   case ('CH4VMR')
         452     3107424 :      chem_surfvals_get = ch4vmr
         453             :   case ('F11VMR')
         454     3107424 :      chem_surfvals_get = f11vmr
         455             :   case ('F12VMR')
         456     3107424 :      chem_surfvals_get = f12vmr
         457             :   case ('O2MMR')
         458     1495368 :      chem_surfvals_get = o2mmr
         459             :   case default
         460    15420432 :      call endrun('chem_surfvals_get does not know name')
         461             :   end select
         462             :
         463    15420432 : end function chem_surfvals_get
         464             :
         465             :
         466             : !=========================================================================================
         467             :
         468     1618250 : function chem_surfvals_co2_rad(vmr_in)
         469             :
         470             :    ! Return the value of CO2 (as mmr) that is radiatively active.
         471             :
         472             :    ! This method is used by ghg_data to set the prescribed value of CO2 in
         473             :    ! the physics buffer.  If the user has set the co2vmr_rad namelist
         474             :    ! variable then that value will override either the value set by the
         475             :    ! co2vmr namelist variable, or the values time interpolated from a
         476             :    ! dataset.
         477             :
         478             :    ! This method is also used by cam_history to write the radiatively active
         479             :    ! CO2 to the history file.  The optional argument allows returning the
         480             :    ! value as vmr.
         481             :
         482             :    use physconst,    only: mwdry, mwco2
         483             :
         484             :    ! Arguments
         485             :    logical, intent(in), optional :: vmr_in  ! return CO2 as vmr
         486             :
         487             :    ! Return value
         488             :    real(r8) :: chem_surfvals_co2_rad
         489             :
         490             :    ! Local variables
         491             :    real(r8) :: convert_vmr      ! convert vmr to desired output
         492             :    !-----------------------------------------------------------------------
         493             :
         494             :    ! by default convert vmr to mmr
         495     1618250 :    convert_vmr = mwco2/mwdry    ! ratio of molecular weights of co2 to dry air
         496     1618250 :    if (present(vmr_in)) then
         497             :       ! if request return vmr
         498      122882 :       if (vmr_in) convert_vmr = 1.0_r8
         499             :    end if
         500             :
         501     1618250 :    if (co2vmr_rad > 0._r8) then
         502           0 :       chem_surfvals_co2_rad = convert_vmr * co2vmr_rad
         503             :    else
         504     1618250 :       chem_surfvals_co2_rad = convert_vmr * co2vmr
         505             :    end if
         506             :
         507     1618250 : end function chem_surfvals_co2_rad
         508             :
         509             : !=========================================================================================
         510             :
         511      372480 : subroutine chem_surfvals_set()
         512             :
         513             :    use ppgrid,         only: begchunk, endchunk
         514             :    use mo_flbc,        only: flbc_gmean_vmr, flbc_chk
         515             :    use scamMod,        only: single_column, scmiop_flbc_inti, use_camiop
         516             :
         517             : !---------------------------Local variables-----------------------------
         518             :
         519             :    integer  :: yr, mon, day, ncsec ! components of a date
         520             :    integer  :: ncdate              ! current date in integer format [yyyymmdd]
         521             :
         522      372480 :    if ( doRamp_ghg ) then
         523           0 :       if(ramp_just_co2) then
         524           0 :          call chem_surfvals_set_co2()
         525             :       else
         526           0 :          call chem_surfvals_set_all()
         527             :       end if
         528      372480 :    elseif (scenario_ghg == 'CHEM_LBC_FILE') then
         529             :       ! set mixing ratios from cam-chem/waccm lbc file
         530      372480 :       call flbc_chk()
         531      372480 :       if (single_column .and. use_camiop) then
         532           0 :          call scmiop_flbc_inti( co2vmr, ch4vmr, n2ovmr, f11vmr, f12vmr )
         533             :       else
         534             :          ! set by lower boundary conditions file
         535      372480 :          call flbc_gmean_vmr(co2vmr,ch4vmr,n2ovmr,f11vmr,f12vmr)
         536             :       endif
         537             :    endif
         538             :
         539      372480 :    if (masterproc .and. is_end_curr_day()) then
         540          12 :       call get_curr_date(yr, mon, day, ncsec)
         541          12 :       ncdate = yr*10000 + mon*100 + day
         542          12 :       write(iulog,*) 'chem_surfvals_set: ncdate= ',ncdate,' co2vmr=',co2vmr
         543             :
         544          12 :       if (.not. ramp_just_co2 .and. mon==1 .and. day==1) then
         545           2 :          write(iulog,*) 'chem_surfvals_set: ch4vmr=', ch4vmr, ' n2ovmr=', n2ovmr, &
         546           4 :                         ' f11vmr=', f11vmr, ' f12vmr=', f12vmr
         547             :       end if
         548             :
         549             :    end if
         550             :
         551      372480 :    return
         552      372480 : end subroutine chem_surfvals_set
         553             :
         554             : !=========================================================================================
         555             :
         556           0 : subroutine chem_surfvals_set_all()
         557             : !-----------------------------------------------------------------------
         558             : !
         559             : ! Purpose:
         560             : ! Computes greenhouse gas volume mixing ratios via interpolation of
         561             : ! yearly input data.
         562             : !
         563             : ! Author: B. Eaton - updated ramp_ghg for use in chem_surfvals module
         564             : !
         565             : !-----------------------------------------------------------------------
         566      372480 :    use interpolate_data, only: get_timeinterp_factors
         567             :
         568             : !---------------------------Local variables-----------------------------
         569             :
         570             :    integer yrmodel           ! model year
         571             :    integer nyrm              ! year index
         572             :    integer nyrp              ! year index
         573             :    integer :: yr, mon, day   ! components of a date
         574             :    integer :: ncdate         ! current date in integer format [yyyymmdd]
         575             :    integer :: ncsec          ! current time of day [seconds]
         576             :
         577             :    real(r8) :: calday            ! current calendar day
         578             :    real(r8) doymodel             ! model day of year
         579             :    real(r8) doydatam             ! day of year for input data yrdata(nyrm)
         580             :    real(r8) doydatap             ! day or year for input data yrdata(nyrp)
         581             :    real(r8) deltat               ! delta time
         582             :    real(r8) fact1, fact2         ! time interpolation factors
         583             :    real(r8) cfcscl               ! cfc scale factor for f11
         584             :
         585             :    integer yearRan_model         ! model ran year
         586             : !
         587             : ! ---------------------------------------------------------------------
         588             : !
         589           0 :    calday = get_curr_calday()
         590           0 :    call get_curr_date(yr, mon, day, ncsec)
         591           0 :    ncdate = yr*10000 + mon*100 + day
         592             : !
         593             : ! determine ghg_use_calendar
         594             : !
         595           0 :    if ( ghg_yearStart_model > 0 .and. ghg_yearStart_data > 0 ) then
         596           0 :       ghg_use_calendar = .false.
         597             :    else
         598           0 :       ghg_use_calendar = .true.
         599             :    end if
         600             : !
         601             : ! determine index into input data
         602             : !
         603           0 :    if ( fixYear_ghg > 0) then
         604           0 :       yrmodel  = fixYear_ghg
         605           0 :       nyrm = fixYear_ghg - yrdata(1) + 1
         606             :    else
         607           0 :       if ( ghg_use_calendar) then
         608           0 :          yrmodel  = yr
         609           0 :          nyrm = yr - yrdata(1) + 1
         610             :       else
         611           0 :          yearRan_model = yr - ghg_yearStart_model
         612           0 :          if ( yearRan_model < 0 ) then
         613           0 :             call endrun('chem_surfvals_set_all: incorrect ghg_yearStart_model')
         614             :          endif
         615           0 :          yrmodel  = yearRan_model + ghg_yearStart_data
         616             :
         617           0 :          nyrm = ghg_yearStart_data + yearRan_model - yrdata(1) + 1
         618             :       end if
         619             :    end if
         620             :
         621           0 :    nyrp       = nyrm + 1
         622             : !
         623             : ! if current date is before yrdata(1), quit
         624             : !
         625           0 :    if (nyrm < 1) then
         626           0 :       write(iulog,*)'chem_surfvals_set_all: data time index is out of bounds'
         627           0 :       write(iulog,*)'nyrm = ',nyrm,' nyrp= ',nyrp, ' ncdate= ', ncdate
         628           0 :       call endrun
         629             :    endif
         630             : !
         631             : ! if current date later than yrdata(ntim), call endrun.
         632             : ! if want to use ntim values - uncomment the following lines
         633             : ! below and comment the call to endrun and previous write
         634             : !
         635           0 :    if (nyrp > ntim) then
         636           0 :       call endrun ('chem_surfvals_set_all: error - current date is past the end of valid data')
         637             : !         write(iulog,*)'chem_surfvals_set_all: using ghg data for ',yrdata(ntim)
         638             : !         co2vmr = co2(ntim)*1.e-06
         639             : !         ch4vmr = ch4(ntim)*1.e-09
         640             : !         n2ovmr = n2o(ntim)*1.e-09
         641             : !         f11vmr = f11(ntim)*1.e-12*(1.+cfcscl)
         642             : !         f12vmr = f12(ntim)*1.e-12
         643             : !         co2mmr = rmwco2 * co2vmr
         644             : !         return
         645             :    endif
         646             : !
         647             : ! determine time interpolation factors, check sanity
         648             : ! of interpolation factors to within 32-bit roundoff
         649             : ! assume that day of year is 1 for all input data
         650             : !
         651           0 :    doymodel = yrmodel*365._r8    + calday
         652           0 :    doydatam = yrdata(nyrm)*365._r8 + 1._r8
         653           0 :    doydatap = yrdata(nyrp)*365._r8 + 1._r8
         654             :
         655             :    call get_timeinterp_factors(.false.,2,doydatam,doydatap, doymodel, &
         656           0 :         fact1, fact2,'chem_surfvals')
         657             :
         658             : !
         659             : ! do time interpolation:
         660             : !   co2     in ppmv
         661             : !   n2o,ch4 in ppbv
         662             : !   f11,f12 in pptv
         663             : !
         664           0 :    co2vmr = (co2(nyrm)*fact1 + co2(nyrp)*fact2)*1.e-06_r8
         665           0 :    ch4vmr = (ch4(nyrm)*fact1 + ch4(nyrp)*fact2)*1.e-09_r8
         666           0 :    n2ovmr = (n2o(nyrm)*fact1 + n2o(nyrp)*fact2)*1.e-09_r8
         667             :
         668           0 :    cfcscl = (adj(nyrm)*fact1 + adj(nyrp)*fact2)
         669           0 :    f11vmr = (f11(nyrm)*fact1 + f11(nyrp)*fact2)*1.e-12_r8*(1._r8+cfcscl)
         670           0 :    f12vmr = (f12(nyrm)*fact1 + f12(nyrp)*fact2)*1.e-12_r8
         671             :
         672           0 :    return
         673             : end subroutine chem_surfvals_set_all
         674             :
         675             : !=========================================================================================
         676             :
         677           0 : subroutine chem_surfvals_set_co2()
         678             : !-----------------------------------------------------------------------
         679             : !
         680             : ! Purpose:
         681             : ! Computes co2 greenhouse gas volume mixing ratio via ramping info
         682             : ! provided in namelist var's
         683             : !
         684             : ! Author: B. Eaton - updated ramp_ghg for use in chem_surfvals module
         685             : !
         686             : !-----------------------------------------------------------------------
         687             :    use shr_kind_mod, only: r8 => shr_kind_r8
         688             :
         689             : !---------------------------Local variables-----------------------------
         690             :
         691             :    real(r8) :: daydiff             ! number of days of co2 ramping
         692             :    integer  :: yr, mon, day, ncsec ! components of a date
         693             :    integer  :: ncdate              ! current date in integer format [yyyymmdd]
         694             : !-----------------------------------------------------------------------
         695             :
         696           0 :    call get_curr_date(yr, mon, day, ncsec)
         697           0 :    ncdate = yr*10000 + mon*100 + day
         698             :
         699           0 :    call timemgr_datediff(co2_start, 0, ncdate, ncsec, daydiff)
         700             :
         701           0 :    if (daydiff > 0.0_r8) then
         702             :
         703           0 :       co2vmr = co2_base*(co2_daily_factor)**daydiff
         704             :
         705           0 :       if(co2_daily_factor < 1.0_r8) then
         706           0 :          co2vmr = max(co2vmr,co2_limit)
         707             :       else
         708           0 :          co2vmr = min(co2vmr,co2_limit)
         709             :       end if
         710             :    end if
         711             :
         712           0 :    return
         713             : end subroutine chem_surfvals_set_co2
         714             :
         715             :
         716             : !=========================================================================================
         717             :
         718             : end module chem_surfvals

Generated by: LCOV version 1.14