LCOV - code coverage report
Current view: top level - dynamics/fv - metdata.F90 (source / functions) Hit Total Coverage
Test: coverage.info Lines: 0 1000 0.0 %
Date: 2025-03-14 01:23:43 Functions: 0 35 0.0 %

          Line data    Source code
       1             : module metdata
       2             : !-----------------------------------------------------------------------
       3             : !
       4             : ! BOP
       5             : !
       6             : ! !MODULE: metdata
       7             : !
       8             : ! !DESCRIPTION
       9             : ! Handles reading and interpolating offline meteorological data which
      10             : ! is used to drive the dynamics.
      11             : !
      12             : ! !USES
      13             :   use shr_kind_mod,       only: r8 => shr_kind_r8, r4 => shr_kind_r4
      14             :   use shr_cal_mod,        only: shr_cal_gregorian
      15             :   use time_manager,       only: get_curr_date, get_step_size, timemgr_is_caltype
      16             :   use spmd_utils,         only: masterproc
      17             :   use ppgrid,             only: pcols, pver, begchunk, endchunk
      18             :   use time_manager,       only: get_curr_calday, get_curr_date, get_step_size
      19             :   use cam_abortutils,     only: endrun
      20             :   use dynamics_vars,      only: T_FVDYCORE_GRID
      21             : 
      22             : #if ( defined SPMD )
      23             :   use mpishorthand,       only: mpicom, mpir8, mpiint,mpichar
      24             :   use mod_comm,           only: mp_sendirr,mp_recvirr
      25             : #endif
      26             :   use perf_mod
      27             :   use cam_logfile,        only: iulog
      28             :   use pio, only: file_desc_t, pio_put_att, pio_global, pio_get_att, pio_inq_att, &
      29             :        pio_inq_dimid, pio_inq_dimlen, pio_closefile, pio_get_var, pio_inq_varid, &
      30             :        pio_offset_kind
      31             :   use cam_pio_utils,      only: cam_pio_openfile
      32             : 
      33             : 
      34             :   implicit none
      35             : 
      36             :   private  ! all unless made public
      37             :   save
      38             : 
      39             : ! !PUBLIC MEMBERS
      40             : 
      41             :   public :: metdata_dyn_init  ! subroutine to open files, allocate blocked arrays, etc
      42             :   public :: metdata_phys_init ! subroutine to allocate chunked arrays
      43             :   public :: advance_met       ! subroutine to read more data and interpolate
      44             :   public :: get_met_fields    ! interface to set meteorology fields
      45             :   public :: get_met_srf1
      46             :   public :: get_met_srf2
      47             :   public :: get_us_vs
      48             :   public :: metdata_readnl
      49             :   public :: met_winds_on_walls
      50             :   public :: write_met_restart
      51             :   public :: read_met_restart
      52             :   public :: met_rlx
      53             :   public :: met_fix_mass
      54             :   public :: met_srf_feedback
      55             : 
      56             :   interface write_met_restart
      57             :      Module procedure write_met_restart_bin
      58             :      Module procedure write_met_restart_pio
      59             :   end interface
      60             : 
      61             :   interface read_met_restart
      62             :      Module procedure read_met_restart_bin
      63             :      Module procedure read_met_restart_pio
      64             :   end interface
      65             : 
      66             : 
      67             :   !------------------------------------------------------------------
      68             :   ! Interface to access the meteorology fields.  Possible invocations
      69             :   ! are as follows:
      70             :   !   call get_met_fields( physics_state, us, vs , tend, dt )
      71             :   !   call get_met_fields( u, v )
      72             :   !   call get_met_fields( cam_in_t )
      73             :   !------------------------------------------------------------------
      74             :   Interface get_met_fields                       ! overload accessors
      75             :      Module Procedure get_dyn_flds
      76             :      Module Procedure get_uv_centered
      77             :      Module Procedure get_ps
      78             :      Module Procedure get_ocn_ice_frcs
      79             :   End Interface
      80             : 
      81             :   real(r8), allocatable :: met_ps_next(:,:)   ! PS interpolated to next timestep
      82             :   real(r8), allocatable :: met_ps_curr(:,:)   ! PS interpolated to next timestep
      83             : 
      84             :   logical :: met_cell_wall_winds = .false.  ! true => met data winds are defined on model grid cell walls
      85             :   logical :: met_remove_file = .false.  ! delete metdata file when finished with it
      86             : 
      87             :   character(len=16) :: met_shflx_name = 'SHFLX'
      88             :   character(len=16) :: met_qflx_name = 'QFLX'
      89             :   real(r8) :: met_snowh_factor = 1._r8
      90             :   real(r8) :: met_shflx_factor = 1._r8
      91             :   real(r8) :: met_qflx_factor  = 1._r8
      92             :   logical  :: met_srf_feedback = .true.
      93             :   logical  :: met_srf_nudge_flux = .true. ! wsx, wsy, shf, and cflx nudged rather than forced.
      94             :                                           ! This is done primarily to prevent unrealistic
      95             :                                           ! surface temperatures.
      96             : 
      97             :   logical  :: met_srf_land = .true.       ! nudge surface fields over land (if false ignores
      98             :                                           ! all surface nudging for gridboxes with LANDFRAC=1)
      99             :   logical  :: met_srf_land_scale = .false. ! when met_srf_land is false, nudges proportional
     100             :                                           ! to the non land fraction, rather than just where
     101             :                                           ! LANDFRAC=1
     102             :   logical  :: met_srf_rad  = .false.      ! nudge albedo and lwup?
     103             :   logical  :: met_srf_refs = .false.      ! nudge 2m Q and T and 10m wind
     104             :   logical  :: met_srf_sst  = .false.      ! nudge sea surface temperature
     105             :   logical  :: met_srf_tau  = .true.       ! nudge taux and tauy
     106             :   logical  :: met_nudge_temp = .true.     ! nudge atmospheric temperature
     107             : 
     108             : 
     109             :   ! radiation/albedo surface field fill value (where there is no sunlight) read in from input data file
     110             :   real(r8) :: srf_fill_value
     111             : 
     112             : ! !REVISION HISTORY:
     113             : !   31 Oct 2003  Francis Vitt     Creation
     114             : !   05 Feb 2004  F Vitt   Removed reading/inperpolating PS for current timestep
     115             : !                         -- only met_ps_next is needed
     116             : !   10 Nov 2004  F Vitt   Implemented ability to read from series of files
     117             : !   16 Dec 2004  F Vitt   Added offline_met_defaultopts and offline_met_setopts
     118             : !   14 Jul 2005  W Sawyer Removed pmgrid, spmd_dyn dependencies
     119             : !   12 Apr 2006  W Sawyer Removed unneeded ghosting of met_us, met_vs
     120             : !   08 Apr 2010  J Edwards Replaced serial netcdf calls with pio interface
     121             : !
     122             : ! EOP
     123             : !-----------------------------------------------------------------------
     124             : ! $Id$
     125             : ! $Author$
     126             : !-----------------------------------------------------------------------
     127             : 
     128             :   type input2d
     129             :      real(r8), dimension(:,:), pointer :: data => null()
     130             :   endtype input2d
     131             : 
     132             :   type input3d
     133             :      real(r8), dimension(:,:,:), pointer :: data => null()
     134             :   endtype input3d
     135             : 
     136             :   real(r8), allocatable :: met_t(:,:,:)  ! interpolated temperature
     137             :   real(r8), allocatable :: met_u(:,:,:)  ! interpolated zonal wind
     138             :   real(r8), allocatable :: met_v(:,:,:)  ! interpolated meridional wind
     139             :   real(r8), allocatable :: met_us(:,:,:) ! interpolated zonal wind -staggered
     140             :   real(r8), allocatable :: met_vs(:,:,:) ! interpolated meridional wind -staggered
     141             :   real(r8), allocatable :: met_q(:,:,:)  ! interpolated water vapor
     142             : 
     143             :   real(r8), allocatable :: met_lhflx(:,:)! interpolated latent heat flux
     144             :   real(r8), allocatable :: met_shflx(:,:)! interpolated sensible heat flux
     145             :   real(r8), allocatable :: met_qflx(:,:) ! interpolated water vapor flux
     146             :   real(r8), allocatable :: met_taux(:,:) ! interpolated
     147             :   real(r8), allocatable :: met_tauy(:,:) ! interpolated
     148             :   real(r8), allocatable :: met_snowh(:,:) ! interpolated snow height
     149             : 
     150             :   real(r8), allocatable :: met_ts(:,:) ! interpolated
     151             : 
     152             :   real(r8), allocatable :: met_asdir(:,:) ! interpolated VIS direct albedo
     153             :   real(r8), allocatable :: met_asdif(:,:) ! interpolated VIS diffuse albedo
     154             :   real(r8), allocatable :: met_aldir(:,:) ! interpolated NIR direct albedo
     155             :   real(r8), allocatable :: met_aldif(:,:) ! interpolated NIR diffuse albedo
     156             :   real(r8), allocatable :: met_lwup(:,:)  ! interpolated upwelling LW flux
     157             :   real(r8), allocatable :: met_sst(:,:)   ! interpolated sea surface temperature
     158             :   real(r8), allocatable :: met_icefrac(:,:)  ! interpolated ice fraction
     159             :   real(r8), allocatable :: met_qref(:,:)  ! interpolated reference (2m) specific humidity
     160             :   real(r8), allocatable :: met_tref(:,:)  ! interpolated reference (2m) temperature
     161             :   real(r8), allocatable :: met_u10(:,:)   ! interpolated 10m wind speed
     162             : 
     163             :   type(input3d) :: met_ti(2)
     164             :   type(input3d) :: met_ui(2)
     165             :   type(input3d) :: met_vi(2)
     166             :   type(input3d) :: met_usi(2)
     167             :   type(input3d) :: met_vsi(2)
     168             :   type(input3d) :: met_qi(2)
     169             : 
     170             :   type(input2d) :: met_psi_next(2)
     171             :   type(input2d) :: met_psi_curr(2)
     172             :   type(input2d) :: met_lhflxi(2)
     173             :   type(input2d) :: met_shflxi(2)
     174             :   type(input2d) :: met_qflxi(2)
     175             :   type(input2d) :: met_tauxi(2)
     176             :   type(input2d) :: met_tauyi(2)
     177             :   type(input2d) :: met_tsi(2)
     178             :   type(input2d) :: met_snowhi(2)
     179             :   type(input2d) :: met_asdiri(2)
     180             :   type(input2d) :: met_asdifi(2)
     181             :   type(input2d) :: met_aldiri(2)
     182             :   type(input2d) :: met_aldifi(2)
     183             :   type(input2d) :: met_lwupi(2)
     184             :   type(input2d) :: met_ssti(2)
     185             :   type(input2d) :: met_icefraci(2)
     186             :   type(input2d) :: met_qrefi(2)
     187             :   type(input2d) :: met_trefi(2)
     188             :   type(input2d) :: met_u10i(2)
     189             : 
     190             :   integer :: dateid           ! var id of the date in the netCDF
     191             :   integer :: secid            ! var id of the sec data
     192             :   real(r8) :: datatimem = -1.e36_r8     ! time of prv. values read in
     193             :   real(r8) :: datatimep = -1.e36_r8     ! time of nxt. values read in
     194             :   real(r8) :: datatimemn = -1.e36_r8    ! time of prv. values read in for next timestep
     195             :   real(r8) :: datatimepn  = -1.e36_r8   ! time of nxt. values read in for next timestep
     196             : 
     197             :   integer, parameter :: nm=1    ! array index for previous (minus) data
     198             :   integer, parameter :: np=2    ! array indes for next (plus) data
     199             : 
     200             :   real(r8) :: curr_mod_time ! model time - calendar day
     201             :   real(r8) :: next_mod_time ! model time - calendar day - next time step
     202             : 
     203             :   character(len=256) :: curr_filename, next_filename, met_data_file
     204             :   character(len=256) :: met_filenames_list = ''
     205             :   character(len=256) :: met_data_path = ''
     206             :   type(file_desc_t) :: curr_fileid, next_fileid     ! the id of the NetCDF file
     207             :   real(r8), pointer, dimension(:) :: curr_data_times => null()
     208             :   real(r8), pointer, dimension(:) :: next_data_times => null()
     209             : 
     210             :   real(r8) :: alpha = 1.0_r8 ! don't read in water vapor
     211             :   !   real(r8), private :: alpha = 0.0  ! read in water vaper each time step
     212             : 
     213             :   real(r8), parameter ::  D0_0                    =     0.0_r8
     214             :   real(r8), parameter ::  D0_5                    =     0.5_r8
     215             :   real(r8), parameter ::  D0_75                   =     0.75_r8
     216             :   real(r8), parameter ::  D1_0                    =     1.0_r8
     217             :   real(r8), parameter ::  days_per_month          =    30.6_r8
     218             :   real(r8), parameter ::  days_per_non_leapyear   =   365.0_r8
     219             :   real(r8), parameter ::  days_per_year           =   365.25_r8
     220             :   real(r8), parameter ::  seconds_per_day         = 86400.0_r8
     221             :   real(r8), parameter ::  fill_value              = -9999.0_r8
     222             : 
     223             :   logical :: online_test = .false.
     224             :   logical :: debug = .false.
     225             : 
     226             :   real(r8) :: met_rlx(pver) = 0._r8
     227             :   integer  :: met_levels
     228             :   integer  :: num_met_levels
     229             : 
     230             :   real(r8) :: met_rlx_top = 60._r8
     231             :   real(r8) :: met_rlx_bot = 50._r8
     232             : 
     233             :   real(r8) :: met_rlx_bot_top = 0._r8
     234             :   real(r8) :: met_rlx_bot_bot = 0._r8
     235             : 
     236             :   real(r8) :: met_rlx_time = 0._r8
     237             : 
     238             : #if ( defined OFFLINE_DYN )
     239             :   logical  :: met_fix_mass = .true.
     240             : #else
     241             :   logical  :: met_fix_mass = .false.
     242             : #endif
     243             :   logical  :: has_ts = .false.
     244             :   logical  :: has_lhflx = .false.      ! Is LHFLX present in the met file?
     245             : 
     246             : contains
     247             : 
     248             : !-------------------------------------------------------------------------
     249             : ! meteorology data options
     250             : !-------------------------------------------------------------------------
     251           0 :   subroutine metdata_readnl(nlfile)
     252             : 
     253             :    use namelist_utils,  only: find_group_name
     254             :    use units,           only: getunit, freeunit
     255             :    use mpishorthand
     256             : 
     257             :    character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
     258             : 
     259             :    ! Local variables
     260             :    integer :: unitn, ierr
     261             :    character(len=*), parameter :: subname = 'metdata_readnl'
     262             : 
     263             :    namelist /metdata_nl/ &
     264             :         met_data_file, &
     265             :         met_data_path, &
     266             :         met_remove_file, &
     267             :         met_cell_wall_winds, &
     268             :         met_filenames_list, &
     269             :         met_rlx_top, &
     270             :         met_rlx_bot, &
     271             :         met_rlx_bot_top, &
     272             :         met_rlx_bot_bot, &
     273             :         met_rlx_time, &
     274             :         met_fix_mass, &
     275             :         met_shflx_name, &
     276             :         met_shflx_factor, &
     277             :         met_snowh_factor, &
     278             :         met_qflx_name, &
     279             :         met_qflx_factor, &
     280             :         met_srf_feedback, &
     281             :         met_srf_nudge_flux, &
     282             :         met_srf_land, &
     283             :         met_srf_land_scale, &
     284             :         met_srf_rad, &
     285             :         met_srf_refs, &
     286             :         met_srf_sst, &
     287             :         met_srf_tau, &
     288             :         met_nudge_temp
     289             : 
     290             :    ! Read namelist
     291           0 :    if (masterproc) then
     292           0 :       unitn = getunit()
     293           0 :       open( unitn, file=trim(nlfile), status='old' )
     294           0 :       call find_group_name(unitn, 'metdata_nl', status=ierr)
     295           0 :       if (ierr == 0) then
     296           0 :          read(unitn, metdata_nl, iostat=ierr)
     297           0 :          if (ierr /= 0) then
     298           0 :             call endrun(subname // ':: ERROR reading namelist')
     299             :          end if
     300             :       end if
     301           0 :       close(unitn)
     302           0 :       call freeunit(unitn)
     303             :    end if
     304             : 
     305             : #if ( defined SPMD )
     306             : 
     307             :    ! Broadcast namelist variables
     308             : 
     309           0 :    call mpibcast (met_data_file  ,len(met_data_file) ,mpichar,0,mpicom)
     310           0 :    call mpibcast (met_data_path  ,len(met_data_path) ,mpichar,0,mpicom)
     311           0 :    call mpibcast (met_remove_file    ,1 ,mpilog, 0, mpicom )
     312           0 :    call mpibcast (met_cell_wall_winds,1 ,mpilog, 0, mpicom )
     313           0 :    call mpibcast (met_filenames_list ,len(met_filenames_list),mpichar,0,mpicom)
     314           0 :    call mpibcast (met_rlx_top,        1 ,mpir8,  0, mpicom )
     315           0 :    call mpibcast (met_rlx_bot,        1 ,mpir8,  0, mpicom )
     316           0 :    call mpibcast (met_rlx_bot_top,    1 ,mpir8,  0, mpicom )
     317           0 :    call mpibcast (met_rlx_bot_bot,    1 ,mpir8,  0, mpicom )
     318           0 :    call mpibcast (met_rlx_time,       1 ,mpir8,  0, mpicom )
     319           0 :    call mpibcast (met_fix_mass,       1 ,mpilog, 0, mpicom )
     320           0 :    call mpibcast (met_qflx_name      ,len(met_qflx_name),     mpichar,0,mpicom)
     321           0 :    call mpibcast (met_shflx_name     ,len(met_shflx_name),    mpichar,0,mpicom)
     322           0 :    call mpibcast (met_qflx_factor    ,1, mpir8,  0, mpicom )
     323           0 :    call mpibcast (met_shflx_factor   ,1, mpir8,  0, mpicom )
     324           0 :    call mpibcast (met_snowh_factor   ,1, mpir8,  0, mpicom )
     325           0 :    call mpibcast (met_srf_feedback   ,1 ,mpilog, 0, mpicom )
     326           0 :    call mpibcast (met_srf_nudge_flux ,1 ,mpilog, 0, mpicom )
     327           0 :    call mpibcast (met_srf_land       ,1 ,mpilog, 0, mpicom )
     328           0 :    call mpibcast (met_srf_land_scale ,1 ,mpilog, 0, mpicom )
     329           0 :    call mpibcast (met_srf_rad        ,1 ,mpilog, 0, mpicom )
     330           0 :    call mpibcast (met_srf_refs       ,1 ,mpilog, 0, mpicom )
     331           0 :    call mpibcast (met_srf_sst        ,1 ,mpilog, 0, mpicom )
     332           0 :    call mpibcast (met_srf_tau        ,1 ,mpilog, 0, mpicom )
     333           0 :    call mpibcast (met_nudge_temp     ,1 ,mpilog, 0, mpicom )
     334             : #endif
     335             : 
     336           0 :    if (masterproc) then
     337           0 :        write(iulog,*)'Time-variant meteorological dataset (met_data_file) is: ', trim(met_data_file)
     338           0 :        write(iulog,*)'Meteorological data file will be removed (met_remove_file): ', met_remove_file
     339           0 :        write(iulog,*)'Meteorological winds are on cell walls (met_cell_wall_winds): ', met_cell_wall_winds
     340           0 :        write(iulog,*)'Meteorological file names list file: ', trim(met_filenames_list)
     341           0 :        write(iulog,*)'Meteorological relax ramp region top at top is (km): ', met_rlx_top
     342           0 :        write(iulog,*)'Meteorological relax ramp region bottom at top is (km): ', met_rlx_bot
     343           0 :        write(iulog,*)'Meteorological relax ramp region top at bottom is (km): ', met_rlx_bot_top
     344           0 :        write(iulog,*)'Meteorological relax ramp region bottom at bottom is (km): ', met_rlx_bot_bot
     345           0 :        write(iulog,*)'Meteorological relaxation time (hours): ',met_rlx_time
     346           0 :        write(iulog,*)'Offline driver mass fixer is trurned on (met_fix_mass): ',met_fix_mass
     347           0 :        write(iulog,*)'Meteorological shflx field name : ', trim(met_shflx_name)
     348           0 :        write(iulog,*)'Meteorological shflx multiplication factor : ', met_shflx_factor
     349           0 :        write(iulog,*)'Meteorological qflx field name : ', trim(met_qflx_name)
     350           0 :        write(iulog,*)'Meteorological qflx multiplication factor : ', met_qflx_factor
     351           0 :        write(iulog,*)'Meteorological snowh multiplication factor : ', met_snowh_factor
     352           0 :        write(iulog,*)'Meteorological allow srf models feedbacks : ', met_srf_feedback
     353           0 :        write(iulog,*)'Meteorological allow srf land nudging : ', met_srf_land
     354           0 :        write(iulog,*)'Meteorological scale srf land nudging : ', met_srf_land_scale
     355           0 :        write(iulog,*)'Meteorological allow srf radiation nudging : ', met_srf_rad
     356           0 :        write(iulog,*)'Meteorological allow srf reference field nudging : ', met_srf_refs
     357           0 :        write(iulog,*)'Meteorological allow srf sst nudging : ', met_srf_sst
     358           0 :        write(iulog,*)'Meteorological allow srf tau nudging : ', met_srf_tau
     359           0 :        write(iulog,*)'Meteorological allow atm tempature nudging : ',met_nudge_temp
     360             :     endif
     361             : 
     362           0 :  end subroutine metdata_readnl
     363             : 
     364             : !--------------------------------------------------------------------------
     365             : ! Opens file, allocates arrays
     366             : !--------------------------------------------------------------------------
     367           0 :   subroutine metdata_dyn_init(grid)
     368             :     use infnan,          only : nan, assignment(=)
     369             :     use cam_control_mod, only : restart_run
     370             :     implicit none
     371             : 
     372             : ! !INPUT PARAMETERS:
     373             :       type (T_FVDYCORE_GRID), intent(in) :: grid
     374             : 
     375             : 
     376             :     integer            :: im, km, jfirst, jlast, kfirst, klast
     377             :     integer            :: ng_d, ng_s
     378             : 
     379           0 :     im        = grid%im
     380           0 :     km        = grid%km
     381           0 :     jfirst    = grid%jfirst
     382           0 :     jlast     = grid%jlast
     383           0 :     kfirst    = grid%kfirst
     384           0 :     klast     = grid%klast
     385           0 :     ng_d      = grid%ng_d
     386           0 :     ng_s      = grid%ng_s
     387             : 
     388             : 
     389           0 :     if (.not. restart_run) then ! initial run or branch run
     390           0 :        curr_filename = met_data_file
     391           0 :        next_filename = ''
     392             :     else
     393             :        ! restart run
     394             :        ! curr_filename & next_filename already set by restart_dynamics
     395             :     endif
     396             : 
     397           0 :     call open_met_datafile( curr_filename, curr_fileid, curr_data_times, met_data_path, check_dims=.true., grid=grid)
     398             : 
     399           0 :     if ( len_trim(next_filename) > 0 ) &
     400           0 :          call open_met_datafile( next_filename, next_fileid, next_data_times, met_data_path )
     401             : 
     402             : !
     403             : ! allocate space for data arrays ...
     404             : !
     405             :     ! dynamics grid
     406             : 
     407           0 :     allocate( met_psi_next(nm)%data(im, jfirst:jlast) )
     408           0 :     allocate( met_psi_next(np)%data(im, jfirst:jlast) )
     409           0 :     allocate( met_psi_curr(nm)%data(im, jfirst:jlast) )
     410           0 :     allocate( met_psi_curr(np)%data(im, jfirst:jlast) )
     411           0 :     allocate( met_ps_next(im, jfirst:jlast) )
     412           0 :     allocate( met_ps_curr(im, jfirst:jlast) )
     413             : 
     414           0 :     allocate( met_us(im, jfirst-ng_d:jlast+ng_s, kfirst:klast) )
     415           0 :     allocate( met_vs(im, jfirst-ng_s:jlast+ng_d, kfirst:klast) )
     416             : 
     417           0 :     met_us = nan
     418           0 :     met_vs = nan
     419             : 
     420           0 :     if (met_cell_wall_winds) then
     421           0 :        allocate( met_usi(nm)%data(im, jfirst:jlast, kfirst:klast) )
     422           0 :        allocate( met_usi(np)%data(im, jfirst:jlast, kfirst:klast) )
     423           0 :        allocate( met_vsi(nm)%data(im, jfirst:jlast, kfirst:klast) )
     424           0 :        allocate( met_vsi(np)%data(im, jfirst:jlast, kfirst:klast) )
     425             :     endif
     426             : 
     427           0 :     if (.not. met_cell_wall_winds) then
     428             : 
     429           0 :        allocate( met_u(im, jfirst-ng_d:jlast+ng_d, kfirst:klast) )
     430           0 :        allocate( met_ui(nm)%data(im, jfirst:jlast, kfirst:klast) )
     431           0 :        allocate( met_ui(np)%data(im, jfirst:jlast, kfirst:klast) )
     432             : 
     433           0 :        allocate( met_v(im, jfirst-ng_s:jlast+ng_d, kfirst:klast) )
     434           0 :        allocate( met_vi(nm)%data(im, jfirst:jlast, kfirst:klast) )
     435           0 :        allocate( met_vi(np)%data(im, jfirst:jlast, kfirst:klast) )
     436             : 
     437             :     endif
     438             : 
     439           0 :  end subroutine metdata_dyn_init
     440             : 
     441             : !=================================================================================
     442             : 
     443           0 :  subroutine metdata_phys_init
     444             :    use infnan,      only : nan, assignment(=)
     445             :    use cam_history, only : addfld, horiz_only
     446             : 
     447           0 :    call addfld ('MET_RLX',     (/ 'lev' /),  'A', ' ', 'Meteorology relax function',              gridname='fv_centers')
     448           0 :    call addfld ('MET_TAUX',    horiz_only,   'A', ' ', 'Meteorology taux',                        gridname='physgrid')
     449           0 :    call addfld ('MET_TAUY',    horiz_only,   'A', ' ', 'Meteorology tauy',                        gridname='physgrid')
     450           0 :    call addfld ('MET_LHFX',    horiz_only,   'A', ' ', 'Meteorology lhflx',                       gridname='physgrid')
     451           0 :    call addfld ('MET_SHFX',    horiz_only,   'A', ' ', 'Meteorology shflx',                       gridname='physgrid')
     452           0 :    call addfld ('MET_QFLX',    horiz_only,   'A', ' ', 'Meteorology qflx',                        gridname='physgrid')
     453           0 :    call addfld ('MET_PS',      horiz_only,   'A', ' ', 'Meteorology PS',                          gridname='fv_centers')
     454           0 :    call addfld ('MET_T',       (/ 'lev' /),  'A', ' ', 'Meteorology T',                           gridname='physgrid')
     455           0 :    call addfld ('MET_U',       (/ 'lev' /),  'A', ' ', 'Meteorology U',                           gridname='fv_centers')
     456           0 :    call addfld ('MET_V',       (/ 'lev' /),  'A', ' ', 'Meteorology V',                           gridname='fv_centers')
     457           0 :    call addfld ('MET_SNOWH',   horiz_only,   'A', ' ', 'Meteorology snow height',                 gridname='physgrid')
     458             : 
     459           0 :    call addfld ('MET_TS',      horiz_only,   'A', 'K', 'Meteorology TS',                          gridname='physgrid')
     460           0 :    call addfld ('MET_OCNFRC',  horiz_only,   'A', 'fraction', 'Ocean frac derived from met TS',   gridname='physgrid')
     461           0 :    call addfld ('MET_ICEFRC',  horiz_only,   'A', 'fraction', 'Sea ice frac derived from met TS', gridname='physgrid')
     462             : 
     463           0 :    call addfld ('MET_ASDIR',   horiz_only,   'A', '1', 'Meteorology ASDIR',                       gridname='physgrid')
     464           0 :    call addfld ('MET_ASDIF',   horiz_only,   'A', '1', 'Meteorology ASDIF',                       gridname='physgrid')
     465           0 :    call addfld ('MET_ALDIR',   horiz_only,   'A', '1', 'Meteorology ALDIR',                       gridname='physgrid')
     466           0 :    call addfld ('MET_ALDIF',   horiz_only,   'A', '1', 'Meteorology ALDIF',                       gridname='physgrid')
     467           0 :    call addfld ('MET_LWUP',    horiz_only,   'A', 'Wm-2', 'Meteorology LWUP',                     gridname='physgrid')
     468           0 :    call addfld ('MET_SST',     horiz_only,   'A', 'K', 'Meteorology SST',                         gridname='physgrid')
     469           0 :    call addfld ('MET_ICEFRAC', horiz_only,   'A', 'K', 'Meteorology ICEFRAC',                     gridname='physgrid')
     470           0 :    call addfld ('MET_QREF',    horiz_only,   'A', 'kg/kg', 'Meteorology QREF',                    gridname='physgrid')
     471           0 :    call addfld ('MET_TREF',    horiz_only,   'A', 'K', 'Meteorology TREF',                        gridname='physgrid')
     472           0 :    call addfld ('MET_U10',     horiz_only,   'A', 'ms-1', 'Meteorology U10',                      gridname='physgrid')
     473             : 
     474             : ! allocate chunked arrays
     475             : 
     476           0 :    allocate( met_ti(nm)%data(pcols,pver,begchunk:endchunk) )
     477           0 :    allocate( met_ti(np)%data(pcols,pver,begchunk:endchunk) )
     478           0 :    allocate( met_t(pcols,pver,begchunk:endchunk) )
     479             : 
     480           0 :    allocate( met_qi(nm)%data(pcols,pver,begchunk:endchunk) )
     481           0 :    allocate( met_qi(np)%data(pcols,pver,begchunk:endchunk) )
     482           0 :    allocate( met_q(pcols,pver,begchunk:endchunk) )
     483             : 
     484           0 :    allocate( met_lhflxi(nm)%data(pcols,begchunk:endchunk) )
     485           0 :    allocate( met_lhflxi(np)%data(pcols,begchunk:endchunk) )
     486           0 :    allocate( met_lhflx(pcols,begchunk:endchunk) )
     487             : 
     488           0 :    allocate( met_shflxi(nm)%data(pcols,begchunk:endchunk) )
     489           0 :    allocate( met_shflxi(np)%data(pcols,begchunk:endchunk) )
     490           0 :    allocate( met_shflx(pcols,begchunk:endchunk) )
     491             : 
     492           0 :    allocate( met_qflxi(nm)%data(pcols,begchunk:endchunk) )
     493           0 :    allocate( met_qflxi(np)%data(pcols,begchunk:endchunk) )
     494           0 :    allocate( met_qflx(pcols,begchunk:endchunk) )
     495             : 
     496           0 :    allocate( met_tauxi(nm)%data(pcols,begchunk:endchunk) )
     497           0 :    allocate( met_tauxi(np)%data(pcols,begchunk:endchunk) )
     498           0 :    allocate( met_taux(pcols,begchunk:endchunk) )
     499             : 
     500           0 :    allocate( met_tauyi(nm)%data(pcols,begchunk:endchunk) )
     501           0 :    allocate( met_tauyi(np)%data(pcols,begchunk:endchunk) )
     502           0 :    allocate( met_tauy(pcols,begchunk:endchunk) )
     503             : 
     504           0 :    allocate( met_tsi(nm)%data(pcols,begchunk:endchunk) )
     505           0 :    allocate( met_tsi(np)%data(pcols,begchunk:endchunk) )
     506           0 :    allocate( met_ts(pcols,begchunk:endchunk) )
     507           0 :    met_ts(:,:) = nan
     508             : 
     509           0 :    if(.not.met_srf_feedback) then
     510           0 :       allocate( met_snowhi(nm)%data(pcols,begchunk:endchunk) )
     511           0 :       allocate( met_snowhi(np)%data(pcols,begchunk:endchunk) )
     512           0 :       allocate( met_snowh(pcols,begchunk:endchunk) )
     513           0 :       met_snowh(:,:) = nan
     514             :    endif
     515             : 
     516           0 :    if(met_srf_rad) then
     517           0 :       allocate( met_asdiri(nm)%data(pcols,begchunk:endchunk) )
     518           0 :       allocate( met_asdiri(np)%data(pcols,begchunk:endchunk) )
     519           0 :       allocate( met_asdir(pcols,begchunk:endchunk) )
     520             : 
     521           0 :       allocate( met_asdifi(nm)%data(pcols,begchunk:endchunk) )
     522           0 :       allocate( met_asdifi(np)%data(pcols,begchunk:endchunk) )
     523           0 :       allocate( met_asdif(pcols,begchunk:endchunk) )
     524             : 
     525           0 :       allocate( met_aldiri(nm)%data(pcols,begchunk:endchunk) )
     526           0 :       allocate( met_aldiri(np)%data(pcols,begchunk:endchunk) )
     527           0 :       allocate( met_aldir(pcols,begchunk:endchunk) )
     528             : 
     529           0 :       allocate( met_aldifi(nm)%data(pcols,begchunk:endchunk) )
     530           0 :       allocate( met_aldifi(np)%data(pcols,begchunk:endchunk) )
     531           0 :       allocate( met_aldif(pcols,begchunk:endchunk) )
     532             : 
     533           0 :       allocate( met_lwupi(nm)%data(pcols,begchunk:endchunk) )
     534           0 :       allocate( met_lwupi(np)%data(pcols,begchunk:endchunk) )
     535           0 :       allocate( met_lwup(pcols,begchunk:endchunk) )
     536             :    end if
     537             : 
     538           0 :    if(met_srf_refs) then
     539           0 :       allocate( met_qrefi(nm)%data(pcols,begchunk:endchunk) )
     540           0 :       allocate( met_qrefi(np)%data(pcols,begchunk:endchunk) )
     541           0 :       allocate( met_qref(pcols,begchunk:endchunk) )
     542             : 
     543           0 :       allocate( met_trefi(nm)%data(pcols,begchunk:endchunk) )
     544           0 :       allocate( met_trefi(np)%data(pcols,begchunk:endchunk) )
     545           0 :       allocate( met_tref(pcols,begchunk:endchunk) )
     546             : 
     547           0 :       allocate( met_u10i(nm)%data(pcols,begchunk:endchunk) )
     548           0 :       allocate( met_u10i(np)%data(pcols,begchunk:endchunk) )
     549           0 :       allocate( met_u10(pcols,begchunk:endchunk) )
     550             :    end if
     551             : 
     552           0 :    if(met_srf_sst) then
     553           0 :       allocate( met_ssti(nm)%data(pcols,begchunk:endchunk) )
     554           0 :       allocate( met_ssti(np)%data(pcols,begchunk:endchunk) )
     555           0 :       allocate( met_sst(pcols,begchunk:endchunk) )
     556             : 
     557           0 :       allocate( met_icefraci(nm)%data(pcols,begchunk:endchunk) )
     558           0 :       allocate( met_icefraci(np)%data(pcols,begchunk:endchunk) )
     559           0 :       allocate( met_icefrac(pcols,begchunk:endchunk) )
     560             :     end if
     561             : 
     562           0 :    call set_met_rlx()
     563             : 
     564             :    ! initialize phys surface fields...
     565           0 :    call get_model_time()
     566           0 :    call check_files()
     567           0 :    call read_phys_srf_flds()
     568           0 :    call interp_phys_srf_flds()
     569           0 :    datatimem = -1.e36_r8
     570           0 :    datatimep = -1.e36_r8
     571           0 :  end subroutine metdata_phys_init
     572             : 
     573             : 
     574             : !-----------------------------------------------------------------------
     575             : ! Reads more data if needed and interpolates data to current model time
     576             : !-----------------------------------------------------------------------
     577           0 :  subroutine advance_met(grid)
     578           0 :    use cam_history, only : outfld
     579             :    implicit none
     580             : 
     581             :    type (T_FVDYCORE_GRID), intent(in) :: grid
     582             : 
     583           0 :    real(r8) ::  met_rlx_2d(grid%ifirstxy:grid%ilastxy,grid%km)
     584             :    integer :: i,j,k, idim
     585             : 
     586           0 :    call t_startf('MET__advance')
     587             : 
     588             : !
     589             : !
     590           0 :     call get_model_time()
     591             : 
     592           0 :     if ( ( curr_mod_time > datatimep ) .or. &
     593             :          ( next_mod_time > datatimepn ) ) then
     594           0 :        call check_files()
     595             :     endif
     596             : 
     597           0 :     if ( curr_mod_time > datatimep ) then
     598           0 :        call read_next_metdata(grid)
     599             :     end if
     600             : 
     601           0 :     if ( next_mod_time > datatimepn ) then
     602           0 :        call read_next_ps(grid)
     603             :     end if
     604             : 
     605             : ! need to inperpolate the data, regardless !
     606             : ! each mpi tasks needs to interpolate
     607           0 :     call interpolate_metdata(grid)
     608             : 
     609           0 :     call t_stopf('MET__advance')
     610             : 
     611           0 :     idim = grid%ilastxy - grid%ifirstxy + 1
     612           0 :     do j = grid%jfirstxy, grid%jlastxy
     613           0 :        do k = 1, grid%km
     614           0 :           do i = grid%ifirstxy, grid%ilastxy
     615           0 :              met_rlx_2d(i,k) = met_rlx(k)
     616             :           enddo
     617             :        enddo
     618           0 :        call outfld('MET_RLX',met_rlx_2d, idim, j)
     619             :     enddo
     620           0 :   end subroutine advance_met
     621             : 
     622             : !-------------------------------------------------------------------
     623             : ! Method to get some the meteorology data.
     624             : ! Sets the following cam_in_t member fields to the
     625             : ! meteorology data :
     626             : !   qflx
     627             : !   lhflx
     628             : !   shflx
     629             : !   taux
     630             : !   tauy
     631             : !   snowh
     632             : !-------------------------------------------------------------------
     633           0 :   subroutine get_met_srf2( cam_in )
     634           0 :     use camsrfexch,          only: cam_in_t
     635             :     use phys_grid,           only: get_ncols_p
     636             :     use cam_history,         only: outfld
     637             :     use shr_const_mod,       only: shr_const_stebol
     638             :     use physconst,           only: latvap, latice
     639             : 
     640             :     implicit none
     641             : 
     642             :     type(cam_in_t), intent(inout), dimension(begchunk:endchunk) :: cam_in
     643             : 
     644             :     integer :: c,ncol,i
     645             :     real(r8) :: met_rlx_sfc(pcols)
     646             :     real(r8) :: lcl_rlx(pcols)
     647             : 
     648           0 :     do c=begchunk,endchunk
     649           0 :        ncol = get_ncols_p(c)
     650             : 
     651             :        ! Nudge or force the surface fields?
     652           0 :        if (met_srf_nudge_flux) then
     653           0 :          met_rlx_sfc(:ncol) = met_rlx(pver)
     654             :        else
     655           0 :          met_rlx_sfc(:ncol) = 1._r8
     656             :        end if
     657             : 
     658             :        ! Don't nudge the surface for locations that are entirely land?
     659           0 :        if (.not. met_srf_land) then
     660             : 
     661             :           ! Nudging land and forcing ocean.
     662           0 :           if (met_srf_land_scale) then
     663           0 :              met_rlx_sfc(:ncol) = (1._r8 - cam_in(c)%landfrac(:ncol)) * &
     664             :                                   met_rlx_sfc(:ncol) + &
     665           0 :                                   cam_in(c)%landfrac(:ncol) * met_rlx(pver)
     666             :           else
     667           0 :             where(cam_in(c)%landfrac(:ncol) == 1._r8) met_rlx_sfc(:ncol) = 0._r8
     668             :           end if
     669             :        end if
     670             : 
     671           0 :        if (met_srf_tau) then
     672           0 :          cam_in(c)%wsx(:ncol)     = (1._r8-met_rlx_sfc(:ncol)) * cam_in(c)%wsx(:ncol)    + met_rlx_sfc(:ncol) * met_taux(:ncol,c)
     673           0 :          cam_in(c)%wsy(:ncol)     = (1._r8-met_rlx_sfc(:ncol)) * cam_in(c)%wsy(:ncol)    + met_rlx_sfc(:ncol) * met_tauy(:ncol,c)
     674             :        end if
     675             : 
     676           0 :        cam_in(c)%shf(:ncol)     = (1._r8-met_rlx_sfc(:ncol)) * cam_in(c)%shf(:ncol)    + &
     677           0 :              met_rlx_sfc(:ncol) * (met_shflx(:ncol,c) * met_shflx_factor)
     678             :        cam_in(c)%cflx(:ncol,1)  = (1._r8-met_rlx_sfc(:ncol)) * cam_in(c)%cflx(:ncol,1) + &
     679           0 :              met_rlx_sfc(:ncol) * (met_qflx(:ncol,c)  * met_qflx_factor)
     680             : 
     681             :        ! If present, nudge the latent heat; otherwise, make an approximation by scaling the
     682             :        ! water vapor flux.
     683           0 :        if (has_lhflx) then
     684             :          cam_in(c)%lhf(:ncol)     = (1._r8-met_rlx_sfc(:ncol)) * cam_in(c)%lhf(:ncol)    + &
     685           0 :                met_rlx_sfc(:ncol) * (met_lhflx(:ncol,c) * met_qflx_factor)
     686             :        else
     687             :          cam_in(c)%lhf(:ncol)     = (1._r8-met_rlx_sfc(:ncol)) * cam_in(c)%lhf(:ncol)    + &
     688           0 :                met_rlx_sfc(:ncol) * (met_qflx(:ncol,c) * met_qflx_factor * latvap)
     689             :        end if
     690             : 
     691           0 :        if (met_srf_rad) then
     692             : 
     693             :           ! There can be fill values in the albedos from the met fields, so use the cam_in value
     694             :           ! if fill is detected. This could be jumpy if that value gets used, but should be in
     695             :           ! an area with no downwelling solar. Time interpolate around the terminator could cause
     696             :           ! problems, but the interpolation provides a non-fill value if either endpoint of the
     697             :           ! interpolation is not fill.
     698           0 :           where(met_asdir(:ncol,c) == srf_fill_value)
     699             :              lcl_rlx(:ncol) = 0._r8
     700             :           elsewhere
     701             :              lcl_rlx(:ncol) = met_rlx_sfc(:ncol)
     702             :           end where
     703           0 :           cam_in(c)%asdir(:ncol) = (1._r8-lcl_rlx(:ncol)) * cam_in(c)%asdir(:ncol) + lcl_rlx(:ncol) * met_asdir(:ncol,c)
     704             : 
     705           0 :           where(met_asdif(:ncol,c) == srf_fill_value)
     706             :              lcl_rlx(:ncol) = 0._r8
     707             :           elsewhere
     708             :              lcl_rlx(:ncol) = met_rlx_sfc(:ncol)
     709             :           end where
     710           0 :           cam_in(c)%asdif(:ncol) = (1._r8-lcl_rlx(:ncol)) * cam_in(c)%asdif(:ncol) + lcl_rlx(:ncol) * met_asdif(:ncol,c)
     711             : 
     712           0 :           where(met_aldir(:ncol,c) == srf_fill_value)
     713             :              lcl_rlx(:ncol) = 0._r8
     714             :           elsewhere
     715             :              lcl_rlx(:ncol) = met_rlx_sfc(:ncol)
     716             :           end where
     717           0 :           cam_in(c)%aldir(:ncol) = (1._r8-lcl_rlx(:ncol)) * cam_in(c)%aldir(:ncol) + lcl_rlx(:ncol) * met_aldir(:ncol,c)
     718             : 
     719           0 :           where(met_aldif(:ncol,c) == srf_fill_value)
     720             :              lcl_rlx(:ncol) = 0._r8
     721             :           elsewhere
     722             :              lcl_rlx(:ncol) = met_rlx_sfc(:ncol)
     723             :           end where
     724           0 :           cam_in(c)%aldif(:ncol) = (1._r8-lcl_rlx(:ncol)) * cam_in(c)%aldif(:ncol) + lcl_rlx(:ncol) * met_aldif(:ncol,c)
     725             : 
     726           0 :           cam_in(c)%lwup(:ncol) = (1._r8-met_rlx_sfc(:ncol)) * cam_in(c)%lwup(:ncol) + met_rlx_sfc(:ncol) * met_lwup(:ncol,c)
     727             :        end if
     728             : 
     729           0 :        if (met_srf_refs) then
     730           0 :           cam_in(c)%qref(:ncol)  = (1._r8-met_rlx_sfc(:ncol)) * cam_in(c)%qref(:ncol)    + met_rlx_sfc(:ncol) * met_qref(:ncol,c)
     731           0 :           cam_in(c)%tref(:ncol)  = (1._r8-met_rlx_sfc(:ncol)) * cam_in(c)%tref(:ncol)    + met_rlx_sfc(:ncol) * met_tref(:ncol,c)
     732           0 :           cam_in(c)%u10(:ncol)   = (1._r8-met_rlx_sfc(:ncol)) * cam_in(c)%u10(:ncol)     + met_rlx_sfc(:ncol) * met_u10(:ncol,c)
     733             :        end if
     734             : 
     735           0 :        if (met_srf_sst) then
     736             : 
     737             :           ! Meteorological sst is 0 over 100% land, so use the cam_in value if the meteorology thinks
     738             :           ! it is land.
     739           0 :           where(met_sst(:ncol,c) == srf_fill_value)
     740             :              lcl_rlx(:ncol) = 0._r8
     741             :           elsewhere
     742             :              lcl_rlx(:ncol) = met_rlx_sfc(:ncol)
     743             :           end where
     744           0 :           cam_in(c)%sst(:ncol) = (1._r8-lcl_rlx(:ncol)) * cam_in(c)%sst(:ncol) + lcl_rlx(:ncol) * met_sst(:ncol,c)
     745             : 
     746           0 :           where(met_icefrac(:ncol,c) == srf_fill_value)
     747             :              lcl_rlx(:ncol) = 0._r8
     748             :           elsewhere
     749             :              lcl_rlx(:ncol) = met_rlx_sfc(:ncol)
     750             :           end where
     751           0 :           cam_in(c)%icefrac(:ncol) = (1._r8-lcl_rlx(:ncol)) * cam_in(c)%icefrac(:ncol) + lcl_rlx(:ncol) * met_icefrac(:ncol,c)
     752             : 
     753             :        end if
     754             :      end do                    ! Chunk loop
     755             : 
     756           0 :     if (debug) then
     757           0 :        if (masterproc) then
     758           0 :           write(iulog,*)'METDATA   maxval(met_taux),minval(met_taux): ',maxval(met_taux),minval(met_taux)
     759           0 :           write(iulog,*)'METDATA   maxval(met_tauy),minval(met_tauy): ',maxval(met_tauy),minval(met_tauy)
     760           0 :           write(iulog,*)'METDATA   maxval(met_lhflx),minval(met_lhflx): ',maxval(met_lhflx),minval(met_lhflx)
     761           0 :           write(iulog,*)'METDATA   maxval(met_shflx),minval(met_shflx): ',maxval(met_shflx),minval(met_shflx)
     762           0 :           write(iulog,*)'METDATA   maxval(met_qflx),minval(met_qflx): ',maxval(met_qflx),minval(met_qflx)
     763           0 :           write(iulog,*)'METDATA   maxval(met_asdir),minval(met_asdir): ',maxval(met_asdir),minval(met_asdir)
     764           0 :           write(iulog,*)'METDATA   maxval(met_asdif),minval(met_asdif): ',maxval(met_asdif),minval(met_asdif)
     765           0 :           write(iulog,*)'METDATA   maxval(met_aldir),minval(met_aldir): ',maxval(met_aldir),minval(met_aldir)
     766           0 :           write(iulog,*)'METDATA   maxval(met_aldif),minval(met_aldif): ',maxval(met_aldif),minval(met_aldif)
     767           0 :           write(iulog,*)'METDATA   maxval(met_lwup),minval(met_lwup): ',maxval(met_lwup),minval(met_lwup)
     768           0 :           write(iulog,*)'METDATA   maxval(met_qref),minval(met_qref): ',maxval(met_qref),minval(met_qref)
     769           0 :           write(iulog,*)'METDATA   maxval(met_tref),minval(met_tref): ',maxval(met_tref),minval(met_tref)
     770           0 :           write(iulog,*)'METDATA   maxval(met_u10),minval(met_u10): ',maxval(met_u10),minval(met_u10)
     771           0 :           write(iulog,*)'METDATA   maxval(met_sst),minval(met_sst): ',maxval(met_sst),minval(met_sst)
     772           0 :           write(iulog,*)'METDATA   maxval(met_icefrac),minval(met_icefrac): ',maxval(met_icefrac),minval(met_icefrac)
     773             :        endif
     774             :     endif
     775             : 
     776           0 :     do c = begchunk, endchunk
     777           0 :        call outfld('MET_TAUX',cam_in(c)%wsx , pcols   ,c   )
     778           0 :        call outfld('MET_TAUY',cam_in(c)%wsy , pcols   ,c   )
     779           0 :        call outfld('MET_LHFX',cam_in(c)%lhf , pcols   ,c   )
     780           0 :        call outfld('MET_SHFX',cam_in(c)%shf , pcols   ,c   )
     781           0 :        call outfld('MET_QFLX',cam_in(c)%cflx(:,1) , pcols   ,c   )
     782           0 :        call outfld('MET_ASDIR',cam_in(c)%asdir , pcols   ,c   )
     783           0 :        call outfld('MET_ASDIF',cam_in(c)%asdif , pcols   ,c   )
     784           0 :        call outfld('MET_ALDIR',cam_in(c)%aldir , pcols   ,c   )
     785           0 :        call outfld('MET_ALDIF',cam_in(c)%aldif , pcols   ,c   )
     786           0 :        call outfld('MET_LWUP',cam_in(c)%lwup , pcols   ,c   )
     787           0 :        call outfld('MET_QREF',cam_in(c)%qref , pcols   ,c   )
     788           0 :        call outfld('MET_TREF',cam_in(c)%tref , pcols   ,c   )
     789           0 :        call outfld('MET_U10',cam_in(c)%u10 , pcols   ,c   )
     790           0 :        call outfld('MET_SST',cam_in(c)%sst , pcols   ,c   )
     791           0 :        call outfld('MET_ICEFRAC',cam_in(c)%icefrac , pcols   ,c   )
     792             :     enddo
     793             : 
     794           0 :   end subroutine get_met_srf2
     795             : 
     796             : !-------------------------------------------------------------------
     797             : !-------------------------------------------------------------------
     798           0 :   subroutine get_met_srf1( cam_in )
     799           0 :     use camsrfexch,          only: cam_in_t
     800             :     use phys_grid,           only: get_ncols_p
     801             :     use cam_history,         only: outfld
     802             :     use shr_const_mod,       only: shr_const_stebol
     803             : 
     804             :     implicit none
     805             : 
     806             :     type(cam_in_t), intent(inout), dimension(begchunk:endchunk) :: cam_in
     807             : 
     808             :     integer :: c,ncol,i
     809             : 
     810           0 :     if (met_srf_feedback) return
     811           0 :     if (.not.has_ts) then
     812           0 :        call endrun('The meteorolgy input must have TS to run with met_srf_feedback set to FALSE')
     813             :     endif
     814             : 
     815           0 :     do c=begchunk,endchunk
     816           0 :        ncol = get_ncols_p(c)
     817           0 :        cam_in(c)%ts(:ncol)     = met_ts(:ncol,c)
     818           0 :        do i = 1,ncol
     819           0 :           cam_in(c)%snowhland(i) = met_snowh(i,c)*cam_in(c)%landfrac(i) * met_snowh_factor
     820             :        enddo
     821             :     end do ! Chunk loop
     822             : 
     823           0 :     if (debug) then
     824           0 :        if (masterproc) then
     825           0 :           write(iulog,*)'METDATA       maxval(met_ts),minval(met_ts): ',maxval(met_ts),minval(met_ts)
     826           0 :           write(iulog,*)'METDATA maxval(met_snowh),minval(met_snowh): ',maxval(met_snowh),minval(met_snowh)
     827             :        endif
     828             :     endif
     829             : 
     830           0 :     do c = begchunk, endchunk
     831           0 :        call outfld('MET_SNOWH',cam_in(c)%snowhland, pcols   ,c   )
     832             :     enddo
     833             : 
     834           0 :   end subroutine get_met_srf1
     835             : 
     836             : !-------------------------------------------------------------------
     837             : !-------------------------------------------------------------------
     838           0 :   subroutine get_ocn_ice_frcs( lndfrc, ocnfrc, icefrc, lchnk, ncol )
     839             : 
     840           0 :     use shr_const_mod, only: SHR_CONST_TKFRZSW
     841             :     use shr_const_mod, only: SHR_CONST_TKFRZ
     842             :     use cam_history,   only: outfld
     843             : 
     844             :     ! args
     845             :     real(r8), intent( in) :: lndfrc (pcols)
     846             :     real(r8), intent(out) :: ocnfrc (pcols)
     847             :     real(r8), intent(out) :: icefrc (pcols)
     848             : 
     849             :     integer, intent(in)   :: lchnk
     850             :     integer, intent(in)   :: ncol
     851             : 
     852             :     ! local vars
     853             :     integer :: i
     854             : 
     855           0 :     if (met_srf_sst) then
     856           0 :       do i = 1,ncol
     857             : 
     858             :         ! If configured for using SST, and ICEFRAC, then get icefrc
     859             :         ! directly from the meteorological data.
     860           0 :         if (met_icefrac(i,lchnk) == srf_fill_value) then
     861           0 :           icefrc(i) = 0._r8
     862             :         else
     863           0 :           icefrc(i) = min(met_icefrac(i,lchnk), 1._r8 - lndfrc(i))
     864             :         end if
     865           0 :         ocnfrc(i) = 1._r8 - lndfrc(i) - icefrc(i)
     866             :       enddo
     867             :     else
     868             : 
     869           0 :       if (.not.has_ts) then
     870           0 :          if (masterproc) then
     871           0 :             write(iulog,*) 'get_ocn_ice_frcs: TS is not in the met dataset and cannot set ocnfrc and icefrc'
     872           0 :             call endrun('get_ocn_ice_frcs: TS is not in the met dataset')
     873             :          endif
     874             :       endif
     875             : 
     876           0 :       do i = 1,ncol
     877             : 
     878           0 :          if ( met_ts(i,lchnk) < SHR_CONST_TKFRZ-2._r8 ) then
     879           0 :             ocnfrc(i) = 0._r8
     880           0 :             icefrc(i) = 1._r8 - lndfrc(i)
     881             :          else
     882           0 :             icefrc(i) = 0._r8
     883           0 :             ocnfrc(i) = 1._r8 - lndfrc(i)
     884             :          endif
     885             : 
     886             :       enddo
     887             :     end if
     888             : 
     889           0 :     call outfld('MET_TS', met_ts(:ncol,lchnk) , ncol   ,lchnk   )
     890           0 :     call outfld('MET_OCNFRC', ocnfrc(:ncol) , ncol   ,lchnk   )
     891           0 :     call outfld('MET_ICEFRC', icefrc(:ncol) , ncol   ,lchnk   )
     892             : 
     893           0 :   endsubroutine get_ocn_ice_frcs
     894             : 
     895             : !-------------------------------------------------------------------
     896             : ! allows access to physics state fields
     897             : !   q  : water vapor
     898             : !   ps : surface pressure
     899             : !   t  : temperature
     900             : !-------------------------------------------------------------------
     901           0 :   subroutine get_dyn_flds( state, tend, dt )
     902             : 
     903           0 :     use physics_types,  only: physics_state, physics_tend, physics_dme_adjust
     904             :     use ppgrid,         only: pcols, pver, begchunk, endchunk
     905             :     use phys_grid,      only: get_ncols_p
     906             :     use cam_history,    only: outfld
     907             :     use air_composition,only: thermodynamic_active_species_liq_num, thermodynamic_active_species_ice_num
     908             :     use air_composition,only: thermodynamic_active_species_liq_idx,thermodynamic_active_species_ice_idx
     909             : 
     910             :     implicit none
     911             : 
     912             :     type(physics_state), intent(inout), dimension(begchunk:endchunk) :: state
     913             :     type(physics_tend ), intent(inout), dimension(begchunk:endchunk) :: tend
     914             :     real(r8),            intent(in   ) :: dt                  ! model physics timestep
     915             : 
     916             :     integer :: lats(pcols)           ! array of latitude indices
     917             :     integer :: lons(pcols)           ! array of longitude indices
     918             :     integer :: c, ncol, i,j,k
     919             :     integer :: m_cnst,m
     920             :     real(r8):: qini(pcols,pver)      ! initial specific humidity
     921             :     real(r8):: totliqini(pcols,pver) ! initial total liquid
     922             :     real(r8):: toticeini(pcols,pver) ! initial total ice
     923             : 
     924             :     real(r8) :: tmp(pcols,pver)
     925             : 
     926           0 :     call t_startf('MET__GET_DYN2')
     927             : 
     928           0 :     do c = begchunk, endchunk
     929           0 :        ncol = get_ncols_p(c)
     930             :        !
     931             :        ! update water variables
     932             :        !
     933           0 :        qini(:ncol,:pver) = state(c)%q(:ncol,:pver,1)
     934           0 :        totliqini = 0.0_r8
     935           0 :        do m_cnst=1,thermodynamic_active_species_liq_num
     936           0 :          m = thermodynamic_active_species_liq_idx(m_cnst)
     937           0 :          totliqini(:ncol,:pver) = totliqini(:ncol,:pver)+state(c)%q(:ncol,:pver,m)
     938             :        end do
     939           0 :        toticeini = 0.0_r8
     940           0 :        do m_cnst=1,thermodynamic_active_species_ice_num
     941           0 :          m = thermodynamic_active_species_ice_idx(m_cnst)
     942           0 :          toticeini(:ncol,:pver) = toticeini(:ncol,:pver)+state(c)%q(:ncol,:pver,m)
     943             :        end do
     944             : 
     945           0 :        do k=1,pver
     946           0 :           do i=1,ncol
     947           0 :              if (met_nudge_temp) then
     948           0 :                 state(c)%t(i,k) = (1._r8-met_rlx(k))*state(c)%t(i,k) + met_rlx(k)*met_t(i,k,c)
     949             :              end if
     950             :              ! at this point tracer mixing ratios have already been
     951             :              ! converted from dry to moist
     952           0 :              state(c)%q(i,k,1) = alpha*state(c)%q(i,k,1) + (D1_0-alpha)*met_q(i,k,c)
     953             : 
     954           0 :              if ((state(c)%q(i,k,1) < D0_0).and. (alpha .ne. D1_0 )) state(c)%q(i,k,1) = D0_0
     955             : 
     956             :           end do
     957             : 
     958             :        end do
     959             : 
     960             :        ! now adjust mass of each layer now that water vapor has changed
     961           0 :        if (( .not. online_test ) .and. (alpha .ne. D1_0 )) then
     962             :           call physics_dme_adjust(state(c), tend(c), qini, totliqini, toticeini, dt)
     963             :        endif
     964             : 
     965             :     end do
     966             : 
     967           0 :     if (debug) then
     968           0 :     if (masterproc) then
     969           0 :       write(iulog,*)'METDATA maxval(met_t),minval(met_t): ', maxval(met_t),minval(met_t)
     970           0 :       write(iulog,*)'METDATA maxval(met_ps_next),minval(met_ps_next): ',  maxval(met_ps_next),minval(met_ps_next)
     971             :     endif
     972             :     endif
     973             : 
     974           0 :     do c = begchunk, endchunk
     975           0 :        call outfld('MET_T  ',state(c)%t , pcols   ,c   )
     976             :     enddo
     977           0 :     call t_stopf('MET__GET_DYN2')
     978             : 
     979           0 :   end subroutine get_dyn_flds
     980             : 
     981             : !------------------------------------------------------------------------
     982             : ! get the meteorological winds on the grid cell centers (A-grid)
     983             : !------------------------------------------------------------------------
     984           0 :   subroutine get_uv_centered( grid, u, v )
     985             : 
     986           0 :     use cam_history,    only: outfld
     987             : 
     988             :     implicit none
     989             : 
     990             :     type (T_FVDYCORE_GRID), intent(in) :: grid
     991             :     real(r8), intent(out) :: u(grid%im, grid%jfirst-grid%ng_d:grid%jlast+grid%ng_d, &
     992             :                                grid%kfirst:grid%klast)  ! u-wind on A-grid
     993             :     real(r8), intent(out) :: v(grid%im, grid%jfirst-grid%ng_s:grid%jlast+grid%ng_d, &
     994             :                                grid%kfirst:grid%klast)  ! v-wind on A-grid
     995             : 
     996             :     integer :: i,j,k
     997             : 
     998             :     integer :: jm, jfirst, jlast, jfirstxy, jlastxy, kfirst, klast, ng_d, ng_s, ifirstxy, ilastxy
     999             : 
    1000           0 :     real(r8) :: u3s_tmp(grid%ifirstxy:grid%ilastxy,grid%km)
    1001           0 :     real(r8) :: v3s_tmp(grid%ifirstxy:grid%ilastxy,grid%km)
    1002             : 
    1003           0 :     jm      = grid%jm
    1004           0 :     jfirstxy= grid%jfirstxy
    1005           0 :     jlastxy = grid%jlastxy
    1006           0 :     jfirst  = grid%jfirst
    1007           0 :     jlast   = grid%jlast
    1008           0 :     kfirst  = grid%kfirst
    1009           0 :     klast   = grid%klast
    1010           0 :     ifirstxy= grid%ifirstxy
    1011           0 :     ilastxy = grid%ilastxy
    1012             : 
    1013           0 :     ng_d    = grid%ng_d
    1014           0 :     ng_s    = grid%ng_s
    1015             : 
    1016           0 :     u(:,:,:) = D0_0
    1017           0 :     v(:,:,:) = D0_0
    1018             : 
    1019           0 :     u(         :,  max(1,jfirst-ng_d):min(jm,jlast+ng_d), kfirst:klast ) = &
    1020           0 :          met_u(:,  max(1,jfirst-ng_d):min(jm,jlast+ng_d), kfirst:klast )
    1021             : 
    1022           0 :     v(         :,  max(1,jfirst-ng_s):min(jm,jlast+ng_d), kfirst:klast ) = &
    1023           0 :          met_v(:,  max(1,jfirst-ng_s):min(jm,jlast+ng_d), kfirst:klast )
    1024             : 
    1025           0 :     if (masterproc) then
    1026           0 :     if (debug) write(iulog,*)'METDATA maxval(u),minval(u),maxval(v),minval(v) : ',&
    1027           0 :          maxval(u(:,  max(1,jfirst-ng_d):min(jm,jlast+ng_d), kfirst:klast )),&
    1028           0 :          minval(u(:,  max(1,jfirst-ng_d):min(jm,jlast+ng_d), kfirst:klast )),&
    1029           0 :          maxval(v(:,  max(1,jfirst-ng_s):min(jm,jlast+ng_d), kfirst:klast )),&
    1030           0 :          minval(v(:,  max(1,jfirst-ng_s):min(jm,jlast+ng_d), kfirst:klast ))
    1031             :     endif
    1032             : 
    1033           0 :     if ( grid%twod_decomp == 0 ) then
    1034           0 :        do j = jfirst, jlast
    1035           0 :           do k = kfirst, klast
    1036           0 :              do i = 1, grid%im
    1037           0 :                 u3s_tmp(i,k) = u(i,j,k)
    1038           0 :                 v3s_tmp(i,k) = v(i,j,k)
    1039             :              enddo
    1040             :           enddo
    1041           0 :           call outfld ('MET_U ', u3s_tmp, grid%im, j )
    1042           0 :           call outfld ('MET_V ', v3s_tmp, grid%im, j )
    1043             :        enddo
    1044             :     endif
    1045             : 
    1046           0 :   end subroutine get_uv_centered
    1047             : 
    1048             : !------------------------------------------------------------------------
    1049             : ! get the meteorological surface pressure interp to dyn substep
    1050             : !------------------------------------------------------------------------
    1051           0 :   subroutine get_ps( grid, ps, nsubsteps, n )
    1052             : 
    1053           0 :     use cam_history,    only: outfld
    1054             : 
    1055             :     implicit none
    1056             : 
    1057             :     type (T_FVDYCORE_GRID), intent(in) :: grid
    1058             :     real(r8), intent(out) :: ps(grid%im, grid%jfirst:grid%jlast)
    1059             :     integer, intent(in) :: nsubsteps
    1060             :     integer, intent(in) :: n
    1061             : 
    1062             :     real(r8) :: num1, num2
    1063             :     integer  :: j
    1064             : 
    1065           0 :     num1 = n
    1066           0 :     num2 = nsubsteps
    1067             : 
    1068           0 :     ps(:,:) = met_ps_curr(:,:) + num1*(met_ps_next(:,:)-met_ps_curr(:,:))/num2
    1069             : 
    1070           0 :     if ( grid%twod_decomp == 0 ) then
    1071           0 :        do j = grid%jfirst, grid%jlast
    1072           0 :           call outfld('MET_PS',ps(:,j), grid%im   ,j   )
    1073             :        enddo
    1074             :     endif
    1075           0 :   end subroutine get_ps
    1076             : 
    1077             : !------------------------------------------------------------------------
    1078             : ! get the meteorological winds on the grid cell walls (vorticity winds)
    1079             : !   us : staggered zonal wind
    1080             : !   vs : staggered meridional wind
    1081             : !------------------------------------------------------------------------
    1082           0 :   subroutine get_us_vs( grid, us, vs )
    1083             : 
    1084             :     implicit none
    1085             : 
    1086             :     type (T_FVDYCORE_GRID), intent(in) :: grid
    1087             :     real(r8), intent(inout) :: us(grid%im, grid%jfirst-grid%ng_d:grid%jlast+grid%ng_s, &
    1088             :                                   grid%kfirst:grid%klast)    ! u-wind on d-grid
    1089             :     real(r8), intent(inout) :: vs(grid%im, grid%jfirst-grid%ng_s:grid%jlast+grid%ng_d, &
    1090             :                                   grid%kfirst:grid%klast)    ! v-wind on d-grid
    1091             : 
    1092             :     integer :: i,j,k
    1093             : 
    1094             :     integer :: jm, jfirst, jlast, kfirst, klast, ng_d, ng_s
    1095             : 
    1096           0 :     jm      = grid%jm
    1097           0 :     jfirst  = grid%jfirst
    1098           0 :     jlast   = grid%jlast
    1099           0 :     kfirst  = grid%kfirst
    1100           0 :     klast   = grid%klast
    1101           0 :     ng_d    = grid%ng_d
    1102           0 :     ng_s    = grid%ng_s
    1103             : 
    1104           0 :     call t_startf('MET__get_us_vs')
    1105             : 
    1106             :     ! vertical relaxation (blending) occurs in dyn_run (dyn_comp.F90)
    1107             : 
    1108           0 :     us(:,:,:) = 1.e36_r8
    1109           0 :     vs(:,:,:) = 1.e36_r8
    1110           0 :     us(          :, max(2,jfirst):  min(jm,jlast), kfirst:klast) =      &
    1111           0 :          met_us( :, max(2,jfirst):  min(jm,jlast), kfirst:klast)
    1112           0 :     vs(          :, max(1,jfirst):  min(jm,jlast), kfirst:klast) =      &
    1113           0 :          met_vs( :, max(1,jfirst):  min(jm,jlast), kfirst:klast)
    1114           0 :     if (masterproc) then
    1115           0 :     if (debug) write(iulog,*)grid%iam,': METDATA maxval(us),minval(us),maxval(vs),minval(vs) : ',&
    1116           0 :          maxval(us(          :, max(2,jfirst):  min(jm,jlast),    kfirst:klast)),&
    1117           0 :          minval(us(          :, max(2,jfirst):  min(jm,jlast),    kfirst:klast)),&
    1118           0 :          maxval(vs(          :, max(1,jfirst):  min(jm,jlast),    kfirst:klast)),&
    1119           0 :          minval(vs(          :, max(1,jfirst):  min(jm,jlast),    kfirst:klast))
    1120             :     endif
    1121             : 
    1122             : !!$    if (debug) then
    1123             : !!$        u3s_tmp = 1.e36
    1124             : !!$       do j = jfirst, jlast
    1125             : !!$          do k = kfirst, klast
    1126             : !!$             do i = 1, im
    1127             : !!$                if (j >= 2) u3s_tmp(i,k) = us(i,j,k)
    1128             : !!$                v3s_tmp(i,k) = vs(i,j,k)
    1129             : !!$             enddo
    1130             : !!$          enddo
    1131             : !!$          call outfld ('MET_US ', u3s_tmp, im, j )
    1132             : !!$          call outfld ('MET_VS ', v3s_tmp, im, j )
    1133             : !!$       enddo
    1134             : !!$    endif
    1135             : !!$
    1136           0 :     call t_stopf('MET__get_us_vs')
    1137             : 
    1138           0 :   end subroutine get_us_vs
    1139             : 
    1140             : !-------------------------------------------------------------------------
    1141             : ! writes file names to restart file
    1142             : !-------------------------------------------------------------------------
    1143             : 
    1144           0 :   subroutine write_met_restart_pio(File)
    1145             :     type(file_desc_t), intent(inout) :: File
    1146             :     integer :: ierr
    1147           0 :     ierr =  pio_put_att(File, PIO_GLOBAL, 'current_metdata_filename', curr_filename)
    1148           0 :     ierr =  pio_put_att(File, PIO_GLOBAL, 'next_metdata_filename', next_filename)
    1149             : 
    1150           0 :   end subroutine write_met_restart_pio
    1151           0 :   subroutine read_met_restart_pio(File)
    1152             :     type(file_desc_t), intent(inout) :: File
    1153             : 
    1154             :     integer :: ierr, xtype
    1155             :     integer(pio_offset_kind) :: slen
    1156             : 
    1157           0 :     ierr = pio_inq_att(File, PIO_GLOBAL, 'current_metdata_filename',xtype, slen)
    1158           0 :     ierr = pio_get_att(File, PIO_GLOBAL, 'current_metdata_filename', curr_filename)
    1159           0 :     curr_filename(slen+1:256) = ''
    1160             : 
    1161           0 :     ierr = pio_inq_att(File, PIO_GLOBAL, 'next_metdata_filename',xtype, slen)
    1162           0 :     ierr = pio_get_att(File, PIO_GLOBAL, 'next_metdata_filename', next_filename)
    1163           0 :     next_filename(slen+1:256) = ''
    1164             : 
    1165           0 :   end subroutine read_met_restart_pio
    1166             : 
    1167           0 :   subroutine write_met_restart_bin( nrg )
    1168             :     implicit none
    1169             :     integer,intent(in) :: nrg     ! Unit number
    1170             :     integer :: ioerr   ! error status
    1171             : 
    1172           0 :     if (masterproc) then
    1173           0 :        write(nrg, iostat=ioerr) curr_filename
    1174           0 :        if (ioerr /= 0 ) then
    1175           0 :           write(iulog,*) 'WRITE ioerror ',ioerr,' on i/o unit = ',nrg
    1176           0 :           call endrun ('WRITE_RESTART_DYNAMICS')
    1177             :        end if
    1178           0 :        write(nrg, iostat=ioerr) next_filename
    1179           0 :        if (ioerr /= 0 ) then
    1180           0 :           write(iulog,*) 'WRITE ioerror ',ioerr,' on i/o unit = ',nrg
    1181           0 :           call endrun ('WRITE_RESTART_DYNAMICS')
    1182             :        end if
    1183             :     end if
    1184           0 :   end subroutine write_met_restart_bin
    1185             : 
    1186             : !-------------------------------------------------------------------------
    1187             : ! reads file names from restart file
    1188             : !-------------------------------------------------------------------------
    1189           0 :   subroutine read_met_restart_bin( nrg )
    1190             :     implicit none
    1191             :     integer,intent(in) :: nrg     ! Unit number
    1192             :     integer :: ioerr   ! error status
    1193             : 
    1194           0 :     if (masterproc) then
    1195           0 :        read(nrg, iostat=ioerr) curr_filename
    1196           0 :        if (ioerr /= 0 ) then
    1197           0 :           write(iulog,*) 'READ ioerror ',ioerr,' on i/o unit = ',nrg
    1198           0 :           call endrun ('READ_RESTART_DYNAMICS')
    1199             :        end if
    1200           0 :        read(nrg, iostat=ioerr) next_filename
    1201           0 :        if (ioerr /= 0 ) then
    1202           0 :           write(iulog,*) 'READ ioerror ',ioerr,' on i/o unit = ',nrg
    1203           0 :           call endrun ('READ_RESTART_DYNAMICS')
    1204             :        end if
    1205             :     end if
    1206             : 
    1207             : #if ( defined SPMD )
    1208           0 :     call mpibcast ( curr_filename ,len(curr_filename) ,mpichar,0,mpicom)
    1209           0 :     call mpibcast ( next_filename ,len(next_filename) ,mpichar,0,mpicom)
    1210             : #endif
    1211           0 :   end subroutine read_met_restart_bin
    1212             : 
    1213             : !-------------------------------------------------------------------------
    1214             : ! returns true if the met winds are defined on cell walls
    1215             : !-------------------------------------------------------------------------
    1216           0 :   function met_winds_on_walls()
    1217             :     logical :: met_winds_on_walls
    1218             : 
    1219           0 :     met_winds_on_walls = met_cell_wall_winds
    1220           0 :   end function met_winds_on_walls
    1221             : 
    1222             : ! internal methods :
    1223             : 
    1224             : !-------------------------------------------------------------------------
    1225             : ! transfers cell-centered winds to cell walls
    1226             : !-------------------------------------------------------------------------
    1227           0 :   subroutine transfer_windsToWalls(grid)
    1228             : 
    1229             :     implicit none
    1230             : 
    1231             :     type (T_FVDYCORE_GRID), intent(in) :: grid
    1232             :     integer :: i,j,k
    1233             :     integer :: im, jfirst, jlast, kfirst, klast
    1234             : 
    1235           0 :     im      = grid%im
    1236           0 :     jfirst  = grid%jfirst
    1237           0 :     jlast   = grid%jlast
    1238           0 :     kfirst  = grid%kfirst
    1239           0 :     klast   = grid%klast
    1240             : 
    1241           0 :     call t_startf('MET__transfer_windsToWalls')
    1242             : 
    1243             : !$omp parallel do private (i, j, k)
    1244           0 :     do k = kfirst, klast
    1245             : 
    1246           0 :        do j = jfirst+1,jlast
    1247           0 :           do i = 1,im
    1248           0 :              met_us(i,j,k) = ( met_u(i,j,k) + met_u(i,j-1,k) )*D0_5
    1249             :           end do
    1250             :        end do
    1251             : 
    1252             : #if defined( SPMD )
    1253           0 :        if ( jfirst .gt. 1 ) then
    1254           0 :           do i = 1, im
    1255             :              ! met_u is alread ghosted at this point
    1256           0 :              met_us(i,jfirst,k) = ( met_u(i,jfirst,k) + met_u(i,jfirst-1,k) )*D0_5
    1257             :           enddo
    1258             :        endif
    1259             : #endif
    1260             : 
    1261           0 :        do j = jfirst,jlast
    1262           0 :           met_vs(1,j,k) = ( met_v(1,j,k) + met_v(im,j,k) )*D0_5
    1263           0 :           do i = 2,im
    1264           0 :              met_vs(i,j,k) = ( met_v(i,j,k) + met_v(i-1,j,k) )*D0_5
    1265             :           end do
    1266             :        end do
    1267             :     end do
    1268             : 
    1269           0 :     call t_stopf('MET__transfer_windsToWalls')
    1270             : 
    1271           0 :   end subroutine transfer_windsToWalls
    1272             : 
    1273           0 :   subroutine get_model_time()
    1274             :     implicit none
    1275             :     integer yr, mon, day, ncsec  ! components of a date
    1276             : 
    1277           0 :     call t_startf('MET__get_model_time')
    1278             : 
    1279           0 :     call get_curr_date(yr, mon, day, ncsec)
    1280             : 
    1281           0 :     curr_mod_time = get_time_float( yr, mon, day, ncsec )
    1282           0 :     next_mod_time = curr_mod_time + get_step_size()/seconds_per_day
    1283             : 
    1284           0 :     call t_stopf('MET__get_model_time')
    1285             : 
    1286           0 :   end subroutine get_model_time
    1287             : 
    1288             : !------------------------------------------------------------------------------
    1289             : !------------------------------------------------------------------------------
    1290           0 :   subroutine check_files()
    1291             : 
    1292             :     use shr_sys_mod, only: shr_sys_system
    1293             :     use ioFileMod, only: getfil
    1294             : 
    1295             :     implicit none
    1296             : 
    1297             : !-----------------------------------------------------------------------
    1298             : !       ... local variables
    1299             : !-----------------------------------------------------------------------
    1300             :     character(len=256) :: ctmp
    1301             :     character(len=256) :: loc_fname
    1302             :     integer            :: istat
    1303             : 
    1304             : 
    1305           0 :     if (next_mod_time > curr_data_times(size(curr_data_times))) then
    1306           0 :        if ( .not. associated(next_data_times) ) then
    1307             :           ! open next file...
    1308           0 :           next_filename = incr_filename( curr_filename )
    1309           0 :           call open_met_datafile( next_filename, next_fileid, next_data_times, met_data_path )
    1310             :        endif
    1311             :     endif
    1312             : 
    1313           0 :     if ( associated(next_data_times) ) then
    1314           0 :        if (curr_mod_time >= next_data_times(1)) then
    1315             : 
    1316             :           ! close current file ...
    1317           0 :           call pio_closefile( curr_fileid )
    1318           0 :           if (masterproc) then
    1319             :              ! remove if requested
    1320           0 :              if( met_remove_file ) then
    1321           0 :                 call getfil( curr_filename, loc_fname, 0 )
    1322           0 :                 write(iulog,*) 'check_files: removing file = ',trim(loc_fname)
    1323           0 :                 ctmp = 'rm -f ' // trim(loc_fname)
    1324           0 :                 write(iulog,*) 'check_files: fsystem issuing command - '
    1325           0 :                 write(iulog,*) trim(ctmp)
    1326           0 :                 call shr_sys_system( ctmp, istat )
    1327             :              end if
    1328             :           endif
    1329             : 
    1330           0 :           curr_filename = next_filename
    1331           0 :           curr_fileid = next_fileid
    1332             : 
    1333           0 :           deallocate( curr_data_times )
    1334           0 :           allocate( curr_data_times( size( next_data_times ) ) )
    1335           0 :           curr_data_times(:) = next_data_times(:)
    1336             : 
    1337           0 :           next_filename = ''
    1338             : 
    1339           0 :           deallocate( next_data_times )
    1340             :           nullify(  next_data_times )
    1341             : 
    1342             :        endif
    1343             :     endif
    1344             : 
    1345           0 :   end subroutine check_files
    1346             : 
    1347             : !-----------------------------------------------------------------------
    1348             : !-----------------------------------------------------------------------
    1349           0 :   function incr_filename( filename )
    1350             : 
    1351             :     !-----------------------------------------------------------------------
    1352             :     !   ... Increment or decrement a date string withing a filename
    1353             :     !           the filename date section is assumed to be of the form
    1354             :     !           yyyy-dd-mm
    1355             :     !-----------------------------------------------------------------------
    1356             : 
    1357             :     use string_utils, only : incstr
    1358             :     use shr_file_mod, only : shr_file_getunit, shr_file_freeunit
    1359             : 
    1360             :     implicit none
    1361             : 
    1362             : 
    1363             :     character(len=*), intent(in) :: filename ! present dynamical dataset filename
    1364             :     character(len=256) :: incr_filename      ! next filename in the sequence
    1365             :     character(len=*), parameter :: subname = 'incr_filename'
    1366             : 
    1367             :     ! set new next_filename ...
    1368             : 
    1369             :     !-----------------------------------------------------------------------
    1370             :     !   ... local variables
    1371             :     !-----------------------------------------------------------------------
    1372             :     integer :: pos, pos1, istat
    1373             :     character(len=256) :: fn_new, line
    1374             :     character(len=6)   :: seconds
    1375             :     character(len=5)   :: num
    1376             :     integer :: ios,unitnumber
    1377             : 
    1378           0 :     if ( len_trim(met_filenames_list) == 0) then
    1379             :        !-----------------------------------------------------------------------
    1380             :        !        ... ccm type filename
    1381             :        !-----------------------------------------------------------------------
    1382           0 :        pos = len_trim( filename )
    1383           0 :        fn_new = filename(:pos)
    1384           0 :        if (masterproc) write(iulog,*) 'incr_flnm: old filename = ',trim(fn_new)
    1385           0 :        if( fn_new(pos-2:) == '.nc' ) then
    1386           0 :           pos = pos - 3
    1387             :        end if
    1388           0 :        istat = incstr( fn_new(:pos), 1 )
    1389           0 :        if( istat /= 0 ) then
    1390           0 :           write(iulog,*) 'incr_flnm: incstr returned ', istat
    1391           0 :           write(iulog,*) '           while trying to decrement ',trim( fn_new )
    1392           0 :           call endrun (subname // ':: ERRROR - check atm.log for error message')
    1393             :        end if
    1394             : 
    1395             :     else
    1396             : 
    1397             :        ! open met_filenames_list
    1398           0 :        if (masterproc) write(iulog,*) 'incr_flnm: old filename = ',trim(filename)
    1399           0 :        if (masterproc) write(iulog,*) 'incr_flnm: open met_filenames_list : ',met_filenames_list
    1400           0 :        unitnumber = shr_file_getUnit()
    1401           0 :        open( unit=unitnumber, file=met_filenames_list, iostat=ios, status="OLD")
    1402           0 :        if (ios /= 0) then
    1403           0 :           call endrun('not able to open met_filenames_list file: '//met_filenames_list)
    1404             :        endif
    1405             : 
    1406             :        ! read file names
    1407           0 :        read( unit=unitnumber, fmt='(A)', iostat=ios ) line
    1408           0 :        if (ios /= 0) then
    1409           0 :           call endrun('not able to increment file name from met_filenames_list file: '//met_filenames_list)
    1410             :        endif
    1411           0 :        do while( trim(line) /= trim(filename) )
    1412           0 :           read( unit=unitnumber, fmt='(A)', iostat=ios ) line
    1413           0 :           if (ios /= 0) then
    1414           0 :              call endrun('not able to increment file name from met_filenames_list file: '//met_filenames_list)
    1415             :           endif
    1416             :        enddo
    1417             : 
    1418           0 :        read( unit=unitnumber, fmt='(A)', iostat=ios ) line
    1419           0 :        if (ios /= 0) then
    1420           0 :           call endrun('not able to increment file name from met_filenames_list file: '//met_filenames_list)
    1421             :        endif
    1422           0 :        fn_new = trim(line)
    1423             : 
    1424           0 :        close(unit=unitnumber)
    1425           0 :        call shr_file_freeUnit(unitnumber)
    1426             :     endif
    1427           0 :     incr_filename = trim(fn_new)
    1428           0 :     if (masterproc) write(iulog,*) 'incr_flnm: new filename = ',incr_filename
    1429             : 
    1430           0 :   end function incr_filename
    1431             : 
    1432             : !------------------------------------------------------------------------------
    1433             : !------------------------------------------------------------------------------
    1434           0 :   subroutine find_times( itms, fids, datatm, datatp, time )
    1435             : 
    1436             :     implicit none
    1437             : 
    1438             :     integer, intent(out) :: itms(2) ! record numbers that bracket time
    1439             :     type(file_desc_t), intent(out) :: fids(2) ! ids of files that contains these recs
    1440             :     real(r8), intent(in) :: time    ! time of interest
    1441             :     real(r8), intent(out):: datatm, datatp
    1442             :     character(len=*), parameter :: subname = 'find_times'
    1443             : 
    1444             :     integer np1        ! current forward time index of dataset
    1445             :     integer n,i      !
    1446             :     integer :: curr_tsize, next_tsize, all_tsize
    1447             : 
    1448           0 :     real(r8), allocatable, dimension(:):: all_data_times
    1449             : 
    1450           0 :     curr_tsize = size(curr_data_times)
    1451           0 :     next_tsize = 0
    1452           0 :     if ( associated(next_data_times)) next_tsize = size(next_data_times)
    1453             : 
    1454           0 :     all_tsize = curr_tsize + next_tsize
    1455             : 
    1456           0 :     allocate( all_data_times( all_tsize ) )
    1457             : 
    1458           0 :     all_data_times(:curr_tsize) = curr_data_times(:)
    1459           0 :     if (next_tsize > 0) all_data_times(curr_tsize+1:all_tsize) = next_data_times(:)
    1460             : 
    1461             :     ! find bracketing times
    1462           0 :     do n=1, all_tsize-1
    1463           0 :        np1 = n + 1
    1464           0 :        datatm = all_data_times(n)
    1465           0 :        datatp = all_data_times(np1)
    1466           0 :        if ( (time .ge. datatm) .and. (time .le. datatp) ) then
    1467             :           goto 20
    1468             :        endif
    1469             :     enddo
    1470             : 
    1471           0 :     write(iulog,*)'FIND_TIMES: Failed to find dates bracketing desired time =', time
    1472           0 :     write(iulog,*)' datatm = ',datatm
    1473           0 :     write(iulog,*)' datatp = ',datatp
    1474           0 :     write(iulog,*)' all_data_times = ',all_data_times
    1475             : 
    1476           0 :     call endrun (subname // ':: ERRROR - check atm.log for error message')
    1477             : 
    1478             : 20  continue
    1479             : 
    1480           0 :     deallocate( all_data_times )
    1481             : 
    1482           0 :     itms(1) = n
    1483           0 :     itms(2) = np1
    1484           0 :     fids(:) = curr_fileid
    1485             : 
    1486           0 :     do i=1,2
    1487           0 :        if ( itms(i) > curr_tsize ) then
    1488           0 :           itms(i) = itms(i) - curr_tsize
    1489           0 :           fids(i) = next_fileid
    1490             :        endif
    1491             :     enddo
    1492             : 
    1493           0 :   end subroutine find_times
    1494             : 
    1495             : !------------------------------------------------------------------------------
    1496             : !------------------------------------------------------------------------------
    1497           0 :   subroutine read_next_ps(grid)
    1498             :     use ncdio_atm,          only: infld
    1499             : 
    1500             :     implicit none
    1501             : 
    1502             :     type (T_FVDYCORE_GRID), intent(in) :: grid
    1503             : 
    1504             :     integer :: recnos(2)
    1505           0 :     type(file_desc_t) :: fids(2)
    1506             :     character(len=8) :: varname
    1507             :     integer :: ifirstxy, ilastxy, jfirstxy, jlastxy
    1508             : 
    1509           0 :     real(r8) :: wrk_xy(grid%ifirstxy:grid%ilastxy, grid%jfirstxy:grid%jlastxy )
    1510             : 
    1511             :     logical :: readvar
    1512             : 
    1513           0 :     if(online_test) then
    1514           0 :        varname='arch_PS'
    1515             :     else
    1516           0 :        varname='PS'
    1517             :     end if
    1518             : 
    1519           0 :     jfirstxy= grid%jfirstxy
    1520           0 :     jlastxy = grid%jlastxy
    1521           0 :     ifirstxy= grid%ifirstxy
    1522           0 :     ilastxy = grid%ilastxy
    1523             : 
    1524           0 :     call find_times( recnos, fids, datatimemn, datatimepn, next_mod_time )
    1525             : 
    1526             :     call infld(varname, fids(1), 'lon', 'lat',  ifirstxy, ilastxy, jfirstxy, jlastxy, &
    1527           0 :          wrk_xy, readvar, gridname='fv_centers', timelevel=recnos(1))
    1528             : 
    1529             :     ! transpose xy -> yz decomposition
    1530           0 :     call transpose_xy2yz_2d( wrk_xy, met_psi_next(nm)%data, grid )
    1531             : 
    1532             :     call infld(varname, fids(2), 'lon', 'lat',  ifirstxy, ilastxy, jfirstxy, jlastxy, &
    1533           0 :          wrk_xy, readvar, gridname='fv_centers', timelevel=recnos(2))
    1534             : 
    1535             :     ! transpose xy -> yz decomposition
    1536           0 :     call transpose_xy2yz_2d( wrk_xy, met_psi_next(np)%data, grid )
    1537             : 
    1538           0 :     if(masterproc) write(iulog,*)'READ_NEXT_PS: Read meteorological data  '
    1539             : 
    1540           0 :   end subroutine read_next_ps
    1541             : 
    1542             : !------------------------------------------------------------------------
    1543             : !------------------------------------------------------------------------
    1544           0 :   subroutine transpose_xy2yz_2d( xy_2d, yz_2d, grid )
    1545             : 
    1546             :     implicit none
    1547             :     type (T_FVDYCORE_GRID), intent(in)  :: grid
    1548             :     real(r8),               intent(in)  :: xy_2d(grid%ifirstxy:grid%ilastxy, grid%jfirstxy:grid%jlastxy)
    1549             :     real(r8),               intent(out) :: yz_2d(1:grid%im, grid%jfirst:grid%jlast)
    1550             : 
    1551           0 :     real(r8) :: xy3(grid%ifirstxy:grid%ilastxy, grid%jfirstxy:grid%jlastxy, 1:grid%npr_z )
    1552             :     integer :: i,j,k
    1553             : 
    1554           0 :     if (grid%iam .lt. grid%npes_xy) then
    1555           0 :        if ( grid%twod_decomp == 1 ) then
    1556             : 
    1557             : #if defined( SPMD )
    1558             : !$omp parallel do private(i,j,k)
    1559           0 :           do k=1,grid%npr_z
    1560           0 :              do j=grid%jfirstxy,grid%jlastxy
    1561           0 :                 do i=grid%ifirstxy,grid%ilastxy
    1562           0 :                    xy3(i,j,k) = xy_2d(i,j)
    1563             :                 enddo
    1564             :              enddo
    1565             :           enddo
    1566             : 
    1567             :           call mp_sendirr(grid%commxy, grid%xy2d_to_yz2d%SendDesc, &
    1568             :                grid%xy2d_to_yz2d%RecvDesc, xy3, yz_2d,             &
    1569           0 :                modc=grid%modc_dynrun )
    1570             :           call mp_recvirr(grid%commxy, grid%xy2d_to_yz2d%SendDesc, &
    1571             :                grid%xy2d_to_yz2d%RecvDesc, xy3, yz_2d,             &
    1572           0 :                modc=grid%modc_dynrun )
    1573             : #endif
    1574             : 
    1575             :        else
    1576           0 :           yz_2d(:,:) = xy_2d(:,:)
    1577             :        endif
    1578             :     endif  !  (grid%iam .lt. grid%npes_xy)
    1579             : 
    1580           0 :   end subroutine transpose_xy2yz_2d
    1581             : 
    1582             : !------------------------------------------------------------------------
    1583             : !------------------------------------------------------------------------
    1584           0 :   subroutine transpose_xy2yz_3d( xy_3d, yz_3d, grid )
    1585             : 
    1586             :     implicit none
    1587             :     type (T_FVDYCORE_GRID), intent(in)  :: grid
    1588             :     real(r8),               intent(in)  :: xy_3d(grid%ifirstxy:grid%ilastxy, grid%jfirstxy:grid%jlastxy, 1:grid%km )
    1589             :     real(r8),               intent(out) :: yz_3d(1:grid%im, grid%jfirst:grid%jlast, grid%kfirst:grid%klast)
    1590             : 
    1591           0 :     if (grid%iam .lt. grid%npes_xy) then
    1592           0 :        if ( grid%twod_decomp == 1 ) then
    1593             : #if defined( SPMD )
    1594             :           call mp_sendirr( grid%commxy, grid%ijk_xy_to_yz%SendDesc, &
    1595             :                grid%ijk_xy_to_yz%RecvDesc, xy_3d, yz_3d,            &
    1596           0 :                modc=grid%modc_dynrun )
    1597             :           call mp_recvirr( grid%commxy, grid%ijk_xy_to_yz%SendDesc, &
    1598             :                grid%ijk_xy_to_yz%RecvDesc, xy_3d, yz_3d,            &
    1599           0 :                modc=grid%modc_dynrun )
    1600             : #endif
    1601             :        else
    1602           0 :           yz_3d(:,:,:) = xy_3d(:,:,:)
    1603             :        endif
    1604             :     endif  !  (grid%iam .lt. grid%npes_xy)
    1605             : 
    1606           0 :   end subroutine transpose_xy2yz_3d
    1607             : 
    1608             : 
    1609             : 
    1610             : !------------------------------------------------------------------------
    1611             : !------------------------------------------------------------------------
    1612           0 :   subroutine read_next_metdata(grid)
    1613             :     use ncdio_atm,          only: infld
    1614             :     use cam_grid_support,   only: cam_grid_check, cam_grid_id
    1615             :     use cam_grid_support,   only: cam_grid_get_dim_names
    1616             : 
    1617             :     implicit none
    1618             : 
    1619             :     type (T_FVDYCORE_GRID), intent(in) :: grid
    1620             :     integer recnos(2),  i      !
    1621           0 :     type(file_desc_t) :: fids(2)
    1622             : 
    1623             :     character(len=8) :: Uname, Vname, Tname, Qname, psname
    1624             :     character(len=8) :: dim1name, dim2name
    1625             :     integer :: im, jm, km
    1626             :     logical :: readvar
    1627             :     integer :: ifirstxy, ilastxy, jfirstxy, jlastxy
    1628           0 :     real(r8), allocatable :: wrk2_xy(:,:)
    1629           0 :     real(r8), allocatable :: wrk3_xy(:,:,:)
    1630           0 :     real(r8), allocatable :: tmp_data(:,:,:)
    1631             :     integer :: elev1,blev1, elev2,blev2
    1632             :     integer :: elev3,blev3, elev4,blev4
    1633             :     integer :: grid_id  ! grid ID for data mapping
    1634             : 
    1635           0 :     call t_startf('MET__read_next_metdata')
    1636             : 
    1637           0 :     jfirstxy= grid%jfirstxy
    1638           0 :     jlastxy = grid%jlastxy
    1639           0 :     ifirstxy= grid%ifirstxy
    1640           0 :     ilastxy = grid%ilastxy
    1641             : 
    1642           0 :     im      = grid%im
    1643           0 :     jm      = grid%jm
    1644           0 :     km      = grid%km
    1645             : 
    1646           0 :     call find_times( recnos, fids, datatimem, datatimep, curr_mod_time )
    1647             :     !
    1648             :     ! Set up hyperslab corners
    1649             :     !
    1650             : 
    1651           0 :     if(online_test) then
    1652           0 :        Tname='arch_T'
    1653           0 :        Qname='arch_Q'
    1654           0 :        PSname='arch_PS'
    1655           0 :        if(met_cell_wall_winds) then
    1656           0 :           Uname='arch_US'
    1657           0 :           Vname='arch_VS'
    1658             :        else
    1659           0 :           Uname='arch_U'
    1660           0 :           Vname='arch_V'
    1661             :        end if
    1662             :     else
    1663           0 :        Tname='T'
    1664           0 :        Qname='Q'
    1665           0 :        PSname='PS'
    1666           0 :        if(met_cell_wall_winds) then
    1667           0 :           Uname='US'
    1668           0 :           Vname='VS'
    1669             :        else
    1670           0 :           Uname='U'
    1671           0 :           Vname='V'
    1672             :        end if
    1673             : 
    1674             :     end if
    1675             : 
    1676             : 
    1677           0 :     if ( num_met_levels>km ) then
    1678             : 
    1679           0 :        blev1 = 1
    1680           0 :        elev1 = km
    1681             : 
    1682           0 :        blev2 = num_met_levels-km+1
    1683           0 :        elev2 = num_met_levels
    1684             : 
    1685           0 :        blev3 = num_met_levels-km+1
    1686           0 :        elev3 = num_met_levels
    1687             : 
    1688           0 :        blev4 = 1
    1689           0 :        elev4 = num_met_levels
    1690             : 
    1691             :     else
    1692             : 
    1693           0 :        blev1 = km-num_met_levels+1
    1694           0 :        elev1 = km
    1695             : 
    1696           0 :        blev2 = 1
    1697           0 :        elev2 = num_met_levels
    1698             : 
    1699           0 :        blev3 = 1
    1700           0 :        elev3 = km
    1701             : 
    1702           0 :        blev4 = km-num_met_levels+1
    1703           0 :        elev4 = km
    1704             : 
    1705             :     endif
    1706             : 
    1707           0 :     allocate(tmp_data(pcols, 1:num_met_levels, begchunk:endchunk))
    1708           0 :     allocate(wrk2_xy(ifirstxy:ilastxy, jfirstxy:jlastxy))
    1709           0 :     allocate(wrk3_xy(ifirstxy:ilastxy, jfirstxy:jlastxy, 1:max(km,num_met_levels)))
    1710             : 
    1711             :     ! physgrid intput for FV is probably always lon/lat but let's be pedantic
    1712           0 :     grid_id = cam_grid_id('physgrid')
    1713           0 :     if (.not. cam_grid_check(grid_id)) then
    1714           0 :        call endrun('read_next_metdata: Internal error, no "physgrid" grid')
    1715             :     end if
    1716           0 :     call cam_grid_get_dim_names(grid_id, dim1name, dim2name)
    1717             : 
    1718           0 :     do i=1,2
    1719             : 
    1720           0 :        met_ti(i)%data = 0._r8
    1721             : 
    1722             :        call infld(Tname, fids(i), dim1name, 'lev', dim2name,  1, pcols, 1,num_met_levels , &
    1723           0 :             begchunk, endchunk, tmp_data, readvar, gridname='physgrid',timelevel=recnos(i))
    1724             : 
    1725           0 :        met_ti(i)%data(:,blev1:elev1,:) = tmp_data(:, blev2:elev2, :)
    1726             : 
    1727           0 :        met_qi(i)%data = 0._r8
    1728             : 
    1729             :        call infld(Qname, fids(i), dim1name, 'lev', dim2name,  1, pcols, 1,num_met_levels, &
    1730           0 :             begchunk, endchunk, tmp_data, readvar, gridname='physgrid',timelevel=recnos(i))
    1731             : 
    1732           0 :        met_qi(i)%data(:,blev1:elev1,:) = tmp_data(:, blev2:elev2, :)
    1733             : 
    1734           0 :        if (met_cell_wall_winds) then
    1735             : 
    1736           0 :           wrk3_xy = 0._r8
    1737           0 :           met_usi(i)%data(:,:,:) = 0._r8
    1738             :           call infld(Uname, fids(i), 'lon', 'slat', 'lev',  ifirstxy, ilastxy, jfirstxy, jlastxy, &
    1739           0 :                1,num_met_levels, wrk3_xy(:,:,blev4:elev4), readvar, gridname='fv_u_stagger',timelevel=recnos(i))
    1740             : 
    1741             :           ! transpose xy -> yz decomposition
    1742           0 :           call transpose_xy2yz_3d( wrk3_xy(:,:,blev3:elev3), met_usi(i)%data(:,:,:), grid )
    1743             : 
    1744           0 :           wrk3_xy = 0._r8
    1745           0 :           met_vsi(i)%data(:,:,:) = 0._r8
    1746             :           call infld(Vname, fids(i), 'slon', 'lat', 'lev',  ifirstxy, ilastxy, jfirstxy, jlastxy, &
    1747           0 :                1,num_met_levels, wrk3_xy(:,:,blev4:elev4), readvar, gridname='fv_v_stagger',timelevel=recnos(i))
    1748             : 
    1749             :           ! transpose xy -> yz decomposition
    1750           0 :           call transpose_xy2yz_3d( wrk3_xy(:,:,blev3:elev3), met_vsi(i)%data(:,:,:), grid )
    1751             : 
    1752             :        else
    1753             : 
    1754             :           ! read into lower portion of the array...
    1755             : 
    1756           0 :           wrk3_xy = 0._r8
    1757           0 :           met_ui(i)%data = 0._r8
    1758             :           call infld(Uname, fids(i), 'lon', 'lat', 'lev',  ifirstxy, ilastxy, jfirstxy, jlastxy, &
    1759           0 :                1,num_met_levels, wrk3_xy(:,:,blev4:elev4), readvar, gridname='fv_centers',timelevel=recnos(i))
    1760             : 
    1761             :           ! transpose xy -> yz decomposition
    1762           0 :           call transpose_xy2yz_3d( wrk3_xy(:,:,blev3:elev3), met_ui(i)%data(:,:,:), grid )
    1763             : 
    1764           0 :           wrk3_xy = 0._r8
    1765           0 :           met_vi(i)%data = 0._r8
    1766             :           call infld(Vname, fids(i), 'lon', 'lat', 'lev',  ifirstxy, ilastxy, jfirstxy, jlastxy, &
    1767           0 :                1,num_met_levels, wrk3_xy(:,:,blev4:elev4), readvar, gridname='fv_centers',timelevel=recnos(i))
    1768             : 
    1769             :           ! transpose xy -> yz decomposition
    1770           0 :           call transpose_xy2yz_3d( wrk3_xy(:,:,blev3:elev3), met_vi(i)%data(:,:,:), grid )
    1771             : 
    1772             :        endif ! met_cell_wall_winds
    1773             : 
    1774             :        call infld(PSname, fids(i), 'lon', 'lat',  ifirstxy, ilastxy, jfirstxy, jlastxy, &
    1775           0 :             wrk2_xy, readvar, gridname='fv_centers', timelevel=recnos(i))
    1776             : 
    1777             :        ! transpose xy -> yz decomposition
    1778           0 :         call transpose_xy2yz_2d( wrk2_xy, met_psi_curr(i)%data, grid )
    1779             : 
    1780             :     enddo
    1781             : 
    1782           0 :     deallocate(tmp_data)
    1783           0 :     deallocate(wrk3_xy)
    1784           0 :     deallocate(wrk2_xy)
    1785             : 
    1786             :     ! 2-D feilds
    1787           0 :     call read_phys_srf_flds( )
    1788             : 
    1789           0 :     if(masterproc) write(iulog,*)'READ_NEXT_METDATA: Read meteorological data '
    1790             : 
    1791           0 :     call t_stopf('MET__read_next_metdata')
    1792             : 
    1793           0 :   end subroutine read_next_metdata
    1794             : 
    1795             : !------------------------------------------------------------------------------
    1796             : !------------------------------------------------------------------------------
    1797           0 :   subroutine read_phys_srf_flds( )
    1798           0 :     use ncdio_atm,          only: infld
    1799             : 
    1800             :     integer :: i, recnos(2)
    1801           0 :     type(file_desc_t) :: fids(2)
    1802             :     logical :: readvar
    1803             : 
    1804           0 :     call find_times( recnos, fids, datatimem, datatimep, curr_mod_time )
    1805           0 :     do i=1,2
    1806             : 
    1807           0 :        call infld(met_shflx_name, fids(i), 'lon', 'lat',  1, pcols, begchunk, endchunk, &
    1808           0 :             met_shflxi(i)%data, readvar, gridname='physgrid',timelevel=recnos(i))
    1809             : 
    1810           0 :        if (has_lhflx) then
    1811             :          call infld('LHFLX', fids(i), 'lon', 'lat',  1, pcols, begchunk, endchunk, &
    1812           0 :               met_lhflxi(i)%data, readvar, gridname='physgrid',timelevel=recnos(i))
    1813             :        end if
    1814             : 
    1815             :        call infld(met_qflx_name, fids(i), 'lon', 'lat',  1, pcols, begchunk, endchunk, &
    1816           0 :             met_qflxi(i)%data, readvar, gridname='physgrid',timelevel=recnos(i))
    1817             :        call infld('TAUX', fids(i), 'lon', 'lat',  1, pcols, begchunk, endchunk, &
    1818           0 :             met_tauxi(i)%data, readvar, gridname='physgrid',timelevel=recnos(i))
    1819             :        call infld('TAUY', fids(i), 'lon', 'lat',  1, pcols, begchunk, endchunk, &
    1820           0 :             met_tauyi(i)%data, readvar, gridname='physgrid',timelevel=recnos(i))
    1821             : 
    1822           0 :        if ( .not.met_srf_feedback ) then
    1823             :           call infld('SNOWH', fids(i), 'lon', 'lat',  1, pcols, begchunk, endchunk, &
    1824           0 :                met_snowhi(i)%data, readvar, gridname='physgrid',timelevel=recnos(i))
    1825             :        endif
    1826             : 
    1827           0 :        if (has_ts) then
    1828             :           call infld('TS', fids(i), 'lon', 'lat',  1, pcols, begchunk, endchunk, &
    1829           0 :                met_tsi(i)%data, readvar, gridname='physgrid',timelevel=recnos(i))
    1830             :        endif
    1831             : 
    1832           0 :        if (met_srf_rad) then
    1833             :           call infld('ASDIR', fids(i), 'lon', 'lat',  1, pcols, begchunk, endchunk, &
    1834           0 :                met_asdiri(i)%data, readvar, gridname='physgrid',timelevel=recnos(i))
    1835             :           call infld('ASDIF', fids(i), 'lon', 'lat',  1, pcols, begchunk, endchunk, &
    1836           0 :                met_asdifi(i)%data, readvar, gridname='physgrid',timelevel=recnos(i))
    1837             :           call infld('ALDIR', fids(i), 'lon', 'lat',  1, pcols, begchunk, endchunk, &
    1838           0 :                met_aldiri(i)%data, readvar, gridname='physgrid',timelevel=recnos(i))
    1839             :           call infld('ALDIF', fids(i), 'lon', 'lat',  1, pcols, begchunk, endchunk, &
    1840           0 :                met_aldifi(i)%data, readvar, gridname='physgrid',timelevel=recnos(i))
    1841             :           call infld('LWUP', fids(i), 'lon', 'lat',  1, pcols, begchunk, endchunk, &
    1842           0 :                met_lwupi(i)%data, readvar, gridname='physgrid',timelevel=recnos(i))
    1843             :        endif
    1844             : 
    1845           0 :        if (met_srf_refs) then
    1846             :           call infld('QREF', fids(i), 'lon', 'lat',  1, pcols, begchunk, endchunk, &
    1847           0 :                met_qrefi(i)%data, readvar, gridname='physgrid',timelevel=recnos(i))
    1848             :           call infld('TREF', fids(i), 'lon', 'lat',  1, pcols, begchunk, endchunk, &
    1849           0 :                met_trefi(i)%data, readvar, gridname='physgrid',timelevel=recnos(i))
    1850             :           call infld('U10', fids(i), 'lon', 'lat',  1, pcols, begchunk, endchunk, &
    1851           0 :                met_u10i(i)%data, readvar, gridname='physgrid',timelevel=recnos(i))
    1852             :        endif
    1853             : 
    1854           0 :        if (met_srf_sst) then
    1855             :           call infld('SST', fids(i), 'lon', 'lat',  1, pcols, begchunk, endchunk, &
    1856           0 :                met_ssti(i)%data, readvar, gridname='physgrid',timelevel=recnos(i))
    1857             :           call infld('ICEFRAC', fids(i), 'lon', 'lat',  1, pcols, begchunk, endchunk, &
    1858           0 :                met_icefraci(i)%data, readvar, gridname='physgrid',timelevel=recnos(i))
    1859             :        endif
    1860             :     enddo
    1861           0 :   end subroutine read_phys_srf_flds
    1862             : 
    1863             : !------------------------------------------------------------------------------
    1864             : !------------------------------------------------------------------------------
    1865           0 :   subroutine interp_phys_srf_flds( )
    1866           0 :     use phys_grid,           only: get_ncols_p
    1867             :     real(r4) :: fact1, fact2
    1868             :     real(r8) :: deltat
    1869             :     integer :: i, c, ncol
    1870             : 
    1871           0 :     deltat = datatimep - datatimem
    1872           0 :     fact1 = (datatimep - curr_mod_time)/deltat
    1873           0 :     fact2 = D1_0-fact1
    1874             : 
    1875             : 
    1876           0 :     do c=begchunk,endchunk
    1877           0 :        ncol = get_ncols_p(c)
    1878           0 :        do i=1,ncol
    1879           0 :           if (has_lhflx) then
    1880           0 :             met_lhflx(i,c) = fact1*met_lhflxi(nm)%data(i,c) + fact2*met_lhflxi(np)%data(i,c)
    1881             :           end if
    1882           0 :           met_shflx(i,c) = fact1*met_shflxi(nm)%data(i,c) + fact2*met_shflxi(np)%data(i,c)
    1883           0 :           met_qflx(i,c) = fact1*met_qflxi(nm)%data(i,c) + fact2*met_qflxi(np)%data(i,c)
    1884           0 :           met_taux(i,c) = fact1*met_tauxi(nm)%data(i,c) + fact2*met_tauxi(np)%data(i,c)
    1885           0 :           met_tauy(i,c) = fact1*met_tauyi(nm)%data(i,c) + fact2*met_tauyi(np)%data(i,c)
    1886             :        enddo
    1887             :     enddo
    1888           0 :     if ( .not.met_srf_feedback ) then
    1889           0 :        do c=begchunk,endchunk
    1890           0 :           ncol = get_ncols_p(c)
    1891           0 :           do i=1,ncol
    1892           0 :              met_snowh(i,c) = fact1*met_snowhi(nm)%data(i,c) + fact2*met_snowhi(np)%data(i,c)
    1893             :           enddo
    1894             :        enddo
    1895             :     endif
    1896           0 :     if (has_ts) then
    1897           0 :        do c=begchunk,endchunk
    1898           0 :           ncol = get_ncols_p(c)
    1899           0 :           do i=1,ncol
    1900           0 :              met_ts(i,c) = fact1*met_tsi(nm)%data(i,c) + fact2*met_tsi(np)%data(i,c)
    1901             :           enddo
    1902             :        enddo
    1903             :     endif
    1904           0 :     if (met_srf_rad) then
    1905           0 :        do c=begchunk,endchunk
    1906           0 :           ncol = get_ncols_p(c)
    1907           0 :           do i=1,ncol
    1908             :              ! Albedo can be fillValue from the meteorology where the is no downwelling
    1909             :              ! solar. However, this changes slowly, so for interpolation use either end-point
    1910             :              ! if nothing is present. If there is no solar, then the albedo won't matter, so
    1911             :              ! should not cause problems.
    1912           0 :              if (met_asdiri(nm)%data(i,c) == srf_fill_value) then
    1913           0 :                 met_asdir(i,c) = met_asdiri(np)%data(i,c)
    1914           0 :              else if (met_asdiri(np)%data(i,c) == srf_fill_value) then
    1915           0 :                 met_asdir(i,c) = met_asdiri(nm)%data(i,c)
    1916             :              else
    1917           0 :                 met_asdir(i,c) = fact1*met_asdiri(nm)%data(i,c) + fact2*met_asdiri(np)%data(i,c)
    1918             :              endif
    1919             : 
    1920           0 :              if (met_asdifi(nm)%data(i,c) == srf_fill_value) then
    1921           0 :                 met_asdif(i,c) = met_asdifi(np)%data(i,c)
    1922           0 :              else if (met_asdifi(np)%data(i,c) == srf_fill_value) then
    1923           0 :                 met_asdif(i,c) = met_asdifi(nm)%data(i,c)
    1924             :              else
    1925           0 :                 met_asdif(i,c) = fact1*met_asdifi(nm)%data(i,c) + fact2*met_asdifi(np)%data(i,c)
    1926             :              endif
    1927             : 
    1928           0 :              if (met_aldiri(nm)%data(i,c) == srf_fill_value) then
    1929           0 :                 met_aldir(i,c) = met_aldiri(np)%data(i,c)
    1930           0 :              else if (met_aldiri(np)%data(i,c) == srf_fill_value) then
    1931           0 :                 met_aldir(i,c) = met_aldiri(nm)%data(i,c)
    1932             :              else
    1933           0 :                 met_aldir(i,c) = fact1*met_aldiri(nm)%data(i,c) + fact2*met_aldiri(np)%data(i,c)
    1934             :              endif
    1935             : 
    1936           0 :              if (met_aldifi(nm)%data(i,c) == srf_fill_value) then
    1937           0 :                 met_aldif(i,c) = met_aldifi(np)%data(i,c)
    1938           0 :              else if (met_aldifi(np)%data(i,c) == srf_fill_value) then
    1939           0 :                 met_aldif(i,c) = met_aldifi(nm)%data(i,c)
    1940             :              else
    1941           0 :                 met_aldif(i,c) = fact1*met_aldifi(nm)%data(i,c) + fact2*met_aldifi(np)%data(i,c)
    1942             :              endif
    1943             : 
    1944           0 :              met_lwup(i,c)  = fact1*met_lwupi(nm)%data(i,c)  + fact2*met_lwupi(np)%data(i,c)
    1945             :           enddo
    1946             :        enddo
    1947             :     endif
    1948           0 :     if (met_srf_refs) then
    1949           0 :        do c=begchunk,endchunk
    1950           0 :           ncol = get_ncols_p(c)
    1951           0 :           do i=1,ncol
    1952           0 :              met_qref(i,c) = fact1*met_qrefi(nm)%data(i,c) + fact2*met_qrefi(np)%data(i,c)
    1953           0 :              met_tref(i,c) = fact1*met_trefi(nm)%data(i,c) + fact2*met_trefi(np)%data(i,c)
    1954           0 :              met_u10(i,c)  = fact1*met_u10i(nm)%data(i,c)  + fact2*met_u10i(np)%data(i,c)
    1955             :           enddo
    1956             :        enddo
    1957             :     endif
    1958           0 :     if (met_srf_sst) then
    1959           0 :        do c=begchunk,endchunk
    1960           0 :           ncol = get_ncols_p(c)
    1961           0 :           do i=1,ncol
    1962             :              ! The sst is fill value over land, which should not change from timestep to
    1963             :              ! timestep, but just in case use the sst value if only one is present.
    1964           0 :              if (met_ssti(nm)%data(i,c) == srf_fill_value) then
    1965           0 :                 met_sst(i,c) = met_ssti(np)%data(i,c)
    1966           0 :              else if (met_ssti(np)%data(i,c) == srf_fill_value) then
    1967           0 :                 met_sst(i,c) = met_ssti(nm)%data(i,c)
    1968             :              else
    1969           0 :                 met_sst(i,c) = fact1*met_ssti(nm)%data(i,c) + fact2*met_ssti(np)%data(i,c)
    1970             :              endif
    1971             : 
    1972           0 :              if (met_icefraci(nm)%data(i,c) == srf_fill_value) then
    1973           0 :                 met_icefrac(i,c) = met_icefraci(np)%data(i,c)
    1974           0 :              else if (met_ssti(np)%data(i,c) == srf_fill_value) then
    1975           0 :                 met_icefrac(i,c) = met_icefraci(nm)%data(i,c)
    1976             :              else
    1977           0 :                 met_icefrac(i,c) = fact1*met_icefraci(nm)%data(i,c) + fact2*met_icefraci(np)%data(i,c)
    1978             :              endif
    1979             :           enddo
    1980             :        enddo
    1981             :     endif
    1982             : 
    1983           0 :   end subroutine interp_phys_srf_flds
    1984             : !------------------------------------------------------------------------------
    1985             : !------------------------------------------------------------------------------
    1986           0 :   subroutine interpolate_metdata(grid)
    1987           0 :     use phys_grid,           only: get_ncols_p
    1988             : #if defined( SPMD )
    1989             :     use mod_comm, only : mp_send4d_ns, mp_recv4d_ns
    1990             : #endif
    1991             : 
    1992             :     implicit none
    1993             :     type (T_FVDYCORE_GRID), intent(in) :: grid
    1994             : 
    1995             :     real(r4) fact1, fact2
    1996             :     real(r4) nfact1, nfact2
    1997             :     real(r8) deltat,deltatn
    1998             :     integer :: i,c,k, ncol
    1999             :     integer :: im, jm, km, jfirst, jlast, kfirst, klast, ng_d, ng_s
    2000             : 
    2001           0 :     im      = grid%im
    2002           0 :     jm      = grid%jm
    2003           0 :     km      = grid%km
    2004           0 :     jfirst  = grid%jfirst
    2005           0 :     jlast   = grid%jlast
    2006           0 :     kfirst  = grid%kfirst
    2007           0 :     klast   = grid%klast
    2008           0 :     ng_d    = grid%ng_d
    2009           0 :     ng_s    = grid%ng_s
    2010             : 
    2011           0 :     call t_startf('MET__interpolate_metdata')
    2012             : 
    2013           0 :     deltat = datatimep - datatimem
    2014           0 :     deltatn = datatimepn - datatimemn
    2015             : 
    2016           0 :     fact1 = (datatimep - curr_mod_time)/deltat
    2017             : !    fact2 = (curr_mod_time - datatimem)/deltat
    2018           0 :     fact2 = D1_0-fact1
    2019             : 
    2020           0 :     nfact1 = (datatimepn - next_mod_time)/deltatn
    2021             : !    nfact2 = (next_mod_time - datatimemn)/deltatn
    2022           0 :     nfact2 = D1_0-nfact1
    2023             : 
    2024           0 :     met_q = 0.0_r8
    2025           0 :     do c=begchunk,endchunk
    2026           0 :        ncol = get_ncols_p(c)
    2027           0 :        do k=1,pver
    2028           0 :           do i=1,ncol
    2029           0 :              met_t(i,k,c) = fact1*met_ti(nm)%data(i,k,c) + fact2*met_ti(np)%data(i,k,c)
    2030           0 :              met_q(i,k,c) = fact1*met_qi(nm)%data(i,k,c) + fact2*met_qi(np)%data(i,k,c)
    2031             :           enddo
    2032             :        enddo
    2033             :     enddo
    2034             : 
    2035           0 :     if (.not. online_test) where (met_q .lt. D0_0) met_q = D0_0
    2036             : 
    2037           0 :     met_ps_next(:,:) =  nfact1*met_psi_next(nm)%data(:,:) + nfact2*met_psi_next(np)%data(:,:)
    2038           0 :     met_ps_curr(:,:) =   fact1*met_psi_curr(nm)%data(:,:) +  fact2*met_psi_curr(np)%data(:,:)
    2039             : 
    2040           0 :     call interp_phys_srf_flds()
    2041             : 
    2042           0 :     if (has_ts) then
    2043           0 :        do c=begchunk,endchunk
    2044           0 :           ncol = get_ncols_p(c)
    2045           0 :           do i=1,ncol
    2046           0 :              met_ts(i,c) = fact1*met_tsi(nm)%data(i,c) + fact2*met_tsi(np)%data(i,c)
    2047             :           enddo
    2048             :        enddo
    2049             :     endif
    2050             : 
    2051           0 :     if ( .not. met_cell_wall_winds ) then
    2052             : 
    2053           0 :        met_u(1:im,jfirst:jlast,kfirst:klast) = fact1*met_ui(nm)%data(1:im,jfirst:jlast,kfirst:klast) &
    2054           0 :                                                  + fact2*met_ui(np)%data(1:im,jfirst:jlast,kfirst:klast)
    2055           0 :        met_v(1:im,jfirst:jlast,kfirst:klast) = fact1*met_vi(nm)%data(1:im,jfirst:jlast,kfirst:klast) &
    2056           0 :                                                  + fact2*met_vi(np)%data(1:im,jfirst:jlast,kfirst:klast)
    2057             : 
    2058             :        ! ghost u,v
    2059             : #if defined( SPMD )
    2060             :        call mp_send4d_ns( grid%commxy, im, jm, km, 1, jfirst, jlast,       &
    2061           0 :             kfirst, klast, ng_d, ng_d, met_u )
    2062             :        call mp_send4d_ns( grid%commxy, im, jm, km, 1, jfirst, jlast,       &
    2063           0 :             kfirst, klast, ng_d, ng_s, met_v )
    2064             :        call mp_recv4d_ns( grid%commxy, im, jm, km, 1, jfirst, jlast,       &
    2065           0 :             kfirst, klast, ng_d, ng_d, met_u )
    2066             :        call mp_recv4d_ns( grid%commxy, im, jm, km, 1, jfirst, jlast,       &
    2067           0 :             kfirst, klast, ng_d, ng_s, met_v )
    2068             : #endif
    2069             : 
    2070             :        ! average to cell walls (vorticity winds)
    2071           0 :        call transfer_windsToWalls(grid)
    2072             :     else
    2073           0 :        met_us(:,jfirst:jlast,kfirst:klast) = fact1*met_usi(nm)%data(:,jfirst:jlast,kfirst:klast) + &
    2074           0 :             fact2*met_usi(np)%data(:,jfirst:jlast,kfirst:klast)
    2075           0 :        met_vs(:,jfirst:jlast,kfirst:klast) = fact1*met_vsi(nm)%data(:,jfirst:jlast,kfirst:klast) + &
    2076           0 :             fact2*met_vsi(np)%data(:,jfirst:jlast,kfirst:klast)
    2077             : 
    2078             :     endif
    2079             : 
    2080             :     ! ghost staggered u,v
    2081             :     ! WS 2006.04.11: not necessary here since it will be done in cd_core
    2082             : 
    2083             : !    write(iulog,*)'INTERPOLATE_METDATA: complete.'
    2084             : 
    2085           0 :     call t_stopf('MET__interpolate_metdata')
    2086             : 
    2087           0 :   end subroutine interpolate_metdata
    2088             : 
    2089             : !-----------------------------------------------------------------------
    2090             : !-----------------------------------------------------------------------
    2091           0 :   subroutine get_dimension( fid, dname, dsize )
    2092             :     implicit none
    2093             :     type(file_desc_t), intent(in) :: fid
    2094             :     character(*), intent(in) :: dname
    2095             :     integer, intent(out) :: dsize
    2096             : 
    2097             :     integer :: dimid, ierr
    2098             : 
    2099           0 :     ierr = pio_inq_dimid( fid, dname, dimid )
    2100           0 :     ierr = pio_inq_dimlen( fid, dimid, dsize )
    2101             : 
    2102           0 :   end subroutine get_dimension
    2103             : 
    2104             : !-----------------------------------------------------------------------
    2105             : !-----------------------------------------------------------------------
    2106           0 :   subroutine open_met_datafile( fname, fileid, times, datapath, check_dims, grid )
    2107             : 
    2108             :     use ioFileMod, only: getfil
    2109             :     use pio,       only: pio_seterrorhandling, PIO_INTERNAL_ERROR, PIO_BCAST_ERROR, PIO_NOERR
    2110             : 
    2111             :     implicit none
    2112             : 
    2113             :     character(*), intent(in) :: fname
    2114             :     type(file_desc_t), intent(inout) :: fileid
    2115             :     real(r8), pointer, intent(inout) :: times(:)
    2116             :     character(*), intent(in) :: datapath
    2117             :     logical, optional, intent(in) :: check_dims
    2118             :     type (T_FVDYCORE_GRID), optional, intent(in) :: grid
    2119             : 
    2120             :     character(len=256) :: filepath
    2121             :     character(len=256) :: filen
    2122             :     character(len=*), parameter :: subname = 'open_met_datafile'
    2123             :     integer :: year, month, day, dsize, i, timesize
    2124             :     integer :: dateid,secid
    2125           0 :     integer, allocatable , dimension(:) :: dates, datesecs
    2126             :     integer :: ierr
    2127             :     integer :: im, jm, km
    2128             :     integer :: varid
    2129             : 
    2130             :     !
    2131             :     ! open file and get fileid
    2132             :     !
    2133           0 :     if (len_trim( datapath )>0) then
    2134           0 :        filepath = trim(datapath)//'/'//trim(fname)
    2135             :     else
    2136           0 :        filepath = trim(fname)
    2137             :     endif
    2138           0 :     call getfil( filepath, filen, 0 )
    2139             : 
    2140           0 :     call cam_pio_openfile( fileid, filen, 0 )
    2141             : 
    2142           0 :     call pio_seterrorhandling(fileid, PIO_BCAST_ERROR)
    2143             : 
    2144           0 :     ierr = pio_inq_varid( fileid, 'TS', varid )
    2145           0 :     has_ts = ierr==PIO_NOERR
    2146             : 
    2147           0 :     ierr = pio_inq_varid( fileid, 'LHFLX', varid )
    2148           0 :     has_lhflx = ierr==PIO_NOERR
    2149             : 
    2150           0 :     call pio_seterrorhandling(fileid, PIO_INTERNAL_ERROR)
    2151             : 
    2152           0 :     if (masterproc) write(iulog,*) 'open_met_datafile: ',trim(filen)
    2153             : 
    2154           0 :     call get_dimension( fileid, 'time', timesize )
    2155             : 
    2156           0 :     if ( associated(times) ) deallocate(times)
    2157           0 :     allocate( times(timesize) )
    2158             : 
    2159           0 :     allocate( dates(timesize) )
    2160           0 :     allocate( datesecs(timesize) )
    2161             : 
    2162           0 :     ierr = pio_inq_varid( fileid, 'date',    dateid  )
    2163           0 :     ierr = pio_inq_varid( fileid, 'datesec', secid  )
    2164             : 
    2165           0 :     ierr = pio_get_var( fileid, dateid, dates )
    2166           0 :     ierr = pio_get_var( fileid, secid,  datesecs  )
    2167             : 
    2168           0 :     do i=1,timesize
    2169           0 :        year = dates(i) / 10000
    2170           0 :        month = mod(dates(i),10000)/100
    2171           0 :        day = mod(dates(i),100)
    2172           0 :        times(i) = get_time_float( year, month, day, datesecs(i) )
    2173             :     enddo
    2174             : 
    2175           0 :     deallocate( dates )
    2176           0 :     deallocate( datesecs )
    2177             : 
    2178             : 
    2179             : !
    2180             : ! check that the data dim sizes match models dimensions
    2181             : !
    2182           0 :     if (present(check_dims) .and. present(grid)) then
    2183           0 :        im         = grid%im
    2184           0 :        jm         = grid%jm
    2185           0 :        km         = grid%km
    2186             : 
    2187           0 :        if (check_dims) then
    2188             : 
    2189           0 :           call get_dimension( fileid, 'lon', dsize )
    2190           0 :           if (dsize /= im) then
    2191           0 :              write(iulog,*)'open_met_datafile: lonsiz=',dsize,' must = ',im
    2192           0 :              call endrun (subname // ':: ERRROR - check atm.log for error message')
    2193             :           endif
    2194           0 :           call get_dimension( fileid, 'lat', dsize )
    2195           0 :           if (dsize /= jm) then
    2196           0 :              write(iulog,*)'open_met_datafile: latsiz=',dsize,' must = ',jm
    2197           0 :              call endrun (subname // ':: ERRROR - check atm.log for error message')
    2198             :           endif
    2199           0 :           call get_dimension( fileid, 'lev', dsize )
    2200           0 :           met_levels = min( dsize, km )
    2201           0 :           num_met_levels = dsize
    2202             :        endif
    2203             :     endif
    2204             : 
    2205           0 :     if (met_srf_rad) then
    2206           0 :        ierr = pio_inq_varid( fileid, 'ASDIR', varid )
    2207           0 :        ierr = pio_get_att( fileid, varid, '_FillValue', srf_fill_value)
    2208             :     endif
    2209             : 
    2210           0 :   end subroutine open_met_datafile
    2211             : 
    2212             : !------------------------------------------------------------------------------
    2213             : !------------------------------------------------------------------------------
    2214           0 :   function get_time_float( year, month, day, sec )
    2215             : 
    2216             : ! returns float representation of time -- number of days
    2217             : ! since 1 jan 0001 00:00:00.000
    2218             : 
    2219             :     implicit none
    2220             : 
    2221             :     integer, intent(in) :: year, month, day
    2222             :     integer, intent(in) :: sec
    2223             :     real(r8) :: get_time_float
    2224             : 
    2225             : ! ref date is 1 jan 0001
    2226             : 
    2227             :     integer :: refyr, refmn, refdy
    2228             :     real(r8) :: refsc, fltdy
    2229             :     integer :: doy(12)
    2230             : 
    2231             : !              jan feb mar apr may jun jul aug sep oct nov dec
    2232             : !              31  28  31  30  31  30  31  31  31  31  30  31
    2233             :     data doy /  1, 32, 60, 91,121,152,182,213,244,274,305,335 /
    2234             : 
    2235           0 :     refyr = 1
    2236           0 :     refmn = 1
    2237           0 :     refdy = 1
    2238           0 :     refsc = D0_0
    2239             : 
    2240           0 :     if ( timemgr_is_caltype(trim(shr_cal_gregorian))) then
    2241           0 :        fltdy = greg2jday(year, month, day) - greg2jday(refyr,refmn,refdy)
    2242             :     else ! assume no_leap (all years are 365 days)
    2243             :        fltdy = (year - refyr)*days_per_non_leapyear + &
    2244           0 :                (doy(month)-doy(refmn)) + &
    2245           0 :                (day-refdy)
    2246             :     endif
    2247             : 
    2248           0 :     get_time_float = fltdy + ((sec-refsc)/seconds_per_day)
    2249             : 
    2250           0 :   endfunction get_time_float
    2251             : 
    2252             : !-----------------------------------------------------------------------
    2253             : !       ... Return Julian day number given Gregorian date.
    2254             : !
    2255             : ! Algorithm from Hatcher,D.A., Simple Formulae for Julian Day Numbers
    2256             : ! and Calendar Dates, Q.Jl.R.astr.Soc. (1984) v25, pp 53-55.
    2257             : !-----------------------------------------------------------------------
    2258           0 :   function greg2jday( year, month, day )
    2259             : 
    2260             :     implicit none
    2261             : 
    2262             :     integer, intent(in) :: year, month, day
    2263             :     integer :: greg2jday
    2264             : 
    2265             :     !-----------------------------------------------------------------------
    2266             :     !   ... Local variables
    2267             :     !-----------------------------------------------------------------------
    2268             :     integer :: ap, mp
    2269             :     integer :: y, d, n, g
    2270             : 
    2271             :     !-----------------------------------------------------------------------
    2272             :     !           ... Modify year and month numbers
    2273             :     !-----------------------------------------------------------------------
    2274           0 :     ap = year - (12 - month)/10
    2275           0 :     mp = MOD( month-3,12 )
    2276           0 :     if( mp < 0 ) then
    2277           0 :        mp = mp + 12
    2278             :     end if
    2279             : 
    2280             :     !-----------------------------------------------------------------------
    2281             :     !           ... Julian day
    2282             :     !-----------------------------------------------------------------------
    2283           0 :     y = INT( days_per_year*( ap + 4712 ) )
    2284           0 :     d = INT( days_per_month*mp + D0_5 )
    2285           0 :     n = y + d + day  + 59
    2286           0 :     g = INT( D0_75*INT( ap/100 + 49 ) ) - 38
    2287           0 :     greg2jday = n - g
    2288             : 
    2289           0 :   end function greg2jday
    2290             : 
    2291             : !-----------------------------------------------------------------------
    2292             : !-----------------------------------------------------------------------
    2293           0 :   subroutine set_met_rlx( )
    2294             : 
    2295             :     use pmgrid, only: plev
    2296             :     use hycoef, only: hypm, ps0
    2297             : 
    2298             :     integer :: k, k_cnt, k_top
    2299             :     real(r8), parameter :: h0 = 7._r8        ! scale height (km)
    2300             :     real(r8), parameter :: hsec = 3600._r8   ! seconds per hour
    2301             :     real(r8) :: p_top, p_bot
    2302             :     real(r8) :: p_bot_top, p_bot_bot
    2303             :     real(r8) ::  met_max_rlxdt,  dtime_hrs
    2304             : 
    2305             : 996 format( 'set_met_rlx: ',a15, I10 )
    2306             : 997 format( 'set_met_rlx: ',a15, E10.2 )
    2307             : 998 format( 'set_met_rlx: ',a15, PLEV(E10.2))
    2308             : 999 format( 'set_met_rlx: ',a15, PLEV(F10.5))
    2309             : 993 format( 'set_met_rlx: ',a25, E10.2 )
    2310             : 991 format( 'set_met_rlx: ',a25, I4, E10.2, E10.2 )
    2311             : 
    2312           0 :     met_rlx(:) = 999._r8
    2313             : 
    2314           0 :     dtime_hrs = get_step_size()/hsec ! hours
    2315             : 
    2316           0 :     if (met_rlx_time > dtime_hrs) then
    2317           0 :        met_max_rlxdt = dtime_hrs/met_rlx_time
    2318           0 :     elseif (met_rlx_time < 0._r8) then
    2319           0 :        met_max_rlxdt = 0._r8
    2320             :     else
    2321           0 :        met_max_rlxdt = 1._r8
    2322             :     endif
    2323             : 
    2324           0 :     if (masterproc) then
    2325           0 :        write(iulog,fmt=993) ' met_rlx_time in hrs= ', met_rlx_time
    2326           0 :        write(iulog,fmt=993) ' met_max_rlxdt in % = ', met_max_rlxdt*100._r8
    2327             :     endif
    2328             : 
    2329           0 :     p_top = ps0 * exp( - met_rlx_top/h0 )
    2330           0 :     p_bot = ps0 * exp( - met_rlx_bot/h0 )
    2331             : 
    2332           0 :     p_bot_top = ps0 * exp( - met_rlx_bot_top/h0 )
    2333           0 :     p_bot_bot = ps0 * exp( - met_rlx_bot_bot/h0 )
    2334             : 
    2335           0 :     if (masterproc) then
    2336           0 :        write(iulog,fmt=997) 'p_top = ',p_top
    2337           0 :        write(iulog,fmt=997) 'p_bot = ',p_bot
    2338           0 :        write(iulog,fmt=997) 'p_bot_top = ',p_bot_top
    2339           0 :        write(iulog,fmt=997) 'p_bot_bot = ',p_bot_bot
    2340             :     endif
    2341             : 
    2342           0 :     if ( p_bot < hypm( pver-met_levels+1 ) .and. ( met_levels < pver ) ) then
    2343           0 :        call endrun(  'set_met_rlx: met_rlx_bot is too high ' )
    2344             :     endif
    2345             : 
    2346           0 :     met_rlx = 0._r8
    2347             : 
    2348           0 :     where (p_top <= hypm .and. hypm < p_bot)
    2349             :        met_rlx = 999._r8
    2350             :     endwhere
    2351             : 
    2352           0 :     where( p_bot <= hypm .and. hypm <= p_bot_top )
    2353             :        met_rlx = met_max_rlxdt
    2354             :     endwhere
    2355             : 
    2356           0 :     where(p_bot_top < hypm .and. hypm <= p_bot_bot)
    2357             :        met_rlx = -999._r8
    2358             :     endwhere
    2359             : 
    2360           0 :     if ( any( met_rlx(:) /= met_max_rlxdt) ) then
    2361           0 :        k_top = max(plev - met_levels, 1)
    2362             : 
    2363             :        ! ramp region at model top
    2364           0 :        k_cnt = count(p_top < hypm .and. hypm < p_bot)
    2365           0 :        if (k_cnt > 0) then
    2366             :           k_top = max(plev - met_levels, 1)
    2367           0 :           do while ( met_rlx(k_top) /= 999._r8 )
    2368           0 :              k_top = k_top + 1
    2369           0 :              if ( k_top == pver ) then
    2370           0 :                 call endrun ( 'set_met_rlx: cannot find ramped region ')
    2371             :              endif
    2372             :           enddo
    2373             : 
    2374           0 :           met_rlx(1:k_top) = 0._r8
    2375             : 
    2376           0 :           k_cnt = count( met_rlx == 999._r8 )
    2377             : 
    2378           0 :           if (masterproc) then
    2379           0 :              write(iulog,*) 'top of model ramped region:'
    2380           0 :              write(iulog,fmt=996) 'k_cnt = ',k_cnt
    2381           0 :              write(iulog,fmt=996) 'k_top = ',k_top
    2382             :           endif
    2383             : 
    2384           0 :           do k = k_top,k_top+k_cnt
    2385           0 :              met_rlx(k) = met_max_rlxdt*real( k - k_top ) / real(k_cnt)
    2386             :           enddo
    2387             :        endif
    2388             : 
    2389             :        ! ramp region at model bottom
    2390           0 :        k_cnt = count(p_bot_top < hypm .and. hypm < p_bot_bot)
    2391           0 :        if (k_cnt > 0) then
    2392           0 :           k_top = max(plev - met_levels, 1)
    2393           0 :           do while ( met_rlx(k_top) /= -999._r8 )
    2394           0 :              k_top = k_top + 1
    2395           0 :              if ( k_top == pver ) then
    2396           0 :                 call endrun ( 'set_met_rlx: cannot find ramped region ')
    2397             :              endif
    2398             :           enddo
    2399             : 
    2400           0 :           if (masterproc) then
    2401           0 :              write(iulog,*) 'bottom of model ramped region:'
    2402           0 :              write(iulog,fmt=996) 'k_cnt = ',k_cnt
    2403           0 :              write(iulog,fmt=996) 'k_top = ',k_top
    2404             :           endif
    2405             : 
    2406           0 :           do k = k_top,k_top+k_cnt-1
    2407           0 :              met_rlx(k) = met_max_rlxdt*(1._r8 - real( k - k_top +1) / real(k_cnt))
    2408             :           enddo
    2409             :        endif
    2410             : 
    2411             :     endif
    2412             : 
    2413           0 :     if (masterproc) then
    2414           0 :        write(iulog,fmt=996) '  met_levels = ',met_levels
    2415           0 :        write(iulog,fmt=996) 'non-zero terms:',count( met_rlx /= 0._r8 )
    2416             :     endif
    2417             : 
    2418           0 :     if ( met_levels < count( met_rlx /= 0._r8 ) ) then
    2419           0 :        call endrun('set_met_rlx: met_rlx_top is too high for the meteorology data')
    2420             :     endif
    2421             : 
    2422           0 :     if (masterproc) then
    2423           0 :       do k = max(plev - met_levels, 1), pver
    2424           0 :          write(iulog,fmt=991) 'press levels, met_rlx = ', k, hypm(k), met_rlx(k)
    2425             :       end do
    2426             :     endif
    2427             : 
    2428           0 :     if ( any( (met_rlx > 1._r8) .or. (met_rlx < 0._r8) ) ) then
    2429           0 :        if (masterproc) then
    2430           0 :           write(iulog,fmt=993) 'set_met_rlx:   dtime_hrs in hrs = ', dtime_hrs
    2431           0 :           write(iulog,fmt=993) 'set_met_rlx:met_rlx_time in hrs = ', met_rlx_time
    2432           0 :           write(iulog,fmt=993) 'set_met_rlx:      met_max_rlxdt = ', met_max_rlxdt
    2433           0 :           write(iulog,fmt=999) 'set_met_rlx:            met_rlx = ', met_rlx
    2434             :        endif
    2435           0 :        call endrun('Offline meteorology relaxation function not set correctly.')
    2436             :     endif
    2437             : 
    2438           0 :   end subroutine set_met_rlx
    2439             : 
    2440           0 : end module metdata

Generated by: LCOV version 1.14