LCOV - code coverage report
Current view: | top level - control - scamMod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 7 | 621 | 1.1 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 1 | 7 | 14.3 %

              Line data    Source code
           1             : module scamMod
           2             :   !----------------------------------------------------------------------
           3             :   !
           4             :   ! This module provides control variables and namelist functionality
           5             :   ! for SCAM.
           6             :   !
           7             :   ! As with global CAM, SCAM is initialized with state information
           8             :   ! from the initial and boundary files. For each succeeding timestep
           9             :   ! SCAM calculates the physics tendencies and combines these with
          10             :   ! the advective tendencies provided by the IOP file to produce
          11             :   ! a forcast. Generally, the control variables in this module
          12             :   ! determine what data and parameterizations will be used to make
          13             :   ! the forecast. For instance, many of the control variables in
          14             :   ! this module provide flexibility to affect the forecast by overriding
          15             :   ! parameterization prognosed tendencies with observed tendencies
          16             :   ! of a particular field program recorded on the IOP file.
          17             :   !
          18             :   ! Public functions/subroutines:
          19             :   !   scam_readnl
          20             :   !-----------------------------------------------------------------------
          21             :
          22             : use shr_kind_mod,   only: r8 => shr_kind_r8, cl => shr_kind_cl
          23             : use spmd_utils,     only: masterproc,npes
          24             : use pmgrid,         only: plon, plat, plev, plevp
          25             : use constituents,   only: cnst_get_ind, pcnst, cnst_name
          26             : use netcdf,         only: NF90_NOERR,NF90_CLOSE,NF90_GET_VAR,NF90_INQUIRE_DIMENSION, &
          27             :                           NF90_INQ_DIMID, NF90_INQ_VARID, NF90_NOWRITE, NF90_OPEN, &
          28             :                           NF90_GET_ATT,NF90_GLOBAL,NF90_INQUIRE_ATTRIBUTE, &
          29             :                           NF90_INQUIRE_VARIABLE, NF90_MAX_VAR_DIMS, nf90_get_var
          30             : use shr_scam_mod,   only: shr_scam_getCloseLatLon
          31             : use cam_logfile,    only: iulog
          32             : use cam_abortutils, only: endrun
          33             : use time_manager,   only: get_curr_date, get_nstep,is_first_step,get_start_date,timemgr_time_inc
          34             : use error_messages, only: handle_ncerr
          35             :
          36             :
          37             : implicit none
          38             : private
          39             :
          40             : ! PUBLIC INTERFACES:
          41             :
          42             : public :: scam_readnl         ! read SCAM namelist options
          43             : public :: readiopdata         ! read iop boundary data
          44             : public :: setiopupdate        ! find index in iopboundary data for current time
          45             : public :: plevs0              ! Define the pressures of the interfaces and midpoints
          46             : public :: scmiop_flbc_inti
          47             : public :: setiopupdate_init
          48             :
          49             : ! PUBLIC MODULE DATA:
          50             :
          51             : real(r8), public ::  pressure_levels(plev)
          52             : real(r8), public ::  scmlat   ! input namelist latitude for scam
          53             : real(r8), public ::  scmlon   ! input namelist longitude for scam
          54             : real(r8), public ::  closeioplat   ! closest iop latitude for scam
          55             : real(r8), public ::  closeioplon   ! closest iop longitude for scam
          56             : integer,  public ::  closeioplatidx   ! file array index of closest iop latitude for scam
          57             : integer,  public ::  closeioplonidx   ! file array index closest iop longitude for scam
          58             :
          59             :
          60             : integer, parameter :: num_switches = 20
          61             : integer, parameter :: max_path_len = 128
          62             :
          63             : logical, public ::  single_column         ! Using IOP file or not
          64             : logical, public ::  use_iop               ! Using IOP file or not
          65             : logical, public ::  use_pert_init         ! perturb initial values
          66             : logical, public ::  use_pert_frc          ! perturb forcing
          67             : logical, public ::  switch(num_switches)  ! Logical flag settings from GUI
          68             : logical, public ::  l_uvphys              ! If true, update u/v after TPHYS
          69             : logical, public ::  l_uvadvect            ! If true, T, U & V will be passed to SLT
          70             : logical, public ::  l_conv                ! use flux divergence terms for T and q?
          71             : logical, public ::  l_divtr               ! use flux divergence terms for constituents?
          72             : logical, public ::  l_diag                ! do we want available diagnostics?
          73             :
          74             : integer, public ::  error_code            ! Error code from netCDF reads
          75             : integer, public ::  initTimeIdx
          76             : integer, public ::  seedval
          77             : integer :: bdate, last_date, last_sec
          78             :
          79             : character(len=max_path_len), public ::  modelfile
          80             : character(len=max_path_len), public ::  analysisfile
          81             : character(len=max_path_len), public ::  sicfile
          82             : character(len=max_path_len), public ::  userfile
          83             : character(len=max_path_len), public ::  sstfile
          84             : character(len=max_path_len), public ::  lsmpftfile
          85             : character(len=max_path_len), public ::  pressfile
          86             : character(len=max_path_len), public ::  topofile
          87             : character(len=max_path_len), public ::  ozonefile
          88             : character(len=max_path_len), public ::  iopfile
          89             : character(len=max_path_len), public ::  absemsfile
          90             : character(len=max_path_len), public ::  aermassfile
          91             : character(len=max_path_len), public ::  aeropticsfile
          92             : character(len=max_path_len), public ::  timeinvfile
          93             : character(len=max_path_len), public ::  lsmsurffile
          94             : character(len=max_path_len), public ::  lsminifile
          95             :
          96             : ! note that scm_zadv_q is set to slt to be consistent with CAM BFB testing
          97             :
          98             :
          99             : character(len=16), public    :: scm_zadv_T  = 'eulc            '
         100             : character(len=16), public    :: scm_zadv_q  = 'slt             '
         101             : character(len=16), public    :: scm_zadv_uv = 'eulc            '
         102             :
         103             : real(r8), public ::  fixmascam
         104             : real(r8), public ::  betacam
         105             : real(r8), public ::  alphacam(pcnst)
         106             : real(r8), public ::  dqfxcam(plon,plev,pcnst)
         107             :
         108             : real(r8), public ::      divq3d(plev,pcnst)  ! 3D q advection
         109             : real(r8), public ::      divt3d(plev)        ! 3D T advection
         110             : real(r8), public ::      divu3d(plev)        ! 3D U advection
         111             : real(r8), public ::      divv3d(plev)        ! 3D V advection
         112             : real(r8), public ::      vertdivq(plev,pcnst)! vertical q advection
         113             : real(r8), public ::      vertdivt(plev)      ! vertical T advection
         114             : real(r8), public ::      vertdivu(plev)      ! vertical T advection
         115             : real(r8), public ::      vertdivv(plev)      ! vertical T advection
         116             : real(r8), public ::      ptend               ! surface pressure tendency
         117             : real(r8), public ::      qdiff(plev)         ! model minus observed humidity
         118             : real(r8), public ::      qobs(plev)          ! actual W.V. Mixing ratio
         119             : real(r8), public ::      qinitobs(plev,pcnst)! initial tracer field
         120             : real(r8), public ::      cldliqobs(plev)     ! actual W.V. Mixing ratio
         121             : real(r8), public ::      cldiceobs(plev)     ! actual W.V. Mixing ratio
         122             : real(r8), public ::      numliqobs(plev)     ! actual
         123             : real(r8), public ::      numiceobs(plev)     ! actual
         124             : real(r8), public ::      precobs(1)          ! observed precipitation
         125             : real(r8), public ::      lhflxobs(1)         ! observed surface latent heat flux
         126             : real(r8), public ::      heat_glob_scm(1)    ! observed heat total
         127             : real(r8), public ::      shflxobs(1)         ! observed surface sensible heat flux
         128             : real(r8), public ::      q1obs(plev)         ! observed apparent heat source
         129             : real(r8), public ::      q2obs(plev)         ! observed apparent heat sink
         130             : real(r8), public ::      tdiff(plev)         ! model minus observed temp
         131             : real(r8), public ::      tground(1)          ! ground temperature
         132             : real(r8), public ::      psobs               ! observed surface pressure
         133             : real(r8), public ::      tobs(plev)          ! observed temperature
         134             : real(r8), public ::      tsair(1)            ! air temperature at the surface
         135             : real(r8), public ::      udiff(plev)         ! model minus observed uwind
         136             : real(r8), public ::      uobs(plev)          ! actual u wind
         137             : real(r8), public ::      vdiff(plev)         ! model minus observed vwind
         138             : real(r8), public ::      vobs(plev)          ! actual v wind
         139             : real(r8), public ::      cldobs(plev)        ! observed cld
         140             : real(r8), public ::      clwpobs(plev)       ! observed clwp
         141             : real(r8), public ::      aldirobs(1)         ! observed aldir
         142             : real(r8), public ::      aldifobs(1)         ! observed aldif
         143             : real(r8), public ::      asdirobs(1)         ! observed asdir
         144             : real(r8), public ::      asdifobs(1)         ! observed asdif
         145             :
         146             : real(r8), public ::      co2vmrobs(1)        ! observed co2vmr
         147             : real(r8), public ::      ch4vmrobs(1)        ! observed ch3vmr
         148             : real(r8), public ::      n2ovmrobs(1)        ! observed n2ovmr
         149             : real(r8), public ::      f11vmrobs(1)        ! observed f11vmr
         150             : real(r8), public ::      f12vmrobs(1)        ! observed f12vmr
         151             : real(r8), public ::      soltsiobs(1)        ! observed solar
         152             :
         153             : real(r8), public ::      wfld(plev)          ! Vertical motion (slt)
         154             : real(r8), public ::      wfldh(plevp)        ! Vertical motion (slt)
         155             : real(r8), public ::      divq(plev,pcnst)    ! Divergence of moisture
         156             : real(r8), public ::      divt(plev)          ! Divergence of temperature
         157             : real(r8), public ::      divu(plev)          ! Horiz Divergence of E/W
         158             : real(r8), public ::      divv(plev)          ! Horiz Divergence of N/S
         159             :                                              ! mo_drydep algorithm
         160             : real(r8), public, pointer :: loniop(:)
         161             : real(r8), public, pointer :: latiop(:)
         162             :
         163             : integer, public ::     iopTimeIdx            ! index into iop dataset
         164             : integer, public ::     steplength            ! Length of time-step
         165             : integer, public ::     base_date             ! Date in (yyyymmdd) of start time
         166             : integer, public ::     base_secs             ! Time of day of start time (sec)
         167             :
         168             : ! SCAM public data defaults
         169             :
         170             : logical, public ::  doiopupdate            = .false. ! do we need to read next iop timepoint
         171             : logical, public ::  have_lhflx             = .false. ! dataset contains lhflx
         172             : logical, public ::  have_shflx             = .false. ! dataset contains shflx
         173             : logical, public ::  have_heat_glob         = .false. ! dataset contains heat total
         174             : logical, public ::  have_tg                = .false. ! dataset contains tg
         175             : logical, public ::  have_tsair             = .false. ! dataset contains tsair
         176             : logical, public ::  have_divq              = .false. ! dataset contains divq
         177             : logical, public ::  have_divt              = .false. ! dataset contains divt
         178             : logical, public ::  have_divq3d            = .false. ! dataset contains divq3d
         179             : logical, public ::  have_vertdivu          = .false. ! dataset contains vertdivu
         180             : logical, public ::  have_vertdivv          = .false. ! dataset contains vertdivv
         181             : logical, public ::  have_vertdivt          = .false. ! dataset contains vertdivt
         182             : logical, public ::  have_vertdivq          = .false. ! dataset contains vertdivq
         183             : logical, public ::  have_divt3d            = .false. ! dataset contains divt3d
         184             : logical, public ::  have_divu3d            = .false. ! dataset contains divu3d
         185             : logical, public ::  have_divv3d            = .false. ! dataset contains divv3d
         186             : logical, public ::  have_divu              = .false. ! dataset contains divu
         187             : logical, public ::  have_divv              = .false. ! dataset contains divv
         188             : logical, public ::  have_omega             = .false. ! dataset contains omega
         189             : logical, public ::  have_phis              = .false. ! dataset contains phis
         190             : logical, public ::  have_ptend             = .false. ! dataset contains ptend
         191             : logical, public ::  have_ps                = .false. ! dataset contains ps
         192             : logical, public ::  have_q                 = .false. ! dataset contains q
         193             : logical, public ::  have_q1                = .false. ! dataset contains Q1
         194             : logical, public ::  have_q2                = .false. ! dataset contains Q2
         195             : logical, public ::  have_prec              = .false. ! dataset contains prec
         196             : logical, public ::  have_t                 = .false. ! dataset contains t
         197             : logical, public ::  have_u                 = .false. ! dataset contains u
         198             : logical, public ::  have_v                 = .false. ! dataset contains v
         199             : logical, public ::  have_cld               = .false. ! dataset contains cld
         200             : logical, public ::  have_cldliq            = .false. ! dataset contains cldliq
         201             : logical, public ::  have_cldice            = .false. ! dataset contains cldice
         202             : logical, public ::  have_numliq            = .false. ! dataset contains numliq
         203             : logical, public ::  have_numice            = .false. ! dataset contains numice
         204             : logical, public ::  have_clwp              = .false. ! dataset contains clwp
         205             : logical, public ::  have_aldir             = .false. ! dataset contains aldir
         206             : logical, public ::  have_aldif             = .false. ! dataset contains aldif
         207             : logical, public ::  have_asdir             = .false. ! dataset contains asdir
         208             : logical, public ::  have_asdif             = .false. ! dataset contains asdif
         209             : logical, public ::  use_camiop             = .false. ! use cam generated forcing
         210             : logical, public ::  use_3dfrc              = .false. ! use 3d forcing
         211             : logical, public ::  isrestart              = .false. ! If this is a restart step or not
         212             :
         213             : ! SCAM namelist defaults
         214             :
         215             : logical, public ::  scm_backfill_iop_w_init = .false. ! Backfill missing IOP data from initial file
         216             : logical, public ::  scm_relaxation         = .false. ! Use relaxation
         217             : logical, public ::  scm_crm_mode           = .false. ! Use column radiation mode
         218             : logical, public ::  scm_cambfb_mode        = .false. ! Use extra CAM IOP fields to assure bit for bit match with CAM run
         219             : logical, public ::  scm_use_obs_T          = .false. ! Use the SCAM-IOP observed T at each timestep instead of forecasting.
         220             : logical, public ::  scm_force_latlon       = .false. ! force scam to use the lat lon fields specified in the namelist not closest
         221             : real(r8), public              ::  scm_relaxation_low      ! lowest level to apply relaxation
         222             : real(r8), public              ::  scm_relaxation_high     ! highest level to apply relaxation
         223             : real(r8), public              ::  scm_relax_top_p         = 0._r8       ! upper bound for scm relaxation
         224             : real(r8), public              ::  scm_relax_bot_p         = huge(1._r8) ! lower bound for scm relaxation
         225             : real(r8), public              ::  scm_relax_tau_sec       = 10800._r8   ! relaxation time constant (sec)
         226             :
         227             : ! +++BPM:
         228             : ! modification... allow a linear ramp in relaxation time scale:
         229             : logical, public :: scm_relax_linear = .false.
         230             : real(r8), public    :: scm_relax_tau_bot_sec = 10800._r8
         231             : real(r8), public    :: scm_relax_tau_top_sec = 10800._r8
         232             : character(len=26), public  :: scm_relax_fincl(pcnst)
         233             :
         234             : !
         235             : ! note that scm_use_obs_uv is set to true to be consistent with CAM BFB testing
         236             : !
         237             :
         238             : logical, public ::  scm_use_obs_uv         = .true. ! Use the SCAM-IOP observed u,v at each time step instead of forecasting.
         239             :
         240             : logical, public ::  scm_use_obs_qv         = .false. ! Use the SCAM-IOP observed qv  at each time step instead of forecasting.
         241             : logical, public ::  scm_use_3dfrc          = .false. ! Use CAMIOP 3d forcing if true, else use dycore vertical plus horizontal
         242             : logical, public ::  scm_iop_lhflxshflxTg   = .false. !turn off LW rad
         243             : logical, public ::  scm_iop_Tg             = .false. !turn off LW rad
         244             :
         245             : character(len=200), public ::  scm_clubb_iop_name   ! IOP name for CLUBB
         246             :
         247             : integer, allocatable, public :: tsec(:)
         248             : integer, public :: ntime
         249             :
         250             : !=======================================================================
         251             : contains
         252             : !=======================================================================
         253             :
         254        1536 : subroutine scam_readnl(nlfile,single_column_in,scmlat_in,scmlon_in)
         255             :
         256             :   use namelist_utils,  only: find_group_name
         257             :   use units,           only: getunit, freeunit
         258             :   use dycore,          only: dycore_is
         259             :   use wrap_nf,         only: wrap_open
         260             :
         261             :
         262             : !---------------------------Arguments-----------------------------------
         263             :
         264             :   character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input  (nlfile=atm_in)
         265             :   logical,          intent(in) :: single_column_in
         266             :   real(r8),         intent(in) :: scmlat_in
         267             :   real(r8),         intent(in) :: scmlon_in
         268             :
         269             :   ! Local variables
         270             :   character(len=*), parameter :: sub = 'scam_readnl'
         271             :   integer :: unitn, ierr, i
         272             :   integer  :: ncid
         273             :   integer  :: iatt
         274             :   logical  :: adv
         275             :
         276             : ! this list should include any variable that you might want to include in the namelist
         277             :   namelist /scam_nl/ iopfile, scm_iop_lhflxshflxTg, scm_iop_Tg, scm_relaxation, &
         278             :        scm_relax_top_p,scm_relax_bot_p,scm_relax_tau_sec, &
         279             :        scm_cambfb_mode,scm_crm_mode,scm_zadv_uv,scm_zadv_T,scm_zadv_q,&
         280             :        scm_use_obs_T, scm_use_obs_uv, scm_use_obs_qv, scm_use_3dfrc, &
         281             :        scm_relax_linear, scm_relax_tau_top_sec, &
         282             :        scm_relax_tau_bot_sec, scm_force_latlon, scm_relax_fincl, &
         283             :        scm_backfill_iop_w_init
         284             :
         285        1536 :   single_column=single_column_in
         286             :
         287        1536 :   iopfile            = ' '
         288        1536 :   scm_clubb_iop_name = ' '
         289       64512 :   scm_relax_fincl(:) = ' '
         290        1536 :   if( single_column ) then
         291           0 :      if( npes>1) call endrun('SCAM_READNL: SCAM doesnt support using more than 1 pe.')
         292             :
         293           0 :      if ( .not. (dycore_is('EUL') .or. dycore_is('SE')) .or. plon /= 1 .or. plat /=1 ) then
         294           0 :         call endrun('SCAM_SETOPTS: must compile model for SCAM mode when namelist parameter single_column is .true.')
         295             :      endif
         296             :
         297           0 :      scmlat=scmlat_in
         298           0 :      scmlon=scmlon_in
         299             :
         300           0 :      if( scmlat < -90._r8 .or. scmlat > 90._r8 ) then
         301           0 :         call endrun('SCAM_READNL: SCMLAT must be between -90. and 90. degrees.')
         302           0 :      elseif( scmlon < 0._r8 .or. scmlon > 360._r8 ) then
         303           0 :         call endrun('SCAM_READNL: SCMLON must be between 0. and 360. degrees.')
         304             :      end if
         305             :
         306             :      ! Read namelist
         307           0 :      if (masterproc) then
         308           0 :         unitn = getunit()
         309           0 :         open( unitn, file=trim(nlfile), status='old' )
         310           0 :         call find_group_name(unitn, 'scam_nl', status=ierr)
         311           0 :         if (ierr == 0) then
         312           0 :            read(unitn, scam_nl, iostat=ierr)
         313           0 :            if (ierr /= 0) then
         314           0 :               call endrun(sub // ':: ERROR reading namelist')
         315             :            end if
         316             :         end if
         317           0 :         close(unitn)
         318           0 :         call freeunit(unitn)
         319             :      end if
         320             :
         321             :      ! Error checking:
         322             :
         323           0 :      iopfile = trim(iopfile)
         324           0 :      if( iopfile /= "" ) then
         325           0 :         use_iop = .true.
         326             :      else
         327           0 :         call endrun('SCAM_READNL: must specify IOP file for single column mode')
         328             :      endif
         329             :
         330           0 :      call wrap_open( iopfile, NF90_NOWRITE, ncid )
         331             :
         332           0 :      if( nf90_inquire_attribute( ncid, NF90_GLOBAL, 'CAM_GENERATED_FORCING', iatt ) ==  NF90_NOERR ) then
         333           0 :         use_camiop = .true.
         334             :      else
         335           0 :         use_camiop = .false.
         336             :      endif
         337             :
         338             :      ! If we are not forcing the lat and lon from the namelist use the closest lat and lon that is found in the IOP file.
         339           0 :      if (.not.scm_force_latlon) then
         340           0 :         call shr_scam_GetCloseLatLon( ncid, scmlat, scmlon, closeioplat, closeioplon, closeioplatidx, closeioplonidx )
         341           0 :         write(iulog,*) 'SCAM_READNL: using closest IOP column to lat/lon specified in drv_in'
         342           0 :         write(iulog,*) '   requested lat,lon    =',scmlat,', ',scmlon
         343           0 :         write(iulog,*) '   closest IOP lat,lon  =',closeioplat,', ',closeioplon
         344           0 :         scmlat = closeioplat
         345           0 :         scmlon = closeioplon
         346             :      end if
         347             :
         348           0 :      if (masterproc) then
         349           0 :         write (iulog,*) 'Single Column Model Options: '
         350           0 :         write (iulog,*) '============================='
         351           0 :         write (iulog,*) '  iopfile                     = ',trim(iopfile)
         352           0 :         write (iulog,*) '  scm_backfill_iop_w_init     = ',scm_backfill_iop_w_init
         353           0 :         write (iulog,*) '  scm_cambfb_mode             = ',scm_cambfb_mode
         354           0 :         write (iulog,*) '  scm_crm_mode                = ',scm_crm_mode
         355           0 :         write (iulog,*) '  scm_force_latlon            = ',scm_force_latlon
         356           0 :         write (iulog,*) '  scm_iop_Tg                  = ',scm_iop_Tg
         357           0 :         write (iulog,*) '  scm_iop_lhflxshflxTg        = ',scm_iop_lhflxshflxTg
         358           0 :         write (iulog,*) '  scm_relaxation              = ',scm_relaxation
         359           0 :         write (iulog,*) '  scm_relax_bot_p             = ',scm_relax_bot_p
         360           0 :         write (iulog,*) '  scm_relax_linear            = ',scm_relax_linear
         361           0 :         write (iulog,*) '  scm_relax_tau_bot_sec       = ',scm_relax_tau_bot_sec
         362           0 :         write (iulog,*) '  scm_relax_tau_sec           = ',scm_relax_tau_sec
         363           0 :         write (iulog,*) '  scm_relax_tau_top_sec       = ',scm_relax_tau_top_sec
         364           0 :         write (iulog,*) '  scm_relax_top_p             = ',scm_relax_top_p
         365           0 :         write (iulog,*) '  scm_use_obs_T               = ',scm_use_obs_T
         366           0 :         write (iulog,*) '  scm_use_3dfrc               = ',scm_use_3dfrc
         367           0 :         write (iulog,*) '  scm_use_obs_qv              = ',scm_use_obs_qv
         368           0 :         write (iulog,*) '  scm_use_obs_uv              = ',scm_use_obs_uv
         369           0 :         write (iulog,*) '  scm_zadv_T                  = ',trim(scm_zadv_T)
         370           0 :         write (iulog,*) '  scm_zadv_q                  = ',trim(scm_zadv_q)
         371           0 :         write (iulog,*) '  scm_zadv_uv                 = ',trim(scm_zadv_uv)
         372           0 :         write (iulog,*) '  scm_relax_finc: '
         373             :         ! output scm_relax_fincl character array
         374           0 :         do i=1,pcnst
         375           0 :            if (scm_relax_fincl(i) /= '') then
         376           0 :               adv = mod(i,4)==0
         377           0 :               if (adv) then
         378           0 :                  write (iulog, "(A18)") "'"//trim(scm_relax_fincl(i))//"',"
         379             :               else
         380           0 :                  write (iulog, "(A18)", ADVANCE="NO") "'"//trim(scm_relax_fincl(i))//"',"
         381             :               end if
         382             :            else
         383             :               exit
         384             :            end if
         385             :         end do
         386           0 :         print *
         387             :      end if
         388             :   end if
         389             :
         390        1536 : end subroutine scam_readnl
         391           0 : subroutine readiopdata(hyam, hybm, hyai, hybi, ps0)
         392             : !-----------------------------------------------------------------------
         393             : !
         394             : !     Open and read netCDF file containing initial IOP  conditions
         395             : !
         396             : !---------------------------Code history--------------------------------
         397             : !
         398             : !     Written by J.  Truesdale    August, 1996, revised January, 1998
         399             : !
         400             : !-----------------------------------------------------------------------
         401             :         use getinterpnetcdfdata, only: getinterpncdata
         402             :         use string_utils,        only: to_lower
         403             :         use wrap_nf,             only: wrap_inq_dimid,wrap_get_vara_realx
         404             : !-----------------------------------------------------------------------
         405             :    implicit none
         406             :
         407             :    character(len=*), parameter ::  sub = "read_iop_data"
         408             : !
         409             : !------------------------------Input Arguments--------------------------
         410             : !
         411             :    real(r8),intent(in) :: hyam(plev),hybm(plev),hyai(plevp),hybi(plevp),ps0
         412             : !
         413             : !------------------------------Locals-----------------------------------
         414             : !
         415             :    integer :: NCID, status
         416             :    integer :: time_dimID, lev_dimID,  lev_varID, varid
         417             :    integer :: i,j
         418             :    integer :: nlev
         419             :    integer :: total_levs
         420             :    integer :: u_attlen
         421             :
         422             :    integer :: k, m
         423             :    integer :: icldliq,icldice
         424             :    integer :: inumliq,inumice
         425             :
         426             :    logical :: have_srf              ! value at surface is available
         427             :    logical :: fill_ends             !
         428             :    logical :: have_cnst(pcnst)
         429             :    real(r8) :: dummy
         430             :    real(r8) :: srf(1)                  ! value at surface
         431             :    real(r8) :: hyamiop(plev)  ! a hybrid coef midpoint
         432             :    real(r8) :: hybmiop(plev)  ! b hybrid coef midpoint
         433             :    real(r8) :: pmid(plev)  ! pressure at model levels (time n)
         434             :    real(r8) :: pint(plevp) ! pressure at model interfaces (n  )
         435             :    real(r8) :: pdel(plev)  ! pdel(k)   = pint  (k+1)-pint  (k)
         436             :    real(r8) :: weight
         437             :    real(r8) :: tmpdata(1)
         438             :    real(r8) :: coldata(plev)
         439           0 :    real(r8), allocatable :: dplevs( : )
         440             :    integer :: strt4(4),cnt4(4)
         441             :    integer :: nstep
         442             :    integer :: ios
         443             :    character(len=128) :: units ! Units
         444             :
         445           0 :    nstep = get_nstep()
         446           0 :    fill_ends= .false.
         447             :
         448             : !
         449             : !     Open IOP dataset
         450             : !
         451             :    call handle_ncerr( nf90_open (iopfile, 0, ncid),&
         452           0 :        'ERROR - scamMod.F90:readiopdata', __LINE__)
         453             :
         454             : !
         455             : !     if the dataset is a CAM generated dataset set use_camiop to true
         456             : !       CAM IOP datasets have a global attribute called CAM_GENERATED_IOP
         457             : !
         458           0 :    if ( nf90_inquire_attribute( ncid, NF90_GLOBAL, 'CAM_GENERATED_FORCING', attnum=i )== NF90_NOERR ) then
         459           0 :       use_camiop = .true.
         460             :    else
         461           0 :       use_camiop = .false.
         462             :    endif
         463             :
         464             : !=====================================================================
         465             : !
         466             : !     Read time variables
         467             :
         468             :
         469           0 :    status = nf90_inq_dimid (ncid, 'time', time_dimID )
         470           0 :    if (status /= NF90_NOERR) then
         471           0 :       status = nf90_inq_dimid (ncid, 'tsec', time_dimID )
         472           0 :       if (status /= NF90_NOERR) then
         473           0 :          if (masterproc) write(iulog,*) sub//':ERROR - Could not find dimension ID for time/tsec'
         474           0 :          status = NF90_CLOSE ( ncid )
         475           0 :          call endrun(sub // ':ERROR - time/tsec must be present on the IOP file.')
         476             :       end if
         477             :    end if
         478             :
         479             :    call handle_ncerr( nf90_inquire_dimension( ncid, time_dimID, len=ntime ),&
         480           0 :          'Error - scamMod.F90:readiopdata unable to find time dimension', __LINE__)
         481             :
         482             : !
         483             : !======================================================
         484             : !     read level data
         485             : !
         486           0 :    status = NF90_INQ_DIMID( ncid, 'lev', lev_dimID )
         487           0 :    if ( status /= nf90_noerr ) then
         488           0 :       if (masterproc) write(iulog,*) sub//':ERROR - Could not find variable dim ID  for lev'
         489           0 :       status = NF90_CLOSE ( ncid )
         490           0 :       call endrun(sub // ':ERROR - Could not find variable dim ID for lev')
         491             :    end if
         492             :
         493             :    call handle_ncerr( nf90_inquire_dimension( ncid, lev_dimID, len=nlev ),&
         494           0 :          'Error - scamMod.f90:readiopdata unable to find level dimension', __LINE__)
         495             :
         496           0 :    allocate(dplevs(nlev+1),stat=ios)
         497           0 :    if( ios /= 0 ) then
         498           0 :       write(iulog,*) sub//':ERROR: failed to allocate dplevs; error = ',ios
         499           0 :       call endrun(sub//':ERROR:readiopdata failed to allocate dplevs')
         500             :    end if
         501             :
         502           0 :    status = NF90_INQ_VARID( ncid, 'lev', lev_varID )
         503           0 :    if ( status /= nf90_noerr ) then
         504           0 :       if (masterproc) write(iulog,*) sub//':ERROR - scamMod.F90:readiopdata:Could not find variable ID for lev'
         505           0 :       status = NF90_CLOSE ( ncid )
         506           0 :       call endrun(sub//':ERROR:ould not find variable ID for lev')
         507             :    end if
         508             :
         509           0 :    call handle_ncerr( nf90_get_var (ncid, lev_varID, dplevs(:nlev)),&
         510           0 :                     'Error - scamMod.F90:readiopdata unable to read pressure levels', __LINE__)
         511             : !
         512             : !CAM generated forcing already has pressure on millibars convert standard IOP if needed.
         513             : !
         514             :    call handle_ncerr(nf90_inquire_attribute(ncid, lev_varID, 'units', len=u_attlen),&
         515           0 :                     'Error - scamMod.F90:readiopdata unable to find units attribute', __LINE__)
         516             :    call handle_ncerr(nf90_get_att(ncid, lev_varID, 'units', units),&
         517           0 :                     'Error - scamMod.F90:readiopdata unable to read units attribute', __LINE__)
         518           0 :    units=trim(to_lower(units(1:u_attlen)))
         519             :
         520           0 :    if ( units=='pa' .or. units=='pascal' .or. units=='pascals' ) then
         521             : !
         522             : !     convert pressure from Pascals to Millibars ( lev is expressed in pascals in iop datasets )
         523             : !
         524           0 :       do i=1,nlev
         525           0 :          dplevs( i ) = dplevs( i )/100._r8
         526             :       end do
         527             :    endif
         528             :
         529           0 :    status = nf90_inq_varid( ncid, 'Ps', varid   )
         530           0 :    if ( status /= nf90_noerr ) then
         531           0 :       have_ps= .false.
         532           0 :       if (masterproc) write(iulog,*) sub//':Could not find variable Ps'
         533           0 :       if ( .not. scm_backfill_iop_w_init ) then
         534           0 :          status = NF90_CLOSE( ncid )
         535           0 :          call endrun(sub//':ERROR :IOP file must contain Surface Pressure (Ps) variable')
         536             :       else
         537           0 :          if ( is_first_step() .and. masterproc) write(iulog,*) 'Using surface pressure value from IC file if present'
         538             :       endif
         539             :    else
         540           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
         541           0 :       status = nf90_get_var(ncid, varid, psobs, strt4)
         542           0 :       have_ps = .true.
         543             :    endif
         544             :
         545             :
         546             : !  If the IOP dataset has hyam,hybm,etc it is assumed to be a hybrid level
         547             : !  dataset
         548             :
         549           0 :    status =  nf90_inq_varid( ncid, 'hyam', varid   )
         550           0 :    if ( status == nf90_noerr .and. have_ps) then
         551           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
         552           0 :       status = nf90_get_var(ncid, varid, hyamiop, strt4)
         553           0 :       status =  nf90_inq_varid( ncid, 'hybm', varid   )
         554           0 :       status = nf90_get_var(ncid, varid, hybmiop, strt4)
         555           0 :       do i = 1, nlev
         556           0 :          dplevs( i ) = 1000.0_r8 * hyamiop( i ) + psobs * hybmiop( i ) / 100.0_r8
         557             :       end do
         558             :    endif
         559             :
         560             : !     add the surface pressure to the pressure level data, so that
         561             : !     surface boundary condition will be set properly,
         562             : !     making sure that it is the highest pressure in the array.
         563             : !
         564             :
         565           0 :    total_levs = nlev+1
         566           0 :    dplevs(nlev+1) = psobs/100.0_r8 ! ps is expressed in pascals
         567           0 :    do i= nlev, 1, -1
         568           0 :       if ( dplevs(i) > psobs/100.0_r8) then
         569           0 :          total_levs = i
         570           0 :          dplevs(i) = psobs/100.0_r8
         571             :       end if
         572             :    end do
         573           0 :    if (.not. use_camiop ) then
         574           0 :       nlev = total_levs
         575             :    endif
         576           0 :    if ( nlev == 1 ) then
         577           0 :       if (masterproc) write(iulog,*) sub//':Error - scamMod.F90:readiopdata: Ps too low!'
         578           0 :       call endrun(sub//':ERROR:Ps value on datasets is incongurent with levs data - mismatch in units?')
         579             :    endif
         580             :
         581             : !=====================================================================
         582             : !get global vmrs from camiop file
         583           0 :    status =  nf90_inq_varid( ncid, 'co2vmr', varid   )
         584           0 :    if ( status == nf90_noerr) then
         585           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
         586           0 :       call wrap_get_vara_realx (ncid,varid,strt4,cnt4,co2vmrobs)
         587             :    else
         588           0 :       if (is_first_step()) write(iulog,*)'using column value of co2vmr from boundary data as global volume mixing ratio'
         589             :    end if
         590           0 :    status =  nf90_inq_varid( ncid, 'ch4vmr', varid   )
         591           0 :    if ( status == nf90_noerr) then
         592           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
         593           0 :       call wrap_get_vara_realx (ncid,varid,strt4,cnt4,ch4vmrobs)
         594             :    else
         595           0 :       if (is_first_step()) write(iulog,*)'using column value of ch4vmr from boundary data as global volume mixing ratio'
         596             :    end if
         597           0 :    status =  nf90_inq_varid( ncid, 'n2ovmr', varid   )
         598           0 :    if ( status == nf90_noerr) then
         599           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
         600           0 :       call wrap_get_vara_realx (ncid,varid,strt4,cnt4,n2ovmrobs)
         601             :    else
         602           0 :       if (is_first_step()) write(iulog,*)'using column value of n2ovmr from boundary data as global volume mixing ratio'
         603             :    end if
         604           0 :    status =  nf90_inq_varid( ncid, 'f11vmr', varid   )
         605           0 :    if ( status == nf90_noerr) then
         606           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
         607           0 :       call wrap_get_vara_realx (ncid,varid,strt4,cnt4,f11vmrobs)
         608             :    else
         609           0 :       if (is_first_step()) write(iulog,*)'using column value of f11vmr from boundary data as global volume mixing ratio'
         610             :    end if
         611           0 :    status =  nf90_inq_varid( ncid, 'f12vmr', varid   )
         612           0 :    if ( status == nf90_noerr) then
         613           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
         614           0 :       call wrap_get_vara_realx (ncid,varid,strt4,cnt4,f12vmrobs)
         615             :    else
         616           0 :       if (is_first_step()) write(iulog,*)'using column value of f12vmr from boundary data as global volume mixing ratio'
         617             :    end if
         618           0 :    status =  nf90_inq_varid( ncid, 'soltsi', varid   )
         619           0 :    if ( status == nf90_noerr) then
         620           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
         621           0 :       call wrap_get_vara_realx (ncid,varid,strt4,cnt4,soltsiobs)
         622             :    else
         623           0 :       if (is_first_step()) write(iulog,*)'using column value of soltsi from boundary data as global solar tsi'
         624             :    end if
         625             : !=====================================================================
         626             : !get state variables from camiop file
         627             :
         628           0 :    status =  nf90_inq_varid( ncid, 'Tsair', varid   )
         629           0 :    if ( status /= nf90_noerr ) then
         630           0 :       have_tsair = .false.
         631             :    else
         632           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
         633           0 :       call wrap_get_vara_realx (ncid,varid,strt4,cnt4,tsair)
         634           0 :       have_tsair = .true.
         635             :    endif
         636             : !
         637             : !      read in Tobs  For cam generated iop readin small t to avoid confusion
         638             : !      with capital T defined in cam
         639             : !
         640           0 :    tobs(:)= 0._r8
         641             :
         642           0 :    if ( use_camiop ) then
         643             :      call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx,'t', have_tsair, &
         644             :           tsair(1), fill_ends, scm_crm_mode, &
         645           0 :           dplevs, nlev,psobs, hyam, hybm,tobs, status )
         646             :    else
         647             :      call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx,'T', have_tsair, &
         648             :           tsair(1), fill_ends, scm_crm_mode, &
         649           0 :           dplevs, nlev,psobs, hyam, hybm, tobs, status )
         650             :    endif
         651           0 :    if ( status /= nf90_noerr ) then
         652           0 :       have_t = .false.
         653           0 :       if (masterproc) write(iulog,*) sub//':Could not find variable T on IOP file'
         654           0 :       if ( scm_backfill_iop_w_init ) then
         655           0 :          if (masterproc) write(iulog,*) sub//':Using value of T(tobs) from IC file if it exists'
         656             :       else
         657           0 :          if (masterproc) write(iulog,*) sub//':set tobs to 0.'
         658             :       endif
         659             : !
         660             : !     set T3 to Tobs on first time step
         661             : !
         662             :    else
         663           0 :       have_t = .true.
         664             :    endif
         665             :
         666           0 :    status = nf90_inq_varid( ncid, 'Tg', varid   )
         667           0 :    if (status /= nf90_noerr) then
         668           0 :       if (masterproc) write(iulog,*) sub//':Could not find variable Tg on IOP dataset'
         669           0 :       if ( have_tsair ) then
         670           0 :          if (masterproc) write(iulog,*) sub//':Using Tsair'
         671           0 :          tground = tsair     ! use surface value from T field
         672           0 :          have_Tg = .true.
         673             :       else
         674           0 :          have_Tg = .true.
         675           0 :          if (masterproc) write(iulog,*) sub//':Using T at lowest level from IOP dataset'
         676           0 :          tground = tobs(plev)
         677             :       endif
         678             :    else
         679           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
         680           0 :       call wrap_get_vara_realx (ncid,varid,strt4,cnt4,tground)
         681           0 :       have_Tg = .true.
         682             :    endif
         683             :
         684           0 :    status = nf90_inq_varid( ncid, 'qsrf', varid   )
         685             :
         686           0 :    if ( status /= nf90_noerr ) then
         687           0 :       have_srf = .false.
         688             :    else
         689           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
         690           0 :       status = nf90_get_var(ncid, varid, srf(1), strt4)
         691           0 :       have_srf = .true.
         692             :    endif
         693             :
         694           0 :    qobs(:)= 0._r8
         695             :    call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx,  'q', have_srf, &
         696             :       srf(1), fill_ends, scm_crm_mode, &
         697           0 :       dplevs, nlev,psobs, hyam, hybm, qobs, status )
         698           0 :    if ( status /= nf90_noerr ) then
         699           0 :       have_q = .false.
         700           0 :       if (masterproc) write(iulog,*) sub//':Could not find variable q on IOP file'
         701           0 :       if ( scm_backfill_iop_w_init ) then
         702           0 :          if (masterproc) write(iulog,*) sub//':Using values for q from IC file if available'
         703             :       else
         704           0 :          if (masterproc) write(iulog,*) sub//':Setting qobs to 0.'
         705             :       endif
         706             :    else
         707           0 :       have_q = .true.
         708             :    endif
         709             :
         710           0 :    cldobs = 0._r8
         711             :    call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx,  'cld', .false., &
         712           0 :       dummy, fill_ends, scm_crm_mode, dplevs, nlev,psobs, hyam, hybm, cldobs, status )
         713           0 :    if ( status /= nf90_noerr ) then
         714           0 :       have_cld = .false.
         715             :    else
         716           0 :       have_cld = .true.
         717             :    endif
         718             :
         719           0 :    clwpobs = 0._r8
         720             :    call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx,  'clwp', .false., &
         721           0 :       dummy, fill_ends, scm_crm_mode, dplevs, nlev,psobs, hyam, hybm, clwpobs, status )
         722           0 :    if ( status /= nf90_noerr ) then
         723           0 :       have_clwp = .false.
         724             :    else
         725           0 :       have_clwp = .true.
         726             :    endif
         727             :
         728             : !
         729             : !       read divq (horizontal advection)
         730             : !
         731           0 :    status = nf90_inq_varid( ncid, 'divqsrf', varid   )
         732           0 :    if ( status /= nf90_noerr ) then
         733           0 :       have_srf = .false.
         734             :    else
         735           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
         736           0 :       status = nf90_get_var(ncid, varid, srf(1), strt4)
         737           0 :       have_srf = .true.
         738             :    endif
         739             :
         740           0 :    divq(:,:)=0._r8
         741             :
         742             :    call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, &
         743             :         'divq', have_srf, srf(1), fill_ends, scm_crm_mode, &
         744           0 :         dplevs, nlev,psobs, hyam, hybm, divq(:,1), status )
         745           0 :    if ( status /= nf90_noerr ) then
         746           0 :       have_divq = .false.
         747             :    else
         748           0 :       have_divq = .true.
         749             :    endif
         750             :
         751             : !
         752             : !     read vertdivq if available
         753             : !
         754           0 :    status = nf90_inq_varid( ncid, 'vertdivqsrf', varid   )
         755           0 :    if ( status /= nf90_noerr ) then
         756           0 :       have_srf = .false.
         757             :    else
         758           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
         759           0 :       status = nf90_get_var(ncid, varid, srf(1), strt4)
         760           0 :       have_srf = .true.
         761             :    endif
         762             :
         763           0 :    vertdivq=0._r8
         764             :
         765             :    call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'vertdivq', &
         766             :         have_srf, srf(1), fill_ends, scm_crm_mode, &
         767           0 :         dplevs, nlev,psobs, hyam, hybm, vertdivq(:,1), status )
         768           0 :    if ( status /= nf90_noerr ) then
         769           0 :       have_vertdivq = .false.
         770             :    else
         771           0 :       have_vertdivq = .true.
         772             :    endif
         773             :
         774           0 :    status = nf90_inq_varid( ncid, 'vertdivqsrf', varid   )
         775           0 :    if ( status /= nf90_noerr ) then
         776           0 :       have_srf = .false.
         777             :    else
         778           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
         779           0 :       status = nf90_get_var(ncid, varid, srf(1), strt4)
         780           0 :       have_srf = .true.
         781             :    endif
         782             : !
         783             : !   add calls to get dynamics tendencies for all prognostic consts
         784             : !
         785           0 :    divq3d=0._r8
         786             :
         787           0 :    do m = 1, pcnst
         788           0 :       call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, trim(cnst_name(m))//'_dten', &
         789             :       have_srf, srf(1), fill_ends, scm_crm_mode, &
         790           0 :       dplevs, nlev,psobs, hyam, hybm, divq3d(:,m), status )
         791           0 :       write(iulog,*)'checking ',trim(cnst_name(m))//'_dten',status
         792           0 :       if ( status /= nf90_noerr ) then
         793           0 :          have_cnst(m) = .false.
         794           0 :          divq3d(1:,m)=0._r8
         795             :       else
         796           0 :          if (m==1) have_divq3d = .true.
         797           0 :          have_cnst(m) = .true.
         798             :       endif
         799             :
         800           0 :       coldata = 0._r8
         801             :       call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, trim(cnst_name(m))//'_dqfx', &
         802             :       have_srf, srf(1), fill_ends, scm_crm_mode, &
         803           0 :       dplevs, nlev,psobs, hyam, hybm, coldata, status )
         804           0 :       if ( STATUS /= NF90_NOERR ) then
         805           0 :          dqfxcam(1,:,m)=0._r8
         806             :       else
         807           0 :          dqfxcam(1,:,m)=coldata(:)
         808             :       endif
         809             :
         810           0 :       tmpdata = 0._r8
         811             :       call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, trim(cnst_name(m))//'_alph', &
         812             :       have_srf, srf(1), fill_ends, scm_crm_mode, &
         813           0 :       dplevs, nlev,psobs, hyam, hybm, tmpdata, status )
         814           0 :       if ( status /= nf90_noerr ) then
         815           0 :          alphacam(m)=0._r8
         816             :       else
         817           0 :           alphacam(m)=tmpdata(1)
         818             :       endif
         819             :
         820             :    end do
         821             :
         822             :
         823           0 :    numliqobs = 0._r8
         824           0 :    call cnst_get_ind('NUMLIQ', inumliq, abort=.false.)
         825           0 :    if ( inumliq > 0 ) then
         826           0 :       have_srf = .false.
         827             :       call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'NUMLIQ', &
         828             :            have_srf, srf(1), fill_ends, scm_crm_mode, &
         829           0 :            dplevs, nlev,psobs, hyam, hybm, numliqobs, status )
         830           0 :       if ( status /= nf90_noerr ) then
         831           0 :          have_numliq = .false.
         832             :       else
         833           0 :          have_numliq = .true.
         834             :       endif
         835             :    else
         836           0 :          have_numliq = .false.
         837             :    end if
         838             :
         839           0 :    have_srf = .false.
         840             :
         841           0 :    cldliqobs = 0._r8
         842           0 :    call cnst_get_ind('CLDLIQ', icldliq, abort=.false.)
         843           0 :    if ( icldliq > 0 ) then
         844             :       call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'CLDLIQ', &
         845             :            have_srf, srf(1), fill_ends, scm_crm_mode, &
         846           0 :            dplevs, nlev,psobs, hyam, hybm, cldliqobs, status )
         847           0 :       if ( status /= nf90_noerr ) then
         848           0 :          have_cldliq = .false.
         849             :       else
         850           0 :          have_cldliq = .true.
         851             :       endif
         852             :    else
         853           0 :          have_cldliq = .false.
         854             :    endif
         855             :
         856           0 :    cldiceobs = 0._r8
         857           0 :    call cnst_get_ind('CLDICE', icldice, abort=.false.)
         858           0 :    if ( icldice > 0 ) then
         859             :       call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'CLDICE', &
         860             :            have_srf, srf(1), fill_ends, scm_crm_mode, &
         861           0 :            dplevs, nlev,psobs, hyam, hybm, cldiceobs, status )
         862           0 :       if ( status /= nf90_noerr ) then
         863           0 :          have_cldice = .false.
         864             :       else
         865           0 :          have_cldice = .true.
         866             :       endif
         867             :    else
         868           0 :       have_cldice = .false.
         869             :    endif
         870             :
         871           0 :    numiceobs = 0._r8
         872           0 :    call cnst_get_ind('NUMICE', inumice, abort=.false.)
         873           0 :    if ( inumice > 0 ) then
         874           0 :       have_srf = .false.
         875             :       call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'NUMICE', &
         876             :          have_srf, srf(1), fill_ends, scm_crm_mode, &
         877           0 :          dplevs, nlev,psobs, hyam, hybm, numiceobs, status )
         878           0 :       if ( status /= nf90_noerr ) then
         879           0 :          have_numice = .false.
         880             :       else
         881           0 :          have_numice = .true.
         882             :       endif
         883             :    else
         884           0 :       have_numice = .false.
         885             :    end if
         886             :
         887             : !
         888             : !       read divu (optional field)
         889             : !
         890           0 :    status = nf90_inq_varid( ncid, 'divusrf', varid   )
         891           0 :    if ( status /= nf90_noerr ) then
         892           0 :       have_srf = .false.
         893             :    else
         894           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
         895           0 :       status = nf90_get_var(ncid, varid, srf(1), strt4)
         896           0 :       have_srf = .true.
         897             :    endif
         898             :
         899           0 :    divu = 0._r8
         900             :    call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'divu', &
         901             :       have_srf, srf(1), fill_ends, scm_crm_mode, &
         902           0 :       dplevs, nlev,psobs, hyam, hybm, divu, status )
         903           0 :    if ( status /= nf90_noerr ) then
         904           0 :       have_divu = .false.
         905             :    else
         906           0 :       have_divu = .true.
         907             :    endif
         908             : !
         909             : !       read divv (optional field)
         910             : !
         911           0 :    status = nf90_inq_varid( ncid, 'divvsrf', varid   )
         912           0 :    if ( status /= nf90_noerr ) then
         913           0 :       have_srf = .false.
         914             :    else
         915           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
         916           0 :       status = nf90_get_var(ncid, varid, srf(1), strt4)
         917           0 :       have_srf = .true.
         918             :    endif
         919             :
         920           0 :    divv = 0._r8
         921             :    call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'divv', &
         922             :       have_srf, srf(1), fill_ends, scm_crm_mode, &
         923           0 :       dplevs, nlev,psobs, hyam, hybm, divv, status )
         924           0 :    if ( status /= nf90_noerr ) then
         925           0 :       have_divv = .false.
         926             :    else
         927           0 :       have_divv = .true.
         928             :    endif
         929             : !
         930             : !       read divt (optional field)
         931             : !
         932           0 :    status = nf90_inq_varid( ncid, 'divtsrf', varid   )
         933           0 :    if ( status /= nf90_noerr ) then
         934           0 :       have_srf = .false.
         935             :    else
         936           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
         937           0 :       status = nf90_get_var(ncid, varid, srf(1), strt4)
         938           0 :       have_srf = .true.
         939             :    endif
         940             :
         941           0 :    divt=0._r8
         942             :    call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, &
         943             :       'divT', have_srf, srf(1), fill_ends, scm_crm_mode, &
         944           0 :       dplevs, nlev,psobs, hyam, hybm, divt, status )
         945           0 :    if ( status /= nf90_noerr ) then
         946           0 :       have_divt = .false.
         947             :    else
         948           0 :       have_divt = .true.
         949             :    endif
         950             :
         951             : !
         952             : !     read vertdivt if available
         953             : !
         954           0 :    status = nf90_inq_varid( ncid, 'vertdivTsrf', varid   )
         955           0 :    if ( status /= nf90_noerr ) then
         956           0 :       have_srf = .false.
         957             :    else
         958           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
         959           0 :       status = nf90_get_var(ncid, varid, srf(1), strt4)
         960           0 :       have_srf = .true.
         961             :    endif
         962             :
         963           0 :    vertdivt=0._r8
         964             :    call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'vertdivTx', &
         965             :       have_srf, srf(1), fill_ends, scm_crm_mode, &
         966           0 :       dplevs, nlev,psobs, hyam, hybm, vertdivt, status )
         967           0 :    if ( status /= nf90_noerr ) then
         968             :       call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'vertdivT', &
         969             :            have_srf, srf(1), fill_ends, scm_crm_mode, &
         970           0 :            dplevs, nlev,psobs, hyam, hybm, vertdivt, status )
         971           0 :       if ( status /= nf90_noerr ) then
         972           0 :          have_vertdivt = .false.
         973             :       else
         974           0 :          have_vertdivt = .true.
         975             :       endif
         976             :    else
         977           0 :       have_vertdivt = .true.
         978             :    endif
         979             : !
         980             : !       read divt3d (combined vertical/horizontal advection)
         981             : !      (optional field)
         982             :
         983           0 :    status = nf90_inq_varid( ncid, 'divT3dsrf', varid   )
         984           0 :    if ( status /= nf90_noerr ) then
         985           0 :       have_srf = .false.
         986             :    else
         987           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
         988           0 :       status = nf90_get_var(ncid, varid, srf(1), strt4)
         989           0 :       have_srf = .true.
         990             :    endif
         991             :
         992           0 :    divT3d = 0._r8
         993             :
         994             :    call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'divT3d', &
         995             :       have_srf, srf(1), fill_ends, scm_crm_mode, &
         996           0 :       dplevs, nlev,psobs, hyam, hybm, divt3d, status )
         997           0 :       write(iulog,*)'checking divT3d:',status,nf90_noerr
         998           0 :    if ( status /= nf90_noerr ) then
         999           0 :       have_divt3d = .false.
        1000             :    else
        1001           0 :       have_divt3d = .true.
        1002             :    endif
        1003             :
        1004           0 :    divU3d = 0._r8
        1005             :
        1006             :    call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'divU3d', &
        1007             :       have_srf, srf(1), fill_ends, scm_crm_mode, &
        1008           0 :       dplevs, nlev,psobs, hyam, hybm, divu3d, status )
        1009           0 :    if ( status /= nf90_noerr ) then
        1010           0 :       have_divu3d = .false.
        1011             :    else
        1012           0 :       have_divu3d = .true.
        1013             :    endif
        1014             :
        1015           0 :    divV3d = 0._r8
        1016             :
        1017             :    call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'divV3d', &
        1018             :       have_srf, srf(1), fill_ends, scm_crm_mode, &
        1019           0 :       dplevs, nlev,psobs, hyam, hybm, divv3d, status )
        1020           0 :    if ( status /= nf90_noerr ) then
        1021           0 :       have_divv3d = .false.
        1022             :    else
        1023           0 :       have_divv3d = .true.
        1024             :    endif
        1025             :
        1026           0 :    status = nf90_inq_varid( ncid, 'Ptend', varid   )
        1027           0 :    if ( status /= nf90_noerr ) then
        1028           0 :       have_ptend = .false.
        1029           0 :       if (masterproc) write(iulog,*) sub//':Could not find variable Ptend. Setting to zero'
        1030           0 :       ptend = 0.0_r8
        1031             :    else
        1032           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
        1033           0 :       status = nf90_get_var(ncid, varid, srf(1), strt4)
        1034           0 :       have_ptend = .true.
        1035           0 :       ptend= srf(1)
        1036             :    endif
        1037             :
        1038           0 :    wfld=0._r8
        1039             :
        1040             :    call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, &
        1041             :       'omega', .true., ptend, fill_ends, scm_crm_mode, &
        1042           0 :       dplevs, nlev,psobs, hyam, hybm, wfld, status )
        1043           0 :    if ( status /= nf90_noerr ) then
        1044           0 :       have_omega = .false.
        1045           0 :       if (masterproc) write(iulog,*) sub//':Could not find variable omega on IOP'
        1046           0 :       if ( scm_backfill_iop_w_init ) then
        1047           0 :          if (masterproc) write(iulog,*) sub//'Using omega from IC file'
        1048             :       else
        1049           0 :          if (masterproc) write(iulog,*) sub//'setting Omega to 0. throughout the column'
        1050             :       endif
        1051             :    else
        1052           0 :       have_omega = .true.
        1053             :    endif
        1054           0 :    call plevs0(plev, psobs, ps0, hyam, hybm, hyai, hybi, pint, pmid ,pdel)
        1055             : !
        1056             : ! Build interface vector for the specified omega profile
        1057             : ! (weighted average in pressure of specified level values)
        1058             : !
        1059           0 :    wfldh(:) = 0.0_r8
        1060             :
        1061           0 :    do k=2,plev
        1062           0 :       weight = (pint(k) - pmid(k-1))/(pmid(k) - pmid(k-1))
        1063           0 :       wfldh(k) = (1.0_r8 - weight)*wfld(k-1) + weight*wfld(k)
        1064             :    end do
        1065             :
        1066           0 :    status = nf90_inq_varid( ncid, 'usrf', varid   )
        1067           0 :    if ( status /= nf90_noerr ) then
        1068           0 :       have_srf = .false.
        1069             :    else
        1070           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
        1071           0 :       call wrap_get_vara_realx (ncid,varid,strt4,cnt4,srf)
        1072           0 :       have_srf = .true.
        1073             :    endif
        1074             :
        1075           0 :    uobs=0._r8
        1076             :
        1077             :    call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, &
        1078             :       'u', have_srf, srf(1), fill_ends, scm_crm_mode, &
        1079           0 :       dplevs, nlev,psobs, hyam, hybm, uobs, status )
        1080           0 :    if ( status /= nf90_noerr ) then
        1081           0 :       have_u = .false.
        1082             :    else
        1083           0 :       have_u = .true.
        1084             :    endif
        1085             :
        1086           0 :    status = nf90_inq_varid( ncid, 'vsrf', varid   )
        1087           0 :    if ( status /= nf90_noerr ) then
        1088           0 :       have_srf = .false.
        1089             :    else
        1090           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
        1091           0 :       call wrap_get_vara_realx (ncid,varid,strt4,cnt4,srf)
        1092           0 :       have_srf = .true.
        1093             :    endif
        1094             :
        1095           0 :    vobs=0._r8
        1096             :
        1097             :    call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, &
        1098             :       'v', have_srf, srf(1), fill_ends, scm_crm_mode, &
        1099           0 :       dplevs, nlev,psobs, hyam, hybm, vobs, status )
        1100           0 :    if ( status /= nf90_noerr ) then
        1101           0 :       have_v = .false.
        1102             :    else
        1103           0 :       have_v = .true.
        1104             :    endif
        1105             :
        1106           0 :    status = nf90_inq_varid( ncid, 'Prec', varid   )
        1107           0 :    if ( status /= nf90_noerr ) then
        1108           0 :       have_prec = .false.
        1109             :    else
        1110           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
        1111           0 :       call wrap_get_vara_realx (ncid,varid,strt4,cnt4,precobs)
        1112           0 :       have_prec = .true.
        1113             :    endif
        1114             :
        1115           0 :    q1obs = 0._r8
        1116             :
        1117             :    call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'Q1', &
        1118             :       .false., dummy, fill_ends, scm_crm_mode, & ! datasets don't contain Q1 at surface
        1119           0 :       dplevs, nlev,psobs, hyam, hybm, q1obs, status )
        1120           0 :    if ( status /= nf90_noerr ) then
        1121           0 :       have_q1 = .false.
        1122             :    else
        1123           0 :       have_q1 = .true.
        1124             :    endif
        1125             :
        1126           0 :    q1obs = 0._r8
        1127             :
        1128             :    call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'Q2', &
        1129             :       .false., dummy, fill_ends, scm_crm_mode, & ! datasets don't contain Q2 at surface
        1130           0 :       dplevs, nlev,psobs, hyam, hybm, q1obs, status )
        1131           0 :    if ( status /= nf90_noerr ) then
        1132           0 :       have_q2 = .false.
        1133             :    else
        1134           0 :       have_q2 = .true.
        1135             :    endif
        1136             :
        1137             : !  Test for BOTH 'lhflx' and 'lh' without overwriting 'have_lhflx'.
        1138             : !  Analagous changes made for the surface heat flux
        1139             :
        1140           0 :    status = nf90_inq_varid( ncid, 'lhflx', varid   )
        1141           0 :    if ( status /= nf90_noerr ) then
        1142           0 :       status = nf90_inq_varid( ncid, 'lh', varid   )
        1143           0 :       if ( status /= nf90_noerr ) then
        1144           0 :         have_lhflx = .false.
        1145             :       else
        1146           0 :          call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
        1147           0 :          call wrap_get_vara_realx (ncid,varid,strt4,cnt4,lhflxobs)
        1148           0 :          have_lhflx = .true.
        1149             :       endif
        1150             :    else
        1151           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
        1152           0 :       call wrap_get_vara_realx (ncid,varid,strt4,cnt4,lhflxobs)
        1153           0 :       have_lhflx = .true.
        1154             :    endif
        1155             :
        1156           0 :    status = nf90_inq_varid( ncid, 'shflx', varid   )
        1157           0 :    if ( status /= nf90_noerr ) then
        1158           0 :       status = nf90_inq_varid( ncid, 'sh', varid   )
        1159           0 :       if ( status /= nf90_noerr ) then
        1160           0 :         have_shflx = .false.
        1161             :       else
        1162           0 :          call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
        1163           0 :         call wrap_get_vara_realx (ncid,varid,strt4,cnt4,shflxobs)
        1164           0 :         have_shflx = .true.
        1165             :       endif
        1166             :    else
        1167           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
        1168           0 :       call wrap_get_vara_realx (ncid,varid,strt4,cnt4,shflxobs)
        1169           0 :       have_shflx = .true.
        1170             :    endif
        1171             :
        1172             :    ! If REPLAY is used, then need to read in the global
        1173             :    !   energy fixer
        1174           0 :    status = nf90_inq_varid( ncid, 'heat_glob', varid   )
        1175           0 :    if (status /= nf90_noerr) then
        1176           0 :       have_heat_glob = .false.
        1177             :    else
        1178           0 :       call wrap_get_vara_realx (ncid,varid,strt4,cnt4,heat_glob_scm)
        1179           0 :       have_heat_glob = .true.
        1180             :    endif
        1181             :
        1182             : !
        1183             : !     fill in 3d forcing variables if we have both horizontal
        1184             : !     and vertical components, but not the 3d
        1185             : !
        1186           0 :    if ( .not. have_cnst(1) .and. have_divq .and. have_vertdivq ) then
        1187           0 :       do k=1,plev
        1188           0 :          do m=1,pcnst
        1189           0 :             divq3d(k,m) = divq(k,m) + vertdivq(k,m)
        1190             :          enddo
        1191             :       enddo
        1192           0 :       have_divq3d = .true.
        1193             :    endif
        1194             :
        1195           0 :    if ( .not. have_divt3d .and. have_divt .and. have_vertdivt ) then
        1196           0 :       if (masterproc) write(iulog,*) sub//'Don''t have divt3d - using divt and vertdivt'
        1197           0 :       do k=1,plev
        1198           0 :          divt3d(k) = divt(k) + vertdivt(k)
        1199             :       enddo
        1200           0 :       have_divt3d = .true.
        1201             :    endif
        1202             : !
        1203             : !     make sure that use_3dfrc flag is set to true if we only have
        1204             : !     3d forcing available
        1205             : !
        1206           0 :    if (scm_use_3dfrc) then
        1207           0 :       if (have_divt3d .and. have_divq3d) then
        1208           0 :          use_3dfrc = .true.
        1209             :       else
        1210           0 :          call endrun(sub//':ERROR :IOP file must have both divt3d and divq3d forcing when scm_use_3dfrc is set to .true.')
        1211             :       endif
        1212             :    endif
        1213             :
        1214           0 :    status =  nf90_inq_varid( ncid, 'beta', varid   )
        1215           0 :    if ( status /= nf90_noerr ) then
        1216           0 :       betacam = 0._r8
        1217             :    else
        1218           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
        1219           0 :       status = nf90_get_var(ncid, varid, srf(1), strt4)
        1220           0 :       betacam=srf(1)
        1221             :    endif
        1222             :
        1223           0 :    status =  nf90_inq_varid( ncid, 'fixmas', varid   )
        1224           0 :    if ( status /= nf90_noerr ) then
        1225           0 :       fixmascam=1.0_r8
        1226             :    else
        1227           0 :       call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
        1228           0 :       status = nf90_get_var(ncid, varid, srf(1), strt4)
        1229           0 :       fixmascam=srf(1)
        1230             :    endif
        1231             :
        1232           0 :    status = nf90_close( ncid )
        1233             :
        1234           0 :    deallocate(dplevs)
        1235             :
        1236           0 : end subroutine readiopdata
        1237             :
        1238           0 : subroutine setiopupdate
        1239             :
        1240             : !-----------------------------------------------------------------------
        1241             : !
        1242             : ! Open and read netCDF file to extract time information
        1243             : !
        1244             : !---------------------------Code history--------------------------------
        1245             : !
        1246             : ! Written by John Truesdale    August, 1996
        1247             : !
        1248             : !-----------------------------------------------------------------------
        1249             :   implicit none
        1250             :
        1251             :   character(len=*), parameter ::  sub = "setiopupdate"
        1252             :
        1253             : !------------------------------Locals-----------------------------------
        1254             :
        1255             :    integer :: next_date, next_sec
        1256             :    integer :: ncsec,ncdate                      ! current time of day,date
        1257             :    integer :: yr, mon, day                      ! year, month, and day component
        1258             : !------------------------------------------------------------------------------
        1259             :
        1260           0 :    call get_curr_date(yr,mon,day,ncsec)
        1261           0 :    ncdate=yr*10000 + mon*100 + day
        1262             :
        1263             : !------------------------------------------------------------------------------
        1264             : !     Check if iop data needs to be updated and set doiopupdate accordingly
        1265             : !------------------------------------------------------------------------------
        1266             :
        1267           0 :    if ( is_first_step() ) then
        1268           0 :       doiopupdate = .true.
        1269             :
        1270             :    else
        1271             :
        1272           0 :       call timemgr_time_inc(bdate, 0, next_date, next_sec, inc_s=tsec(iopTimeIdx+1))
        1273           0 :       if ( ncdate > next_date .or. (ncdate == next_date &
        1274             :          .and. ncsec >= next_sec)) then
        1275           0 :          doiopupdate = .true.
        1276             :          ! check to see if we need to move iopindex ahead more than 1 step
        1277           0 :          do while ( ncdate > next_date .or. (ncdate == next_date .and. ncsec >= next_sec))
        1278           0 :             iopTimeIdx = iopTimeIdx + 1
        1279           0 :             call timemgr_time_inc(bdate, 0, next_date, next_sec, inc_s=tsec(iopTimeIdx+1))
        1280             :          end do
        1281             : #if DEBUG > 2
        1282             :          if (masterproc) write(iulog,*) sub//'nstep = ',get_nstep()
        1283             :          if (masterproc) write(iulog,*) sub//'ncdate=',ncdate,' ncsec=',ncsec
        1284             :          if (masterproc) write(iulog,*) sub//'next_date=',next_date,' next_sec=',next_sec
        1285             :          if (masterproc) write(iulog,*) sub//':******* do iop update'
        1286             : #endif
        1287             :       else
        1288           0 :          doiopupdate = .false.
        1289             :       end if
        1290             :    endif                     ! if (endstep = 1 )
        1291             : !
        1292             : !     make sure we're
        1293             : !     not going past end of iop data
        1294             : !
        1295           0 :    if ( ncdate > last_date .or. (ncdate == last_date &
        1296             :       .and. ncsec > last_sec))  then
        1297           0 :       call endrun(sub//':ERROR: Reached the end of the time varient dataset')
        1298             :    endif
        1299             :
        1300             : #if DEBUG > 1
        1301             :    if (masterproc) write(iulog,*) sub//':iop time index = ' , ioptimeidx
        1302             : #endif
        1303             :
        1304           0 : end subroutine setiopupdate
        1305             :
        1306             : !===============================================================================
        1307             :
        1308           0 : subroutine plevs0 (nver, ps, ps0, hyam, hybm, hyai, hybi, pint    ,pmid    ,pdel)
        1309             :
        1310             : !-----------------------------------------------------------------------
        1311             : !
        1312             : ! Purpose:
        1313             : ! Define the pressures of the interfaces and midpoints from the
        1314             : ! coordinate definitions and the surface pressure.
        1315             : !
        1316             : ! Author: B. Boville
        1317             : !
        1318             : !-----------------------------------------------------------------------
        1319             :   implicit none
        1320             :
        1321             :
        1322             : !-----------------------------------------------------------------------
        1323             :   integer , intent(in)  :: nver         ! vertical dimension
        1324             :   real(r8), intent(in)  :: ps           ! Surface pressure (pascals)
        1325             :   real(r8), intent(in)  :: ps0          ! reference pressure (pascals)
        1326             :   real(r8), intent(in)  :: hyam(plev)   ! hybrid midpoint coef
        1327             :   real(r8), intent(in)  :: hybm(plev)   ! hybrid midpoint coef
        1328             :   real(r8), intent(in)  :: hyai(plevp)  ! hybrid interface coef
        1329             :   real(r8), intent(in)  :: hybi(plevp)  ! hybrid interface coef
        1330             :   real(r8), intent(out) :: pint(nver+1) ! Pressure at model interfaces
        1331             :   real(r8), intent(out) :: pmid(nver)   ! Pressure at model levels
        1332             :   real(r8), intent(out) :: pdel(nver)   ! Layer thickness (pint(k+1) - pint(k))
        1333             : !-----------------------------------------------------------------------
        1334             :
        1335             : !---------------------------Local workspace-----------------------------
        1336             :   integer :: k             ! Longitude, level indices
        1337             : !-----------------------------------------------------------------------
        1338             : !
        1339             : ! Set interface pressures
        1340             : !
        1341             : !$OMP PARALLEL DO PRIVATE (K)
        1342           0 :   do k=1,nver+1
        1343           0 :      pint(k) = hyai(k)*ps0 + hybi(k)*ps
        1344             :   end do
        1345             : !
        1346             : ! Set midpoint pressures and layer thicknesses
        1347             : !
        1348             : !$OMP PARALLEL DO PRIVATE (K)
        1349           0 :   do k=1,nver
        1350           0 :      pmid(k) = hyam(k)*ps0 + hybm(k)*ps
        1351           0 :      pdel(k) = pint(k+1) - pint(k)
        1352             :   end do
        1353             :
        1354           0 : end subroutine plevs0
        1355             :
        1356           0 : subroutine scmiop_flbc_inti ( co2vmr, ch4vmr, n2ovmr, f11vmr, f12vmr )
        1357             :   !-----------------------------------------------------------------------
        1358             :   !
        1359             :   ! Purpose:
        1360             :   ! Get start count for variable
        1361             :   !
        1362             :   !-----------------------------------------------------------------------
        1363             :
        1364             :   implicit none
        1365             :
        1366             :   real(r8), intent(out)  :: co2vmr, ch4vmr, n2ovmr, f11vmr, f12vmr
        1367             :
        1368             :   !-----------------------------------------------------------------------
        1369             :
        1370           0 :   co2vmr=co2vmrobs(1)
        1371           0 :   ch4vmr=ch4vmrobs(1)
        1372           0 :   n2ovmr=n2ovmrobs(1)
        1373           0 :   f11vmr=f11vmrobs(1)
        1374           0 :   f12vmr=f12vmrobs(1)
        1375           0 : end subroutine scmiop_flbc_inti
        1376             :
        1377           0 : subroutine get_start_count (ncid    ,varid  ,scmlat, scmlon, timeidx, start    ,count)
        1378             :
        1379             :   !-----------------------------------------------------------------------
        1380             :   !
        1381             :   ! Purpose:
        1382             :   ! set global lower boundary conditions
        1383             :   !
        1384             :   !-----------------------------------------------------------------------
        1385             :
        1386             :   implicit none
        1387             :
        1388             :   character(len=*), parameter ::  sub = "get_start_count"
        1389             :
        1390             : !-----------------------------------------------------------------------
        1391             :   integer , intent(in)    :: ncid         ! file id
        1392             :   integer , intent(in)    :: varid        ! variable id
        1393             :   integer , intent(in)    :: TimeIdx      ! time index
        1394             :   real(r8), intent(in)    :: scmlat,scmlon! scm lat/lon
        1395             :   integer , intent(out) :: start(:),count(:)
        1396             :
        1397             : !---------------------------Local workspace-----------------------------
        1398             :   integer            :: dims_set,nlev,var_ndims
        1399             :   logical            :: usable_var
        1400             :   character(len=cl)  :: dim_name
        1401             :   integer            :: var_dimIDs( NF90_MAX_VAR_DIMS )
        1402             :   real(r8)           :: closelat,closelon
        1403             :   integer            :: latidx,lonidx,status,i
        1404             : !-----------------------------------------------------------------------
        1405             :
        1406           0 :    call shr_scam_GetCloseLatLon(ncid,scmlat,scmlon,closelat,closelon,latidx,lonidx)
        1407             :
        1408           0 :    STATUS = NF90_INQUIRE_VARIABLE( NCID, varID, ndims=var_ndims )
        1409             : !
        1410             : !     surface variables
        1411             : !
        1412           0 :    if ( var_ndims == 0 ) then
        1413           0 :       call endrun(sub//':ERROR: var_ndims is 0 for varid:',varid)
        1414             :    endif
        1415             :
        1416           0 :    STATUS = NF90_INQUIRE_VARIABLE( NCID, varID, dimids=var_dimIDs)
        1417           0 :    if ( STATUS /= NF90_NOERR ) then
        1418           0 :       write(iulog,* ) sub//'ERROR - Cant get dimension IDs for varid', varid
        1419           0 :       call endrun(sub//':ERROR: Cant get dimension IDs for varid',varid)
        1420             :    endif
        1421             : !
        1422             : !     Initialize the start and count arrays
        1423             : !
        1424           0 :    dims_set = 0
        1425           0 :    nlev = 1
        1426           0 :    do i =  var_ndims, 1, -1
        1427             :
        1428           0 :       usable_var = .false.
        1429           0 :       STATUS = NF90_INQUIRE_DIMENSION( NCID, var_dimIDs( i ), dim_name )
        1430             :
        1431           0 :       if ( trim(dim_name) == 'lat' ) then
        1432           0 :          start( i ) =  latIdx
        1433           0 :          count( i ) = 1           ! Extract a single value
        1434           0 :          dims_set = dims_set + 1
        1435           0 :          usable_var = .true.
        1436             :       endif
        1437             :
        1438           0 :       if ( trim(dim_name) == 'lon' .or. trim(dim_name) == 'ncol' .or. trim(dim_name) == 'ncol_d' ) then
        1439           0 :          start( i ) = lonIdx
        1440           0 :          count( i ) = 1           ! Extract a single value
        1441           0 :          dims_set = dims_set + 1
        1442           0 :          usable_var = .true.
        1443             :       endif
        1444             :
        1445           0 :       if ( trim(dim_name) == 'lev' ) then
        1446           0 :          STATUS = NF90_INQUIRE_DIMENSION( NCID, var_dimIDs( i ), len=nlev )
        1447           0 :          start( i ) = 1
        1448           0 :          count( i ) = nlev       ! Extract all levels
        1449           0 :          dims_set = dims_set + 1
        1450           0 :          usable_var = .true.
        1451             :       endif
        1452             :
        1453           0 :       if ( trim(dim_name) == 'ilev' ) then
        1454           0 :          STATUS = NF90_INQUIRE_DIMENSION( NCID, var_dimIDs( i ), len=nlev )
        1455           0 :          start( i ) = 1
        1456           0 :          count( i ) = nlev        ! Extract all levels
        1457           0 :          dims_set = dims_set + 1
        1458           0 :          usable_var = .true.
        1459             :       endif
        1460             :
        1461           0 :       if ( trim(dim_name) == 'time' .OR. trim(dim_name) == 'tsec' ) then
        1462           0 :          start( i ) = TimeIdx
        1463           0 :          count( i ) = 1           ! Extract a single value
        1464           0 :          dims_set = dims_set + 1
        1465           0 :          usable_var = .true.
        1466             :       endif
        1467             :    end do
        1468           0 :  end subroutine get_start_count
        1469             :
        1470             : !=========================================================================
        1471           0 : subroutine setiopupdate_init
        1472             :
        1473             : !-----------------------------------------------------------------------
        1474             : !
        1475             : ! Open and read netCDF file to extract time information
        1476             : !   This subroutine should be called at the first SCM time step
        1477             : !
        1478             : !---------------------------Code history--------------------------------
        1479             : !
        1480             : ! Written by John Truesdale    August, 1996
        1481             : ! Modified for E3SM by  Peter Bogenschutz 2017 - onward
        1482             : !
        1483             : !-----------------------------------------------------------------------
        1484             :   implicit none
        1485             :
        1486             : !------------------------------Locals-----------------------------------
        1487             :
        1488             :    integer :: NCID,i
        1489             :    integer :: tsec_varID, time_dimID
        1490             :    integer :: bdate_varID
        1491             :    integer :: STATUS
        1492             :    integer :: next_date, next_sec
        1493             :    integer :: ncsec,ncdate                      ! current time of day,date
        1494             :    integer :: yr, mon, day                      ! year, month, and day component
        1495             :    integer :: start_ymd,start_tod
        1496             :
        1497             :    character(len=*), parameter ::  sub = "setiopupdate_init"
        1498             : !!------------------------------------------------------------------------------
        1499             :
        1500             :    ! Open and read pertinent information from the IOP file
        1501             :
        1502             :    call handle_ncerr( nf90_open (iopfile, 0, ncid),&
        1503           0 :         'ERROR - scamMod.F90:setiopupdate_init Failed to open iop file', __LINE__)
        1504             :
        1505             :    ! Read time (tsec) variable
        1506             :
        1507           0 :     STATUS = NF90_INQ_VARID( NCID, 'tsec', tsec_varID )
        1508           0 :     if ( STATUS /= NF90_NOERR ) then
        1509           0 :        write(iulog,*)sub//':ERROR: Cant get variable ID for tsec'
        1510           0 :        STATUS = NF90_CLOSE ( NCID )
        1511           0 :        call endrun(sub//':ERROR: Cant get variable ID for tsec')
        1512             :     end if
        1513             :
        1514           0 :     STATUS = NF90_INQ_VARID( NCID, 'bdate', bdate_varID )
        1515           0 :     if ( STATUS /= NF90_NOERR ) then
        1516           0 :        STATUS = NF90_INQ_VARID( NCID, 'basedate', bdate_varID )
        1517           0 :        if ( STATUS /= NF90_NOERR ) then
        1518           0 :           write(iulog,*)'ERROR - setiopupdate:Cant get variable ID for base date'
        1519           0 :           STATUS = NF90_CLOSE ( NCID )
        1520           0 :           call endrun(sub//':ERROR: Cant get variable ID for base date')
        1521             :        endif
        1522             :     endif
        1523             :
        1524           0 :     STATUS = NF90_INQ_DIMID( NCID, 'time', time_dimID )
        1525           0 :     if ( STATUS /= NF90_NOERR )  then
        1526           0 :        STATUS = NF90_INQ_DIMID( NCID, 'tsec', time_dimID )
        1527           0 :        if ( STATUS /= NF90_NOERR )  then
        1528           0 :           write(iulog,* )'ERROR - setiopupdate:Could not find variable dim ID for time'
        1529           0 :           STATUS = NF90_CLOSE ( NCID )
        1530           0 :           call endrun(sub//':ERROR:Could not find variable dim ID for time')
        1531             :        end if
        1532             :     end if
        1533             :
        1534           0 :     if ( STATUS /= NF90_NOERR )  &
        1535           0 :        write(iulog,*)'ERROR - setiopupdate:Cant get variable dim ID for time'
        1536             :
        1537           0 :     STATUS = NF90_INQUIRE_DIMENSION( NCID, time_dimID, len=ntime )
        1538           0 :     if ( STATUS /= NF90_NOERR )then
        1539           0 :        write(iulog,*)'ERROR - setiopupdate:Cant get time dimlen'
        1540             :     endif
        1541             :
        1542           0 :     if (.not.allocated(tsec)) allocate(tsec(ntime))
        1543             :
        1544           0 :     STATUS = NF90_GET_VAR( NCID, tsec_varID, tsec )
        1545           0 :     if ( STATUS /= NF90_NOERR )then
        1546           0 :        write(iulog,*)'ERROR - setiopupdate:Cant get variable tsec'
        1547             :     endif
        1548           0 :     STATUS = NF90_GET_VAR( NCID, bdate_varID, bdate )
        1549           0 :     if ( STATUS /= NF90_NOERR )then
        1550           0 :        write(iulog,*)'ERROR - setiopupdate:Cant get variable bdate'
        1551             :     endif
        1552             :
        1553             :     ! Close the netCDF file
        1554           0 :     STATUS = NF90_CLOSE( NCID )
        1555             :
        1556             :     ! determine the last date in the iop dataset
        1557             :
        1558           0 :     call timemgr_time_inc(bdate, 0, last_date, last_sec, inc_s=tsec(ntime))
        1559             :
        1560             :     ! set the iop dataset index
        1561           0 :     iopTimeIdx=0
        1562           0 :     do i=1,ntime           ! set the first ioptimeidx
        1563           0 :        call timemgr_time_inc(bdate, 0, next_date, next_sec, inc_s=tsec(i))
        1564           0 :        call get_start_date(yr,mon,day,start_tod)
        1565           0 :        start_ymd = yr*10000 + mon*100 + day
        1566             :
        1567           0 :        if ( start_ymd > next_date .or. (start_ymd == next_date &
        1568           0 :           .and. start_tod >= next_sec)) then
        1569           0 :           iopTimeIdx = i
        1570             :        endif
        1571             :     enddo
        1572             :
        1573           0 :     call get_curr_date(yr,mon,day,ncsec)
        1574           0 :     ncdate=yr*10000 + mon*100 + day
        1575             :
        1576           0 :     if (iopTimeIdx == 0.or.iopTimeIdx >= ntime) then
        1577           0 :        call timemgr_time_inc(bdate, 0, next_date, next_sec, inc_s=tsec(1))
        1578           0 :        write(iulog,*) 'Error::setiopupdate: Current model time does not fall within IOP period'
        1579           0 :        write(iulog,*) ' Current CAM Date is ',ncdate,' and ',ncsec,' seconds'
        1580           0 :        write(iulog,*) ' IOP start is        ',next_date,' and ',next_sec,'seconds'
        1581           0 :        write(iulog,*) ' IOP end is          ',last_date,' and ',last_sec,'seconds'
        1582           0 :        call endrun(sub//':ERROR: Current model time does not fall within IOP period')
        1583             :     endif
        1584             :
        1585           0 :     doiopupdate = .true.
        1586             :
        1587           0 : end subroutine setiopupdate_init
        1588             :
        1589             : end module scamMod

Generated by: LCOV version 1.14