fv_arrays_nlm.F90

Go to the documentation of this file.

!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of fvGFS.                                       *
!*                                                                     *
!* fvGFS is free software; you can redistribute it and/or modify it    *
!* and are expected to follow the terms of the GNU General Public      *
!* License as published by the Free Software Foundation; either        *
!* version 2 of the License, or (at your option) any later version.    *
!*                                                                     *
!* fvGFS is distributed in the hope that it will be useful, but        *
!* WITHOUT ANY WARRANTY; without even the implied warranty of          *
!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU   *
!* General Public License for more details.                            *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************
module fv_arrays_nlm_mod
#include <fms_platform.h>
  use mpp_domains_mod,       only: domain2d
  use fms_io_mod,            only: restart_file_type
  use time_manager_mod,      only: time_type
  use horiz_interp_type_mod, only: horiz_interp_type
  use mpp_domains_mod,       only: nest_domain_type
  use mpp_mod,               only: mpp_broadcast
  use platform_mod,          only: r8_kind
  use, intrinsic :: iso_fortran_env, only: real64, real32
  public

  integer, public, parameter :: r_grid = r8_kind
  integer, parameter :: real4 = real32
  integer, parameter :: real8 = real64
#ifdef SINGLE_FV
  integer, parameter :: fvprc = real4  ! Full Build Precision for the model
#else
  integer, parameter :: fvprc = real8  ! Full Build Precision for the model
#endif

  !Several 'auxiliary' structures are introduced here. These are for
  ! the internal use by certain modules, and although fv_atmos_type
  !  contains one of each of these structures all memory management
  !   is performed by the module in question.


  integer, parameter:: max_step = 1000
!--- MAY NEED TO TEST THIS
#ifdef OVERLOAD_R4
  real, parameter:: real_big = 1.e8    ! big enough to cause blowup if used
#else
  real, parameter:: real_big = 1.e30   ! big enough to cause blowup if used
#endif
  type fv_diag_type


 integer ::id_ps, id_slp, id_ua, id_va, id_pt, id_omga, id_vort,  &
           id_tm, id_pv, id_zsurf, id_oro, id_sgh, id_divg, id_w, &
           id_ke, id_te, id_zs, id_ze, id_mq, id_vorts, id_us, id_vs,    &
           id_tq, id_rh, id_c15, id_c25, id_c35, id_c45,          &
                         id_f15, id_f25, id_f35, id_f45, id_ctp,  &
           id_ppt, id_ts, id_tb, id_ctt, id_pmask, id_pmaskv2,    &
           id_delp, id_delz, id_zratio, id_ws, id_iw, id_lw,      &
           id_pfhy, id_pfnh,                                      &
           id_qn, id_qn200, id_qn500, id_qn850, id_qp, id_mdt, id_qdt, id_aam, id_amdt, &
           id_acly, id_acl, id_acl2, id_dbz, id_maxdbz, id_basedbz, id_dbz4km

! Selected p-level fields from 3D variables:
 integer :: id_vort200, id_vort500, id_w500, id_w700
 integer :: id_vort850, id_w850, id_x850, id_srh, id_srh25, id_srh01, &
            id_uh03, id_uh25, id_theta_e,  &
            id_w200, id_s200, id_sl12, id_sl13, id_w5km, id_rain5km, id_w2500m
! NGGPS 31-level diag
 integer, allocatable :: id_u(:), id_v(:), id_t(:), id_h(:), id_q(:), id_omg(:)

 integer:: id_u_plev, id_v_plev, id_t_plev, id_h_plev, id_q_plev, id_omg_plev
! IPCC diag
 integer :: id_rh10,  id_rh50,  id_rh100, id_rh200,  id_rh250, id_rh300, &
            id_rh500, id_rh700, id_rh850, id_rh925,  id_rh1000

 integer :: id_rh1000_cmip, id_rh925_cmip, id_rh850_cmip, id_rh700_cmip, id_rh500_cmip, &
            id_rh300_cmip,  id_rh250_cmip, id_rh100_cmip, id_rh50_cmip,  id_rh10_cmip

 integer :: id_hght
 integer :: id_u100m, id_v100m, id_w100m

     ! For initial conditions:
     integer ic_ps, ic_ua, ic_va, ic_ppt
     integer ic_sphum
     integer, allocatable :: id_tracer(:)
! ESM requested diagnostics  -  dry mass/volume mixing ratios
 integer, allocatable :: id_tracer_dmmr(:)
 integer, allocatable :: id_tracer_dvmr(:)
 real,    allocatable :: w_mr(:)

     real, allocatable :: phalf(:)
     real, allocatable :: zsurf(:,:)
     real, allocatable :: zxg(:,:)
     real, allocatable :: pt1(:)


     logical :: initialized = .false.
     real  sphum, liq_wat, ice_wat       ! GFDL physics
     real  rainwat, snowwat, graupel

     real :: efx(max_step), efx_sum, efx_nest(max_step), efx_sum_nest, mtq(max_step), mtq_sum
     integer :: steps

  end type fv_diag_type


  !fv_grid_type is made up of grid-dependent information from fv_grid_tools and fv_grid_utils.
  ! It should not contain any user options (that goes in a different structure) nor data which
  ! is altered outside of those two modules.
  type fv_grid_type
     real(kind=R_GRID), allocatable, dimension(:,:,:) :: grid_64, agrid_64
     real(kind=R_GRID), allocatable, dimension(:,:) :: area_64, area_c_64
     real(kind=R_GRID), allocatable, dimension(:,:) :: sina_64, cosa_64
     real(kind=R_GRID), allocatable, dimension(:,:) :: dx_64, dy_64
     real(kind=R_GRID), allocatable, dimension(:,:) :: dxc_64, dyc_64
     real(kind=R_GRID), allocatable, dimension(:,:) :: dxa_64, dya_64

     real, allocatable, dimension(:,:,:) :: grid, agrid
     real, allocatable, dimension(:,:) :: area, area_c
     real, allocatable, dimension(:,:) :: rarea, rarea_c     

     real, allocatable, dimension(:,:) :: sina, cosa
     real, allocatable, dimension(:,:,:) :: e1,e2
     real, allocatable, dimension(:,:) :: dx, dy
     real, allocatable, dimension(:,:) :: dxc, dyc
     real, allocatable, dimension(:,:) :: dxa, dya
     real, allocatable, dimension(:,:) :: rdx, rdy
     real, allocatable, dimension(:,:) :: rdxc, rdyc
     real, allocatable, dimension(:,:) :: rdxa, rdya

     ! Scalars:
     real(kind=R_GRID), allocatable :: edge_s(:)
     real(kind=R_GRID), allocatable :: edge_n(:)
     real(kind=R_GRID), allocatable :: edge_w(:)
     real(kind=R_GRID), allocatable :: edge_e(:)
     ! Vector:
     real(kind=R_GRID), allocatable :: edge_vect_s(:)
     real(kind=R_GRID), allocatable :: edge_vect_n(:)
     real(kind=R_GRID), allocatable :: edge_vect_w(:)
     real(kind=R_GRID), allocatable :: edge_vect_e(:)
     ! scalar:
     real(kind=R_GRID), allocatable :: ex_s(:)
     real(kind=R_GRID), allocatable :: ex_n(:)
     real(kind=R_GRID), allocatable :: ex_w(:)
     real(kind=R_GRID), allocatable :: ex_e(:)

     real, allocatable :: l2c_u(:,:), l2c_v(:,:)
     ! divergence Damping:
     real, allocatable :: divg_u(:,:), divg_v(:,:)    !
     ! del6 diffusion:
     real, allocatable :: del6_u(:,:), del6_v(:,:)    !
     ! Cubed_2_latlon:
     real, allocatable :: a11(:,:)
     real, allocatable :: a12(:,:)
     real, allocatable :: a21(:,:)
     real, allocatable :: a22(:,:)
     ! latlon_2_cubed:
     real, allocatable :: z11(:,:)
     real, allocatable :: z12(:,:)
     real, allocatable :: z21(:,:)
     real, allocatable :: z22(:,:)

!    real, allocatable :: w00(:,:)

     real, allocatable :: cosa_u(:,:)
     real, allocatable :: cosa_v(:,:)
     real, allocatable :: cosa_s(:,:)
     real, allocatable :: sina_u(:,:)
     real, allocatable :: sina_v(:,:)
     real, allocatable :: rsin_u(:,:)
     real, allocatable :: rsin_v(:,:)
     real, allocatable ::  rsina(:,:)
     real, allocatable ::  rsin2(:,:)
     real(kind=R_GRID), allocatable :: ee1(:,:,:)
     real(kind=R_GRID), allocatable :: ee2(:,:,:)
     real(kind=R_GRID), allocatable :: ec1(:,:,:)
     real(kind=R_GRID), allocatable :: ec2(:,:,:)
     real(kind=R_GRID), allocatable :: ew(:,:,:,:)
     real(kind=R_GRID), allocatable :: es(:,:,:,:)


     !- 3D Super grid to contain all geometrical factors --
     ! the 3rd dimension is 9
     real, allocatable :: sin_sg(:,:,:)
     real, allocatable :: cos_sg(:,:,:)
     !--------------------------------------------------

     ! Unit Normal vectors at cell edges:
     real(kind=R_GRID), allocatable :: en1(:,:,:)
     real(kind=R_GRID), allocatable :: en2(:,:,:)

     ! Extended Cubed cross-edge winds
     real, allocatable :: eww(:,:)
     real, allocatable :: ess(:,:)

     ! Unit vectors for lat-lon grid
     real(kind=R_GRID), allocatable :: vlon(:,:,:), vlat(:,:,:)
     real, allocatable :: fc(:,:), f0(:,:)

     integer, dimension(:,:,:), allocatable :: iinta, jinta, iintb, jintb
  
     !Scalar data
     
     integer :: npx_g, npy_g, ntiles_g ! global domain

     real(kind=R_GRID) :: global_area
     logical :: g_sum_initialized = .false. !Not currently used but can be useful
     logical:: sw_corner, se_corner, ne_corner, nw_corner

     real(kind=R_GRID) :: da_min, da_max, da_min_c, da_max_c

     real  :: acapn, acaps
     real  :: globalarea  ! total Global Area
     
     logical :: latlon = .false.
     logical :: cubed_sphere = .false.
     logical :: have_south_pole = .false.
     logical :: have_north_pole = .false.
     logical :: stretched_grid = .false.

     logical :: square_domain = .false.


     !! Convenience pointers

     integer, pointer :: grid_type
     logical, pointer :: nested

  end type fv_grid_type

  type fv_flags_type

     !! FOR EACH VARIABLE IN FV_FLAGS:
     !! 1. Must be defined here:
     !! 2. Must be broadcast in fv_atmos_data
     !! 3. If a namelist entry, a pointer must
     !!    be defined and associated in fv_control
     !! 4. Must NOT appear in fv_current_grid_mod.
     !!    (this module will soon be removed)
     !! 5. Must be referenced through Atm%flagstruct,
     !!    not Atm%, unless a convenience
     !!    pointer is defined

!-----------------------------------------------------------------------
! Grid descriptor file setup
!-----------------------------------------------------------------------
   character(len=80) :: grid_name = 'Gnomonic'
   character(len=120):: grid_file = 'Inline'
  integer      :: grid_type = 0     ! -1: read from file; 0: ED Gnomonic
!                                    !  0: the "true" equal-distance Gnomonic grid
!                                    !  1: the traditional equal-distance Gnomonic grid
!                                    !  2: the equal-angular Gnomonic grid
!                                    !  3: the lat-lon grid -- to be implemented
!                                    !  4: double periodic boundary condition on Cartesian grid
!                                    !  5: channel flow on Cartesian grid
!  -> moved to grid_tools

! Momentum (or KE) options:
   integer :: hord_mt = 9    ! the best option for Gnomonic grids  
   integer :: kord_mt = 8    ! vertical mapping option for (u,v)
   integer :: kord_wz = 8    ! vertical mapping option for w

! Vorticity & w transport options:
   integer :: hord_vt = 9    ! 10 not recommended (noisy case-5) 

! Heat & air mass (delp) transport options:
   integer :: hord_tm = 9    ! virtual potential temperature
   integer :: hord_dp = 9    ! delp (positive definite)
   integer :: kord_tm = 8    !

! Tracer transport options:
   integer :: hord_tr = 12   !11: PPM mono constraint (Lin 2004); fast 
                             !12: Huynh 2nd constraint (Lin 2004) +
                             !    positive definite (Lin & Rood 1996); slower
                             !>12: positive definite only (Lin & Rood 1996); fastest
   integer :: kord_tr = 8    ! 
   real    :: scale_z = 0.   ! diff_z = scale_z**2 * 0.25
   real    :: w_max = 75.    ! max w (m/s) threshold for hydostatiic adjustment 
   real    :: z_min = 0.05   ! min ratio of dz_nonhydrostatic/dz_hydrostatic

   integer :: nord=1         ! 0: del-2, 1: del-4, 2: del-6, 3: del-8 divergence damping
                             ! Alternative setting for high-res: nord=1; d4_bg = 0.075
   integer :: nord_tr=0      ! 0: del-2, 1: del-4, 2: del-6
   real    :: dddmp = 0.0    ! coefficient for del-2 divergence damping (0.2)
                             ! for C90 or lower: 0.2
   real    :: d2_bg = 0.0    ! coefficient for background del-2 divergence damping
   real    :: d4_bg = 0.16   ! coefficient for background del-4(6) divergence damping
                             ! for stability, d4_bg must be <=0.16 if nord=3
   real    :: vtdm4 = 0.0    ! coefficient for del-4 vorticity damping
   real    :: trdm2 = 0.0    ! coefficient for del-2 tracer damping
   real    :: d2_bg_k1 = 4.         ! factor for d2_bg (k=1)
   real    :: d2_bg_k2 = 2.         ! factor for d2_bg (k=2)
   real    :: d2_divg_max_k1 = 0.15 ! d2_divg max value (k=1)
   real    :: d2_divg_max_k2 = 0.08 ! d2_divg max value (k=2)
   real    :: damp_k_k1 = 0.2       ! damp_k value (k=1)
   real    :: damp_k_k2 = 0.12      ! damp_k value (k=2)

! Additional (after the fact) terrain filter (to further smooth the terrain after cold start)
   integer ::    n_zs_filter=0      !  number of application of the terrain filter
   integer :: nord_zs_filter=4      !  use del-2 (2) OR del-4 (4)
   logical :: full_zs_filter=.false.! perform full filtering of topography (in external_ic only )

   logical :: consv_am  = .false.   ! Apply Angular Momentum Correction (to zonal wind component)
   logical :: do_sat_adj= .false.   ! 
   logical :: do_f3d    = .false.   ! 
   logical :: no_dycore = .false.   ! skip the dycore
   logical :: convert_ke = .false. 
   logical :: do_vort_damp = .false. 
   logical :: use_old_omega = .true. 
! PG off centering:
   real    :: beta  = 0.0    ! 0.5 is "neutral" but it may not be stable
#ifdef SW_DYNAMICS
   integer :: n_zfilter = 0  ! Number of layers at the top of the atmosphere for dz filter
   integer :: n_sponge = 0   ! Number of sponge layers at the top of the atmosphere
   real    :: d_ext = 0.    
   integer :: nwat  = 0      ! Number of water species
   logical :: warm_start = .false. 
   logical :: inline_q = .true.
   logical :: adiabatic = .true.     ! Run without physics (full or idealized).
#else
   integer :: n_zfilter = 0  ! Number of layers at the top of the atmosphere for dz filter
   integer :: n_sponge = 1   ! Number of sponge layers at the top of the atmosphere
   real    :: d_ext = 0.02   ! External model damping (was 0.02)
   integer :: nwat  = 3      ! Number of water species
   logical :: warm_start = .true. 
                             ! Set to .F. if cold_start is desired (including terrain generation)
   logical :: inline_q = .false.
   logical :: adiabatic = .false.     ! Run without physics (full or idealized).
#endif
!-----------------------------------------------------------
! Grid shifting, rotation, and the Schmidt transformation:
!-----------------------------------------------------------
   real :: shift_fac   =  18.   ! shift west by 180/shift_fac = 10 degrees
! Defaults for Schmidt transformation:
   logical :: do_schmidt = .false. 
   real(kind=R_GRID) :: stretch_fac =   1.   ! No stretching
   real(kind=R_GRID) :: target_lat  = -90.   ! -90: no grid rotation 
   real(kind=R_GRID) :: target_lon  =   0.   ! 

!-----------------------------------------------------------------------------------------------
! Example #1a: US regional climate simulation, center located over Oklahoma city: (262.4, 35.4)
!              stretching factor: 2.5
! Example #1b: US Hurricane model, center over Miami: (279.7, 25.8)
!              stretching factor: 3-5
! Example #2a: South-East Asia Regional Climate H*** (SERACH), Central Taiwan: (121.0, 23.5)
! Example #2b: Typhoon Model: (124.0, 22.0)
!              stretching factor: 5-10
!-----------------------------------------------------------------------------------------------

   logical :: reset_eta = .false. 
   real    :: p_fac  = 0.05
   real    :: a_imp  = 0.75  ! Off center parameter for the implicit solver [0.5,1.0]
   integer :: n_split = 0    ! Number of time splits for the lagrangian dynamics
                             ! Default = 0 (automatic computation of best value)
   integer :: m_split = 0    ! Number of time splits for Riemann solver
   integer :: k_split = 1    ! Number of time splits for Remapping

   logical :: use_logp = .false.

!            For doubly periodic domain with sim_phys
!                     5km        150         20 (7.5 s)  2
!
!                     Estimates for Gnomonic grids:
            !===================================================
            !        dx (km)    dt (sc)    n_split    m_split
            !===================================================
            ! C1000:  ~10        150         16          3
            ! C2000:   ~5         90         18 (5 s)    2
            !===================================================
! The nonhydrostatic algorithm is described in Lin 2006, QJ, (submitted)
! C2000 should easily scale to at least 6 * 100 * 100 = 60,000 CPUs  
! For a 1024 system: try 6 x 13 * 13 = 1014 CPUs
  
   integer :: q_split = 0    ! Number of time splits for tracer transport

   integer :: print_freq = 0 ! Print max/min of selected fields
                             ! 0: off
                             ! positive n: every n hours
                             ! negative n: every time step

!------------------------------------------
! Model Domain parameters
!------------------------------------------
   integer :: npx                     ! Number of Grid Points in X- dir
   integer :: npy                     ! Number of Grid Points in Y- dir
   integer :: npz                     ! Number of Vertical Levels
   integer :: npz_rst = 0             ! Original Vertical Levels (in the restart)
                                      ! 0: no change (default)
   integer :: ncnst = 0               ! Number of advected consituents
   integer :: pnats = 0               ! Number of non-advected consituents
   integer :: dnats = 0               ! Number of non-advected consituents (as seen by dynamics)
   integer :: ntiles = 1                 ! Number or tiles that make up the Grid 
   integer :: ndims = 2     ! Lat-Lon Dims for Grid in Radians
   integer :: nf_omega  = 1           ! Filter omega "nf_omega" times
   integer :: fv_sg_adj = -1          ! Perform grid-scale dry adjustment if > 0
                                      ! Relaxzation time  scale (sec) if positive
   integer :: na_init = 0             ! Perform adiabatic initialization
   real    :: p_ref = 1.e5
   real    :: dry_mass = 98290.
   integer :: nt_prog = 0
   integer :: nt_phys = 0
   real    :: tau_h2o = 0.            ! Time scale (days) for ch4_chem

   real    :: delt_max = 1.           ! limiter for dissipative heating rate
                                      ! large value (~1) essentially imposes no limit
   real    :: d_con = 0.
   real    :: ke_bg = 0.              ! background KE production (m^2/s^3) over a small step
                                      ! Use this to conserve total energy if consv_te=0
   real    :: consv_te = 0.
   real    :: tau = 0.                ! Time scale (days) for Rayleigh friction
   real    :: rf_cutoff = 30.e2       ! cutoff pressure level for RF
   logical :: filter_phys = .false.
   logical :: dwind_2d = .false.
   logical :: breed_vortex_inline = .false.
   logical :: range_warn = .false.
   logical :: fill = .false.
   logical :: fill_dp = .false.
   logical :: fill_wz = .false.
   logical :: check_negative = .false.
   logical :: non_ortho = .true.
   logical :: moist_phys = .true.     ! Run with moist physics
   logical :: do_held_suarez = .false.
   logical :: do_reed_physics = .false.
   logical :: reed_cond_only = .false.
   logical :: reproduce_sum = .true.  ! Make global sum for consv_te reproduce
   logical :: adjust_dry_mass = .false.
   logical :: fv_debug  = .false.
   logical :: srf_init  = .false.
   logical :: mountain  = .true.
   integer :: remap_option  = 0 
   logical :: z_tracer = .false.      ! transport tracers layer by layer with independent
                                      ! time split; use this if tracer number is huge and/or
                                      ! high resolution (nsplt > 1)

   logical :: old_divg_damp = .false. ! parameter to revert damping parameters back to values
                                      ! defined in a previous revision
                                      ! old_values:
                                      !    d2_bg_k1 = 6.           d2_bg_k2 = 4.
                                      !    d2_divg_max_k1 = 0.02   d2_divg_max_k2 = 0.01
                                      !    damp_k_k1 = 0.          damp_k_k2 = 0.
                                      ! current_values:
                                      !    d2_bg_k1 = 4.           d2_bg_k2 = 2.
                                      !    d2_divg_max_k1 = 0.15   d2_divg_max_k2 = 0.08
                                      !    damp_k_k1 = 0.2         damp_k_k2 = 0.12

   logical :: fv_land = .false.       ! To cold starting the model with USGS terrain
!--------------------------------------------------------------------------------------
! The following options are useful for NWP experiments using datasets on the lat-lon grid
!--------------------------------------------------------------------------------------
   logical :: nudge = .false.         ! Perform nudging
   logical :: nudge_ic = .false.      ! Perform nudging on IC
   logical :: ncep_ic = .false.       ! use NCEP ICs 
   logical :: nggps_ic = .false.      ! use NGGPS ICs 
   logical :: ecmwf_ic = .false.      ! use ECMWF ICs 
   logical :: gfs_phil = .false.      ! if .T., compute geopotential inside of GFS physics
   logical :: agrid_vel_rst = .false. ! if .T., include ua/va (agrid winds) in the restarts
   logical :: use_new_ncep = .false.  ! use the NCEP ICs created after 2014/10/22, if want to read CWAT
   logical :: use_ncep_phy = .false.  ! if .T., separate CWAT by weights of liq_wat and liq_ice in FV_IC
   logical :: fv_diag_ic = .false.    ! reconstruct IC from fv_diagnostics on lat-lon grid
   logical :: external_ic = .false.   ! use ICs from external sources; e.g. lat-lon FV core
                                      ! or NCEP re-analysis; both vertical remapping & horizontal
                                      ! (lat-lon to cubed sphere) interpolation will be done
! Default restart files from the "Memphis" latlon FV core:
   character(len=128) :: res_latlon_dynamics = 'INPUT/fv_rst.res.nc'
   character(len=128) :: res_latlon_tracers  = 'INPUT/atmos_tracers.res.nc'
! The user also needs to copy the "cold start" cubed sphere restart files (fv_core.res.tile1-6)
! to the INPUT dir during runtime
!------------------------------------------------
! Parameters related to non-hydrostatic dynamics:
!------------------------------------------------
   logical :: hydrostatic = .true.
   logical :: phys_hydrostatic = .true.    ! heating/cooling term from the physics is hydrostatic
   logical :: use_hydro_pressure = .false. !  GFS control
   logical :: do_uni_zfull = .false.       ! compute zfull as a simply average of two zhalf
   logical :: hybrid_z    = .false.        ! use hybrid_z for remapping
   logical :: make_nh     = .false.        ! Initialize (w, delz) from hydro restart file 
   logical :: make_hybrid_z  = .false.     ! transform hydrostatic eta-coord IC into non-hydrostatic hybrid_z
   logical :: nudge_qv  = .false.          ! Nudge the water vapor (during na_init) above 30 mb towards HALOE climatology
   real    :: add_noise = -1.              !Amplitude of random noise added upon model startup; <=0 means no noise added

   integer :: a2b_ord = 4    ! order for interpolation from A to B Grid (corners)
   integer :: c2l_ord = 4    ! order for interpolation from D to lat-lon A winds for phys & output

  real(kind=R_GRID) :: dx_const = 1000.    ! spatial resolution for double periodic boundary configuration [m]
  real(kind=R_GRID) :: dy_const = 1000.
  real(kind=R_GRID) :: deglat=15.
  !The following deglat_*, deglon_* options are not used.
  real(kind=R_GRID) :: deglon_start = -30., deglon_stop = 30., &  ! boundaries of latlon patch
                       deglat_start = -30., deglat_stop = 30.

  !Convenience pointers
  integer, pointer :: grid_number

  !f1p
  logical  :: adj_mass_vmr = .false. !TER: This is to reproduce answers for verona patch.  This default can be changed
                                     !     to .true. in the next city release if desired
  
  !integer, pointer :: test_case
  !real,    pointer :: alpha

  end type fv_flags_type

  type fv_nest_bc_type_3d

     !!! CLEANUP: could we have pointers to np[xyz], nest_domain, and the index/weight arrays?

     real, allocatable, dimension(:,:,:) :: west_t1, east_t1, south_t1, north_t1
     real, allocatable, dimension(:,:,:) :: west_t0, east_t0, south_t0, north_t0

     integer :: istag, jstag

     logical :: allocated = .false.
     logical :: initialized = .false.

  end type fv_nest_bc_type_3d

  type fv_nest_bc_type_4d

     real, allocatable, dimension(:,:,:,:) :: west_t1, east_t1, south_t1, north_t1
     real, allocatable, dimension(:,:,:,:) :: west_t0, east_t0, south_t0, north_t0

     integer :: istag, jstag

     logical :: allocated = .false.
     logical :: initialized = .false.

  end type fv_nest_bc_type_4d

  type fv_nest_type

!nested grid flags:

     integer :: refinement = 3  !Refinement wrt parent

     integer :: parent_tile = 1     !Tile (of cubed sphere) in which nested grid lies 
     logical :: nested = .false.
     integer :: nestbctype = 1
     integer :: nsponge = 0
     integer :: nestupdate = 0       
     logical :: twowaynest = .false. 
     integer :: ioffset, joffset !Position of nest within parent grid

     integer :: nest_timestep = 0 !Counter for nested-grid timesteps
     integer :: tracer_nest_timestep = 0 !Counter for nested-grid timesteps
     real    :: s_weight = 1.e-6 !sponge weight
     logical :: first_step = .true.
     integer :: refinement_of_global = 1
     integer :: npx_global
     integer :: upoff = 1 ! currently the same for all variables
     integer :: isu = -999, ieu = -1000, jsu = -999, jeu = -1000 ! limits of update regions on coarse grid 

     type(nest_domain_type) :: nest_domain !Structure holding link from this grid to its parent
     type(nest_domain_type), allocatable :: nest_domain_all(:)

     !Interpolation arrays for grid nesting
     integer, allocatable, dimension(:,:,:) :: ind_h, ind_u, ind_v, ind_b
     real, allocatable, dimension(:,:,:) :: wt_h, wt_u, wt_v, wt_b
     integer, allocatable, dimension(:,:,:) :: ind_update_h

     !These arrays are not allocated by allocate_fv_atmos_type; but instead
     !allocated for all grids, regardless of whether the grid is
     !on a PE of a concurrent run.
     logical, allocatable, dimension(:) :: child_grids

     logical :: parent_proc, child_proc
     logical :: parent_of_twoway = .false.
   
     !These are for time-extrapolated BCs
     type(fv_nest_bc_type_3d) :: delp_bc, u_bc, v_bc, uc_bc, vc_bc, divg_bc
     type(fv_nest_bc_type_3d), allocatable, dimension(:) :: q_bc
#ifndef SW_DYNAMICS
     type(fv_nest_bc_type_3d) :: pt_bc, w_bc, delz_bc
#ifdef USE_COND
     type(fv_nest_bc_type_3d) :: q_con_bc
#ifdef MOIST_CAPPA
     type(fv_nest_bc_type_3d) :: cappa_bc
#endif
#endif
#endif

     !These are for tracer flux BCs
     logical :: do_flux_bcs, do_2way_flux_bcs !For a parent grid; determine whether there is a need to send BCs
     type(restart_file_type) :: bcfile_ne, bcfile_sw

  end type fv_nest_type

  interface allocate_fv_nest_bc_type
     module procedure allocate_fv_nest_bc_type_3d
     module procedure allocate_fv_nest_bc_type_3d_atm
  end interface

  interface deallocate_fv_nest_bc_type
     module procedure deallocate_fv_nest_bc_type_3d
  end interface

  type fv_grid_bounds_type

     integer :: is,  ie,  js,  je
     integer :: isd, ied, jsd, jed
     integer :: isc, iec, jsc, jec

     integer :: ng

  end type fv_grid_bounds_type

  type fv_atmos_type

     logical :: allocated = .false.
     logical :: dummy = .false. ! same as grids_on_this_pe(n)
     integer :: grid_number = 1

     !Timestep-related variables.

     type(time_type) :: time_init, time, run_length, time_end, time_step_atmos

#ifdef GFS_PHYS
     !--- used for GFS PHYSICS only
     real, dimension(2048) :: fdiag = 0.
#endif

     logical :: grid_active = .true. !Always active for now

     !This is kept here instead of in neststruct% simply for convenience
     type(fv_atmos_type), pointer :: parent_grid _null

!-----------------------------------------------------------------------
! Five prognostic state variables for the f-v dynamics
!-----------------------------------------------------------------------
! dyn_state:
! D-grid prognostatic variables: u, v, and delp (and other scalars)
!
!     o--------u(i,j+1)----------o
!     |           |              |
!     |           |              |
!  v(i,j)------scalar(i,j)----v(i+1,j)
!     |           |              |
!     |           |              |
!     o--------u(i,j)------------o
!
! The C grid component is "diagnostic" in that it is predicted every time step
! from the D grid variables.
    real, _allocatable :: u(:,:,:)    _null  ! D grid zonal wind (m/s)
    real, _allocatable :: v(:,:,:)    _null  ! D grid meridional wind (m/s)
    real, _allocatable :: pt(:,:,:)   _null  ! temperature (K)
    real, _allocatable :: delp(:,:,:) _null  ! pressure thickness (pascal)
    real, _allocatable :: q(:,:,:,:)  _null  ! specific humidity and prognostic constituents
    real, _allocatable :: qdiag(:,:,:,:)  _null  ! diagnostic tracers

!----------------------
! non-hydrostatic state:
!----------------------------------------------------------------------
    real, _allocatable ::     w(:,:,:)  _null  ! cell center vertical wind (m/s)
    real, _allocatable ::  delz(:,:,:)  _null  ! layer thickness (meters)
    real, _allocatable ::   ze0(:,:,:)  _null  ! height at layer edges for remapping
    real, _allocatable ::  q_con(:,:,:) _null  ! total condensates

!-----------------------------------------------------------------------
! Auxilliary pressure arrays:
! The 5 vars below can be re-computed from delp and ptop.
!-----------------------------------------------------------------------
! dyn_aux:
    real, _allocatable :: ps (:,:)      _null  ! Surface pressure (pascal)
    real, _allocatable :: pe (:,:,: )   _null  ! edge pressure (pascal)
    real, _allocatable :: pk  (:,:,:)   _null  ! pe**cappa
    real, _allocatable :: peln(:,:,:)   _null  ! ln(pe)
    real, _allocatable :: pkz (:,:,:)   _null  ! finite-volume mean pk

! For phys coupling:
    real, _allocatable :: u_srf(:,:)    _null  ! Surface u-wind
    real, _allocatable :: v_srf(:,:)    _null  ! Surface v-wind
    real, _allocatable :: sgh(:,:)      _null  ! Terrain standard deviation
    real, _allocatable :: oro(:,:)      _null  ! land fraction (1: all land; 0: all water)
    real, _allocatable :: ts(:,:)       _null  ! skin temperature (sst) from NCEP/GFS (K) -- tile
 
!-----------------------------------------------------------------------
! Others:
!-----------------------------------------------------------------------
    real, _allocatable :: phis(:,:)     _null  ! Surface geopotential (g*Z_surf)
    real, _allocatable :: omga(:,:,:)   _null  ! Vertical pressure velocity (pa/s)
    real, _allocatable :: ua(:,:,:)     _null  ! (ua, va) are mostly used as the A grid winds
    real, _allocatable :: va(:,:,:)     _null
    real, _allocatable :: uc(:,:,:)     _null  ! (uc, vc) are mostly used as the C grid winds
    real, _allocatable :: vc(:,:,:)     _null

    real, _allocatable :: ak(:)  _null
    real, _allocatable :: bk(:)  _null

   integer :: ks

! Accumulated Mass flux arrays
    real, _allocatable ::  mfx(:,:,:)  _null
    real, _allocatable ::  mfy(:,:,:)  _null
! Accumulated Courant number arrays
    real, _allocatable ::  cx(:,:,:)  _null
    real, _allocatable ::  cy(:,:,:)  _null

    type(fv_flags_type) :: flagstruct
    
    !! Convenience pointers
    integer, pointer :: npx, npy, npz, ncnst, ng

     integer, allocatable, dimension(:) :: pelist

     type(fv_grid_bounds_type) :: bd

     type(domain2d) :: domain
#if defined(SPMD)

     type(domain2d) :: domain_for_coupler ! domain used in coupled model with halo = 1.

     integer :: num_contact, npes_per_tile, tile, npes_this_grid
     integer :: layout(2), io_layout(2) = (/ 1,1 /)

#endif
     !These do not actually belong to the grid, but to the process
     !integer :: masterproc
     !integer :: gid 

!!!!!!!!!!!!!!!!
! From fv_grid_tools
!!!!!!!!!!!!!!!!


     real    :: ptop

  type(fv_grid_type) :: gridstruct
  

!!!!!!!!!!!!!!!!
!fv_diagnostics!
!!!!!!!!!!!!!!!!

     type(fv_diag_type) :: idiag

!!!!!!!!!!!!!!
! From fv_io !
!!!!!!!!!!!!!!
     type(restart_file_type) :: fv_restart, sst_restart, fv_tile_restart, &
          rsf_restart, mg_restart, lnd_restart, tra_restart

     type(fv_nest_type) :: neststruct

     !Hold on to coarse-grid global grid, so we don't have to waste processor time getting it again when starting to do grid nesting
     real(kind=R_GRID), allocatable, dimension(:,:,:,:) :: grid_global

  integer :: atmos_axes(4)


  end type fv_atmos_type

contains

  subroutine allocate_fv_atmos_type(Atm, isd_in, ied_in, jsd_in, jed_in, is_in, ie_in, js_in, je_in, &
       npx_in, npy_in, npz_in, ndims_in, ncnst_in, nq_in, ng_in, dummy, alloc_2d, ngrids_in)

    !WARNING: Before calling this routine, be sure to have set up the
    ! proper domain parameters from the namelists (as is done in
    ! fv_control_nlm.F90)

    implicit none
    type(fv_atmos_type), intent(INOUT), target :: Atm
    integer, intent(IN) :: isd_in, ied_in, jsd_in, jed_in, is_in, ie_in, js_in, je_in
    integer, intent(IN) :: npx_in, npy_in, npz_in, ndims_in, ncnst_in, nq_in, ng_in
    logical, intent(IN) :: dummy, alloc_2d
    integer, intent(IN) :: ngrids_in
    integer:: isd, ied, jsd, jed, is, ie, js, je
    integer:: npx, npy, npz, ndims, ncnst, nq, ng

    !For 2D utility arrays
    integer:: isd_2d, ied_2d, jsd_2d, jed_2d, is_2d, ie_2d, js_2d, je_2d
    integer:: npx_2d, npy_2d, npz_2d, ndims_2d, ncnst_2d, nq_2d, ng_2d

    integer :: i,j,k, ns, n

    if (atm%allocated) return

    if (dummy) then
       isd     =  0   
       ied=   -1   
       jsd=   0   
       jed=   -1   
       is=   0    
       ie=   -1    
       js=   0    
       je=   -1    
       npx=   1   
       npy=   1   
       npz=   1   
       ndims=   1 
       ncnst=   1 
       nq=   1
       ng     =   1   
    else
       isd     =  isd_in   
       ied=   ied_in   
       jsd=   jsd_in   
       jed=   jed_in   
       is=   is_in    
       ie=   ie_in    
       js=   js_in    
       je=   je_in    
       npx=   npx_in   
       npy=   npy_in   
       npz=   npz_in   
       ndims=   ndims_in 
       ncnst=   ncnst_in 
       nq=   nq_in
       ng     =   ng_in    
    endif

    if ((.not. dummy) .or. alloc_2d) then
       isd_2d     =  isd_in   
       ied_2d=   ied_in   
       jsd_2d=   jsd_in   
       jed_2d=   jed_in   
       is_2d=   is_in    
       ie_2d=   ie_in    
       js_2d=   js_in    
       je_2d=   je_in    
       npx_2d=   npx_in   
       npy_2d=   npy_in   
       npz_2d=   npz_in   
       ndims_2d=   ndims_in 
       ncnst_2d=   ncnst_in 
       nq_2d=   nq_in 
       ng_2d     =   ng_in 
    else
       isd_2d     =  0   
       ied_2d=   -1   
       jsd_2d=   0   
       jed_2d=   -1   
       is_2d=   0    
       ie_2d=   -1    
       js_2d=   0    
       je_2d=   -1    
       npx_2d=   1   
       npy_2d=   1   
       npz_2d=   0 !for ak, bk   
       ndims_2d=   1 
       ncnst_2d=   1 
       nq_2d=   1 
       ng_2d     =   1        
    endif

!This should be set up in fv_mp_nlm_mod
!!$    Atm%bd%isd = isd_in
!!$    Atm%bd%ied = ied_in
!!$    Atm%bd%jsd = jsd_in
!!$    Atm%bd%jed = jed_in
!!$
!!$    Atm%bd%is = is_in
!!$    Atm%bd%ie = ie_in
!!$    Atm%bd%js = js_in
!!$    Atm%bd%je = je_in
!!$
!!$    Atm%bd%isc = Atm%bd%is
!!$    Atm%bd%iec = Atm%bd%ie
!!$    Atm%bd%jsc = Atm%bd%js
!!$    Atm%bd%jec = Atm%bd%je

    atm%bd%ng  = ng

    !Convenience pointers
    atm%npx => atm%flagstruct%npx
    atm%npy => atm%flagstruct%npy
    atm%npz => atm%flagstruct%npz
    atm%ncnst => atm%flagstruct%ncnst

    atm%ng => atm%bd%ng

!!$    Atm%npx = npx_in
!!$    Atm%npy = npy_in
!!$    Atm%npz = npz_in
    atm%flagstruct%ndims = ndims_in

    allocate (    atm%u(isd:ied  ,jsd:jed+1,npz) )
    allocate (    atm%v(isd:ied+1,jsd:jed  ,npz) )

    allocate (   atm%pt(isd:ied  ,jsd:jed  ,npz) )
    allocate ( atm%delp(isd:ied  ,jsd:jed  ,npz) )
    allocate (    atm%q(isd:ied  ,jsd:jed  ,npz, nq) )
    allocate (atm%qdiag(isd:ied  ,jsd:jed  ,npz, nq+1:ncnst) )

    ! Allocate Auxilliary pressure arrays
    allocate (   atm%ps(isd:ied  ,jsd:jed) )
    allocate (   atm%pe(is-1:ie+1, npz+1,js-1:je+1) )
    allocate (   atm%pk(is:ie    ,js:je  , npz+1) )
    allocate ( atm%peln(is:ie,npz+1,js:je) )
    allocate (  atm%pkz(is:ie,js:je,npz) )

    allocate ( atm%u_srf(is:ie,js:je) )
    allocate ( atm%v_srf(is:ie,js:je) )

    if ( atm%flagstruct%fv_land ) then
       allocate ( atm%sgh(is:ie,js:je) )
       allocate ( atm%oro(is:ie,js:je) )
    else
       allocate ( atm%oro(1,1) )
    endif

    ! Allocate others
    allocate ( atm%ts(is:ie,js:je) )
    allocate ( atm%phis(isd:ied  ,jsd:jed  ) )
    allocate ( atm%omga(isd:ied  ,jsd:jed  ,npz) ); atm%omga=0.
    allocate (   atm%ua(isd:ied  ,jsd:jed  ,npz) )
    allocate (   atm%va(isd:ied  ,jsd:jed  ,npz) )
    allocate (   atm%uc(isd:ied+1,jsd:jed  ,npz) )
    allocate (   atm%vc(isd:ied  ,jsd:jed+1,npz) )
    ! For tracer transport:
    allocate ( atm%mfx(is:ie+1, js:je,  npz) )
    allocate ( atm%mfy(is:ie  , js:je+1,npz) )
    allocate (  atm%cx(is:ie+1, jsd:jed, npz) )
    allocate (  atm%cy(isd:ied ,js:je+1, npz) )

    allocate (  atm%ak(npz_2d+1) )
    allocate (  atm%bk(npz_2d+1) )

    !--------------------------
    ! Non-hydrostatic dynamics:
    !--------------------------
    if ( atm%flagstruct%hydrostatic ) then
       !Note length-one initialization if hydrostatic = .true.
       allocate (    atm%w(isd:isd, jsd:jsd  ,1) )
       allocate ( atm%delz(isd:isd, jsd:jsd  ,1) )
       allocate (  atm%ze0(is:is, js:js  ,1) )
    else
       allocate (    atm%w(isd:ied, jsd:jed  ,npz  ) )
       allocate ( atm%delz(isd:ied, jsd:jed  ,npz) )
       if( atm%flagstruct%hybrid_z ) then
          allocate (  atm%ze0(is:ie, js:je ,npz+1) )
       else
          allocate (  atm%ze0(is:is, js:js  ,1) )
       endif
       !         allocate ( mono(isd:ied, jsd:jed, npz))
    endif

#ifdef USE_COND
      allocate ( atm%q_con(isd:ied,jsd:jed,1:npz) )
#else
      allocate ( atm%q_con(isd:isd,jsd:jsd,1) )
#endif

#ifndef NO_TOUCH_MEM
! Notes by SJL
! Place the memory in the optimal shared mem space
! This will help the scaling with OpenMP
!$OMP parallel do default(none) shared(isd,ied,jsd,jed,npz,Atm,nq,ncnst)
     do k=1, npz
        do j=jsd, jed
           do i=isd, ied
                atm%ua(i,j,k) = real_big
                atm%va(i,j,k) = real_big
                atm%pt(i,j,k) = real_big
              atm%delp(i,j,k) = real_big
           enddo
        enddo
        do j=jsd, jed+1
           do i=isd, ied
               atm%u(i,j,k) = real_big
              atm%vc(i,j,k) = real_big
           enddo
        enddo
        do j=jsd, jed
           do i=isd, ied+1
               atm%v(i,j,k) = real_big
              atm%uc(i,j,k) = real_big
           enddo
        enddo
        if ( .not. atm%flagstruct%hydrostatic ) then
           do j=jsd, jed
              do i=isd, ied
                    atm%w(i,j,k) = real_big
                 atm%delz(i,j,k) = real_big
              enddo
           enddo
        endif
        do n=1,nq
        do j=jsd, jed
           do i=isd, ied
              atm%q(i,j,k,n) = real_big
           enddo
        enddo
        enddo
        do n=nq+1,ncnst
        do j=jsd, jed
           do i=isd, ied
              atm%qdiag(i,j,k,n) = real_big
           enddo
        enddo
        enddo
     enddo
#endif

    allocate ( atm%gridstruct% area(isd_2d:ied_2d  ,jsd_2d:jed_2d  ) )   ! Cell Centered
    allocate ( atm%gridstruct% area_64(isd_2d:ied_2d  ,jsd_2d:jed_2d  ) ) ! Cell Centered
    allocate ( atm%gridstruct%rarea(isd_2d:ied_2d  ,jsd_2d:jed_2d  ) )   ! Cell Centered
    
    allocate ( atm%gridstruct% area_c(isd_2d:ied_2d+1,jsd_2d:jed_2d+1) )   ! Cell Corners
    allocate ( atm%gridstruct% area_c_64(isd_2d:ied_2d+1,jsd_2d:jed_2d+1) )! Cell Corners
    allocate ( atm%gridstruct%rarea_c(isd_2d:ied_2d+1,jsd_2d:jed_2d+1) )   ! Cell Corners
   
    atm%gridstruct% area = 0.0
    atm%gridstruct% area_64 = 0.0
    atm%gridstruct%rarea = 0.0

    atm%gridstruct% area_c = 0.0
    atm%gridstruct% area_c_64 = 0.0
    atm%gridstruct%rarea_c = 0.0
 
    allocate ( atm%gridstruct% dx(isd_2d:ied_2d  ,jsd_2d:jed_2d+1) )
    allocate ( atm%gridstruct% dx_64(isd_2d:ied_2d  ,jsd_2d:jed_2d+1) )
    allocate ( atm%gridstruct%rdx(isd_2d:ied_2d  ,jsd_2d:jed_2d+1) )
    allocate ( atm%gridstruct% dy(isd_2d:ied_2d+1,jsd_2d:jed_2d  ) )
    allocate ( atm%gridstruct% dy_64(isd_2d:ied_2d+1,jsd_2d:jed_2d  ) )
    allocate ( atm%gridstruct%rdy(isd_2d:ied_2d+1,jsd_2d:jed_2d  ) )
    
    allocate ( atm%gridstruct% dxc(isd_2d:ied_2d+1,jsd_2d:jed_2d  ) )
    allocate ( atm%gridstruct% dxc_64(isd_2d:ied_2d+1,jsd_2d:jed_2d  ) )
    allocate ( atm%gridstruct%rdxc(isd_2d:ied_2d+1,jsd_2d:jed_2d  ) )
    allocate ( atm%gridstruct% dyc(isd_2d:ied_2d  ,jsd_2d:jed_2d+1) )
    allocate ( atm%gridstruct% dyc_64(isd_2d:ied_2d  ,jsd_2d:jed_2d+1) )
    allocate ( atm%gridstruct%rdyc(isd_2d:ied_2d  ,jsd_2d:jed_2d+1) )
    
    allocate ( atm%gridstruct% dxa(isd_2d:ied_2d  ,jsd_2d:jed_2d  ) )
    allocate ( atm%gridstruct% dxa_64(isd_2d:ied_2d  ,jsd_2d:jed_2d  ) )
    allocate ( atm%gridstruct%rdxa(isd_2d:ied_2d  ,jsd_2d:jed_2d  ) )
    allocate ( atm%gridstruct% dya(isd_2d:ied_2d  ,jsd_2d:jed_2d  ) )
    allocate ( atm%gridstruct% dya_64(isd_2d:ied_2d  ,jsd_2d:jed_2d  ) )
    allocate ( atm%gridstruct%rdya(isd_2d:ied_2d  ,jsd_2d:jed_2d  ) )
    
    allocate ( atm%gridstruct%grid (isd_2d:ied_2d+1,jsd_2d:jed_2d+1,1:ndims_2d) )
    allocate ( atm%gridstruct%grid_64 (isd_2d:ied_2d+1,jsd_2d:jed_2d+1,1:ndims_2d) )
    allocate ( atm%gridstruct%agrid(isd_2d:ied_2d  ,jsd_2d:jed_2d  ,1:ndims_2d) )
    allocate ( atm%gridstruct%agrid_64(isd_2d:ied_2d  ,jsd_2d:jed_2d  ,1:ndims_2d) )
    allocate ( atm%gridstruct% sina(isd_2d:ied_2d+1,jsd_2d:jed_2d+1) )   ! SIN(angle of intersection)
    allocate ( atm%gridstruct% sina_64(isd_2d:ied_2d+1,jsd_2d:jed_2d+1) )   ! SIN(angle of intersection)
    allocate ( atm%gridstruct%rsina(is_2d:ie_2d+1,js_2d:je_2d+1) )      ! Why is the size different?
    allocate ( atm%gridstruct% cosa(isd_2d:ied_2d+1,jsd_2d:jed_2d+1) )   ! COS(angle of intersection)
    allocate ( atm%gridstruct% cosa_64(isd_2d:ied_2d+1,jsd_2d:jed_2d+1) )   ! COS(angle of intersection)
    
    allocate ( atm%gridstruct%  e1(3,isd_2d:ied_2d+1,jsd_2d:jed_2d+1) )
    allocate ( atm%gridstruct%  e2(3,isd_2d:ied_2d+1,jsd_2d:jed_2d+1) )

    allocate (atm%gridstruct%iinta(4, isd_2d:ied_2d ,jsd_2d:jed_2d), &
         atm%gridstruct%jinta(4, isd_2d:ied_2d ,jsd_2d:jed_2d),  &
         atm%gridstruct%iintb(4, is_2d:ie_2d+1 ,js_2d:je_2d+1), &
         atm%gridstruct%jintb(4, is_2d:ie_2d+1 ,js_2d:je_2d+1) )

    allocate ( atm%gridstruct%edge_s(npx_2d) )
    allocate ( atm%gridstruct%edge_n(npx_2d) )
    allocate ( atm%gridstruct%edge_w(npy_2d) )
    allocate ( atm%gridstruct%edge_e(npy_2d) )

    allocate ( atm%gridstruct%edge_vect_s(isd_2d:ied_2d) )
    allocate ( atm%gridstruct%edge_vect_n(isd_2d:ied_2d) )
    allocate ( atm%gridstruct%edge_vect_w(jsd_2d:jed_2d) )
    allocate ( atm%gridstruct%edge_vect_e(jsd_2d:jed_2d) )

    allocate ( atm%gridstruct%ex_s(npx_2d) )
    allocate ( atm%gridstruct%ex_n(npx_2d) )
    allocate ( atm%gridstruct%ex_w(npy_2d) )
    allocate ( atm%gridstruct%ex_e(npy_2d) )


    allocate (  atm%gridstruct%l2c_u(is_2d:ie_2d,  js_2d:je_2d+1) )
    allocate (  atm%gridstruct%l2c_v(is_2d:ie_2d+1,js_2d:je_2d) )

    ! For diveregnce damping:
    allocate (  atm%gridstruct%divg_u(isd_2d:ied_2d,  jsd_2d:jed_2d+1) )
    allocate (  atm%gridstruct%divg_v(isd_2d:ied_2d+1,jsd_2d:jed_2d) )
    ! For del6 diffusion:
    allocate (  atm%gridstruct%del6_u(isd_2d:ied_2d,  jsd_2d:jed_2d+1) )
    allocate (  atm%gridstruct%del6_v(isd_2d:ied_2d+1,jsd_2d:jed_2d) )

    allocate (  atm%gridstruct%z11(is_2d-1:ie_2d+1,js_2d-1:je_2d+1) )
    allocate (  atm%gridstruct%z12(is_2d-1:ie_2d+1,js_2d-1:je_2d+1) )
    allocate (  atm%gridstruct%z21(is_2d-1:ie_2d+1,js_2d-1:je_2d+1) )
    allocate (  atm%gridstruct%z22(is_2d-1:ie_2d+1,js_2d-1:je_2d+1) )

!   if (.not.Atm%flagstruct%hydrostatic)    &
!   allocate (  Atm%gridstruct%w00(is_2d-1:ie_2d+1,js_2d-1:je_2d+1) )

    allocate (  atm%gridstruct%a11(is_2d-1:ie_2d+1,js_2d-1:je_2d+1) )
    allocate (  atm%gridstruct%a12(is_2d-1:ie_2d+1,js_2d-1:je_2d+1) )
    allocate (  atm%gridstruct%a21(is_2d-1:ie_2d+1,js_2d-1:je_2d+1) )
    allocate (  atm%gridstruct%a22(is_2d-1:ie_2d+1,js_2d-1:je_2d+1) )
    allocate ( atm%gridstruct%vlon(is_2d-2:ie_2d+2,js_2d-2:je_2d+2,3) )
    allocate ( atm%gridstruct%vlat(is_2d-2:ie_2d+2,js_2d-2:je_2d+2,3) )
    ! Coriolis parameters:
    allocate ( atm%gridstruct%f0(isd_2d:ied_2d  ,jsd_2d:jed_2d  ) )
    allocate ( atm%gridstruct%fC(isd_2d:ied_2d+1,jsd_2d:jed_2d+1) )

    ! Corner unit vectors:
    allocate( atm%gridstruct%ee1(3,isd_2d:ied_2d+1,jsd_2d:jed_2d+1) )
    allocate( atm%gridstruct%ee2(3,isd_2d:ied_2d+1,jsd_2d:jed_2d+1) )

    ! Center unit vectors:
    allocate( atm%gridstruct%ec1(3,isd_2d:ied_2d,jsd_2d:jed_2d) )
    allocate( atm%gridstruct%ec2(3,isd_2d:ied_2d,jsd_2d:jed_2d) )

    ! Edge unit vectors:
    allocate( atm%gridstruct%ew(3,isd_2d:ied_2d+1,jsd_2d:jed_2d,  2) )
    allocate( atm%gridstruct%es(3,isd_2d:ied_2d  ,jsd_2d:jed_2d+1,2) )

    ! Edge unit "Normal" vectors: (for omega computation)
    allocate( atm%gridstruct%en1(3,is_2d:ie_2d,  js_2d:je_2d+1) )   ! E-W edges
    allocate( atm%gridstruct%en2(3,is_2d:ie_2d+1,js_2d:je_2d  ) )   ! N-S egdes

    allocate ( atm%gridstruct%cosa_u(isd_2d:ied_2d+1,jsd_2d:jed_2d) )
    allocate ( atm%gridstruct%sina_u(isd_2d:ied_2d+1,jsd_2d:jed_2d) )
    allocate ( atm%gridstruct%rsin_u(isd_2d:ied_2d+1,jsd_2d:jed_2d) )

    allocate ( atm%gridstruct%cosa_v(isd_2d:ied_2d,jsd_2d:jed_2d+1) )
    allocate ( atm%gridstruct%sina_v(isd_2d:ied_2d,jsd_2d:jed_2d+1) )
    allocate ( atm%gridstruct%rsin_v(isd_2d:ied_2d,jsd_2d:jed_2d+1) )

    allocate ( atm%gridstruct%cosa_s(isd_2d:ied_2d,jsd_2d:jed_2d) )    ! cell center

    allocate (  atm%gridstruct%rsin2(isd_2d:ied_2d,jsd_2d:jed_2d) )    ! cell center


    ! Super (composite) grid:

    !     9---4---8
    !     |       |
    !     1   5   3
    !     |       |
    !     6---2---7

    allocate ( atm%gridstruct%cos_sg(isd_2d:ied_2d,jsd_2d:jed_2d,9) )
    allocate ( atm%gridstruct%sin_sg(isd_2d:ied_2d,jsd_2d:jed_2d,9) )

    allocate( atm%gridstruct%eww(3,4) )
    allocate( atm%gridstruct%ess(3,4) )

    if (atm%neststruct%nested) then

       allocate(atm%neststruct%ind_h(isd:ied,jsd:jed,4))
       allocate(atm%neststruct%ind_u(isd:ied,jsd:jed+1,4))
       allocate(atm%neststruct%ind_v(isd:ied+1,jsd:jed,4))

       allocate(atm%neststruct%wt_h(isd:ied,   jsd:jed,  4))
       allocate(atm%neststruct%wt_u(isd:ied,   jsd:jed+1,4))
       allocate(atm%neststruct%wt_v(isd:ied+1, jsd:jed,  4))
       allocate(atm%neststruct%ind_b(isd:ied+1,jsd:jed+1,4))
       allocate(atm%neststruct%wt_b(isd:ied+1, jsd:jed+1,4))

       ns = atm%neststruct%nsponge

       call allocate_fv_nest_bc_type(atm%neststruct%delp_BC,atm,ns,0,0,dummy)
       call allocate_fv_nest_bc_type(atm%neststruct%u_BC,atm,ns,0,1,dummy)
       call allocate_fv_nest_bc_type(atm%neststruct%v_BC,atm,ns,1,0,dummy)
       call allocate_fv_nest_bc_type(atm%neststruct%uc_BC,atm,ns,1,0,dummy)
       call allocate_fv_nest_bc_type(atm%neststruct%vc_BC,atm,ns,0,1,dummy)
       call allocate_fv_nest_bc_type(atm%neststruct%divg_BC,atm,ns,1,1,dummy)

       if (ncnst > 0) then
          allocate(atm%neststruct%q_BC(ncnst))
          do n=1,ncnst
             call allocate_fv_nest_bc_type(atm%neststruct%q_BC(n),atm,ns,0,0,dummy)
          enddo
       endif

#ifndef SW_DYNAMICS

       call allocate_fv_nest_bc_type(atm%neststruct%pt_BC,atm,ns,0,0,dummy)
#ifdef USE_COND
       call allocate_fv_nest_bc_type(atm%neststruct%q_con_BC,atm,ns,0,0,dummy)
#ifdef MOIST_CAPPA
       call allocate_fv_nest_bc_type(atm%neststruct%cappa_BC,atm,ns,0,0,dummy)
#endif
#endif
       if (.not.atm%flagstruct%hydrostatic) then
          call allocate_fv_nest_bc_type(atm%neststruct%w_BC,atm,ns,0,0,dummy)
          call allocate_fv_nest_bc_type(atm%neststruct%delz_BC,atm,ns,0,0,dummy)
       endif

#endif

       if (atm%neststruct%twowaynest) allocate(&
            atm%neststruct%ind_update_h( &
              atm%parent_grid%bd%isd:atm%parent_grid%bd%ied+1, &
              atm%parent_grid%bd%jsd:atm%parent_grid%bd%jed+1,2))

    end if

    !--- Do the memory allocation only for nested model
    if( ngrids_in > 1 ) then
       if (atm%flagstruct%grid_type < 4) then
          if (atm%neststruct%nested) then
             allocate(atm%grid_global(1-ng_2d:npx_2d  +ng_2d,1-ng_2d:npy_2d  +ng_2d,2,1))
          else
             allocate(atm%grid_global(1-ng_2d:npx_2d  +ng_2d,1-ng_2d:npy_2d  +ng_2d,2,1:6))
          endif
       end if
    endif

    atm%allocated = .true.
    if (dummy) atm%dummy = .true.
    
  end subroutine allocate_fv_atmos_type

  subroutine deallocate_fv_atmos_type(Atm)

    implicit none
    type(fv_atmos_type), intent(INOUT) :: Atm

    integer :: n

    if (.not.atm%allocated) return

    deallocate (    atm%u )
    deallocate (    atm%v )
    deallocate (   atm%pt )
    deallocate ( atm%delp )
    deallocate (    atm%q )
    deallocate (    atm%qdiag )
    deallocate (   atm%ps )
    deallocate (   atm%pe )
    deallocate (   atm%pk )
    deallocate ( atm%peln )
    deallocate (  atm%pkz )
    deallocate ( atm%phis )
    deallocate ( atm%omga )
    deallocate (   atm%ua )
    deallocate (   atm%va )
    deallocate (   atm%uc )
    deallocate (   atm%vc )
    deallocate ( atm%mfx )
    deallocate ( atm%mfy )
    deallocate (  atm%cx )
    deallocate (  atm%cy )
    deallocate (  atm%ak )
    deallocate (  atm%bk )

    deallocate ( atm%u_srf )
    deallocate ( atm%v_srf )
    if( atm%flagstruct%fv_land ) deallocate ( atm%sgh )
    deallocate ( atm%oro )

    deallocate ( atm%w )
    deallocate ( atm%delz  )
    deallocate ( atm%ze0   )
    deallocate ( atm%q_con )

    deallocate ( atm%gridstruct% area )   ! Cell Centered
    deallocate ( atm%gridstruct%rarea )   ! Cell Centered
    
    deallocate ( atm%gridstruct% area_c )  ! Cell Corners
    deallocate ( atm%gridstruct%rarea_c )  ! Cell Corners
    
    deallocate ( atm%gridstruct% dx )
    deallocate ( atm%gridstruct%rdx )
    deallocate ( atm%gridstruct% dy )
    deallocate ( atm%gridstruct%rdy )
    
    deallocate ( atm%gridstruct% dxc )
    deallocate ( atm%gridstruct%rdxc )
    deallocate ( atm%gridstruct% dyc )
    deallocate ( atm%gridstruct%rdyc )
    
    deallocate ( atm%gridstruct% dxa )
    deallocate ( atm%gridstruct%rdxa )
    deallocate ( atm%gridstruct% dya )
    deallocate ( atm%gridstruct%rdya )
    
    deallocate ( atm%gridstruct%grid  )
    deallocate ( atm%gridstruct%agrid )
    deallocate ( atm%gridstruct%sina )   ! SIN(angle of intersection)
    deallocate ( atm%gridstruct%cosa )   ! COS(angle of intersection)
    
    deallocate ( atm%gridstruct%  e1 )
    deallocate ( atm%gridstruct%  e2 )




    deallocate (atm%gridstruct%iinta, &
         atm%gridstruct%jinta,  &
         atm%gridstruct%iintb, &
         atm%gridstruct%jintb )

    deallocate ( atm%gridstruct%edge_s )
    deallocate ( atm%gridstruct%edge_n )
    deallocate ( atm%gridstruct%edge_w )
    deallocate ( atm%gridstruct%edge_e )

    deallocate ( atm%gridstruct%edge_vect_s )
    deallocate ( atm%gridstruct%edge_vect_n )
    deallocate ( atm%gridstruct%edge_vect_w )
    deallocate ( atm%gridstruct%edge_vect_e )

    deallocate ( atm%gridstruct%ex_s )
    deallocate ( atm%gridstruct%ex_n )
    deallocate ( atm%gridstruct%ex_w )
    deallocate ( atm%gridstruct%ex_e )


    deallocate (  atm%gridstruct%l2c_u )
    deallocate (  atm%gridstruct%l2c_v )
    ! For diveregnce damping:
    deallocate (  atm%gridstruct%divg_u )
    deallocate (  atm%gridstruct%divg_v )
    ! For del6 diffusion:

    deallocate (  atm%gridstruct%z11 )
    deallocate (  atm%gridstruct%z12 )
    deallocate (  atm%gridstruct%z21 )
    deallocate (  atm%gridstruct%z22 )

    deallocate (  atm%gridstruct%a11 )
    deallocate (  atm%gridstruct%a12 )
    deallocate (  atm%gridstruct%a21 )
    deallocate (  atm%gridstruct%a22 )
    deallocate ( atm%gridstruct%vlon )
    deallocate ( atm%gridstruct%vlat )
    ! Coriolis parameters:
    deallocate ( atm%gridstruct%f0 )
    deallocate ( atm%gridstruct%fC )

    ! Corner unit vectors:
    deallocate( atm%gridstruct%ee1 )
    deallocate( atm%gridstruct%ee2 )

    ! Center unit vectors:
    deallocate( atm%gridstruct%ec1 )
    deallocate( atm%gridstruct%ec2 )

    ! Edge unit vectors:
    deallocate( atm%gridstruct%ew )
    deallocate( atm%gridstruct%es )

    ! Edge unit "Normal" vectors: (for omega computation)
    deallocate( atm%gridstruct%en1 )   ! E-W edges
    deallocate( atm%gridstruct%en2 )   ! N-S egdes

    deallocate ( atm%gridstruct%cosa_u )
    deallocate ( atm%gridstruct%sina_u )
    deallocate ( atm%gridstruct%rsin_u )

    deallocate ( atm%gridstruct%cosa_v )
    deallocate ( atm%gridstruct%sina_v )
    deallocate ( atm%gridstruct%rsin_v )

    deallocate ( atm%gridstruct%cosa_s )    ! cell center

    deallocate (  atm%gridstruct%rsin2 )    ! cell center


    ! Super (composite) grid:

    !     9---4---8
    !     |       |
    !     1   5   3
    !     |       |
    !     6---2---7

    deallocate ( atm%gridstruct%cos_sg )
    deallocate ( atm%gridstruct%sin_sg )

    deallocate( atm%gridstruct%eww )
    deallocate( atm%gridstruct%ess )

    if (atm%neststruct%nested) then
       deallocate(atm%neststruct%ind_h)
       deallocate(atm%neststruct%ind_u)
       deallocate(atm%neststruct%ind_v)

       deallocate(atm%neststruct%wt_h)
       deallocate(atm%neststruct%wt_u)
       deallocate(atm%neststruct%wt_v)

       deallocate(atm%neststruct%ind_b)
       deallocate(atm%neststruct%wt_b)

       call deallocate_fv_nest_bc_type(atm%neststruct%delp_BC)
       call deallocate_fv_nest_bc_type(atm%neststruct%u_BC)
       call deallocate_fv_nest_bc_type(atm%neststruct%v_BC)
       call deallocate_fv_nest_bc_type(atm%neststruct%uc_BC)
       call deallocate_fv_nest_bc_type(atm%neststruct%vc_BC)
       call deallocate_fv_nest_bc_type(atm%neststruct%divg_BC)

       if (allocated(atm%neststruct%q_BC)) then
          do n=1,size(atm%neststruct%q_BC)
             call deallocate_fv_nest_bc_type(atm%neststruct%q_BC(n))
          enddo
       endif

#ifndef SW_DYNAMICS
       call deallocate_fv_nest_bc_type(atm%neststruct%pt_BC)
#ifdef USE_COND
       call deallocate_fv_nest_bc_type(atm%neststruct%q_con_BC)
#ifdef MOIST_CAPPA
       call deallocate_fv_nest_bc_type(atm%neststruct%cappa_BC)
#endif
#endif
       if (.not.atm%flagstruct%hydrostatic) then
          call deallocate_fv_nest_bc_type(atm%neststruct%w_BC)
          call deallocate_fv_nest_bc_type(atm%neststruct%delz_BC)
       endif
#endif


       if (atm%neststruct%twowaynest) deallocate(atm%neststruct%ind_update_h)

    end if

    if (atm%flagstruct%grid_type < 4) then
       if(allocated(atm%grid_global)) deallocate(atm%grid_global)
    end if
    
    atm%allocated = .false.

  end subroutine deallocate_fv_atmos_type


subroutine allocate_fv_nest_bc_type_3d_atm(BC,Atm,ns,istag,jstag,dummy)

  type(fv_nest_BC_type_3D), intent(INOUT) :: BC
  type(fv_atmos_type), intent(IN) :: Atm
  integer, intent(IN) :: ns, istag, jstag
  logical, intent(IN) :: dummy

  integer :: is, ie, js, je, isd, ied, jsd, jed, npx, npy, npz, ng

  if (bc%allocated) return

  is = atm%bd%is
  ie = atm%bd%ie
  js = atm%bd%js
  je = atm%bd%je

  isd = atm%bd%isd
  ied = atm%bd%ied
  jsd = atm%bd%jsd
  jed = atm%bd%jed

  npx = atm%npx
  npy = atm%npy
  npz = atm%npz

  ng = atm%ng

  call allocate_fv_nest_bc_type_3d(bc,is,ie,js,je,isd,ied,jsd,jed,npx,npy,npz,ng,ns,istag,jstag,dummy)


end subroutine allocate_fv_nest_bc_type_3d_atm

subroutine allocate_fv_nest_bc_type_3d(BC,is,ie,js,je,isd,ied,jsd,jed,npx,npy,npz,ng,ns,istag,jstag,dummy)

  type(fv_nest_BC_type_3D), intent(INOUT) :: BC
  integer, intent(IN) :: ns, istag, jstag
  logical, intent(IN) :: dummy

  integer, intent(IN) :: is, ie, js, je, isd, ied, jsd, jed, npx, npy, npz, ng

  if (bc%allocated) return


  if (ie == npx-1 .and. .not. dummy) then
     allocate(bc%east_t1(ie+1-ns+istag:ied+istag,jsd:jed+jstag,npz))
     allocate(bc%east_t0(ie+1-ns+istag:ied+istag,jsd:jed+jstag,npz))
     do k=1,npz
     do j=jsd,jed+jstag
     do i=ie+1-ns+istag,ied+istag
        bc%east_t1(i,j,k) = 0.
        bc%east_t0(i,j,k) = 0.
     enddo
     enddo
     enddo
  else
     allocate(bc%east_t1(1,1,npz))
     allocate(bc%east_t0(1,1,npz))
  end if

  if (js == 1 .and. .not. dummy) then
     allocate(bc%south_t1(isd:ied+istag,jsd:js-1+ns,npz))
     allocate(bc%south_t0(isd:ied+istag,jsd:js-1+ns,npz))
     do k=1,npz
     do j=jsd,js-1+ns
     do i=isd,ied+istag
        bc%south_t1(i,j,k) = 0.
        bc%south_t0(i,j,k) = 0.
     enddo
     enddo
     enddo
  else
     allocate(bc%south_t1(1,1,npz))
     allocate(bc%south_t0(1,1,npz))
  end if

  if (is == 1 .and. .not. dummy) then
     allocate(bc%west_t1(isd:is-1+ns,jsd:jed+jstag,npz))
     allocate(bc%west_t0(isd:is-1+ns,jsd:jed+jstag,npz))
     do k=1,npz
     do j=jsd,jed+jstag
     do i=isd,is-1+ns
        bc%west_t1(i,j,k) = 0.
        bc%west_t0(i,j,k) = 0.
     enddo
     enddo
     enddo
  else
     allocate(bc%west_t1(1,1,npz))
     allocate(bc%west_t0(1,1,npz))
  end if

  if (je == npy-1 .and. .not. dummy) then
     allocate(bc%north_t1(isd:ied+istag,je+1-ns+jstag:jed+jstag,npz))
     allocate(bc%north_t0(isd:ied+istag,je+1-ns+jstag:jed+jstag,npz))
     do k=1,npz
     do j=je+1-ns+jstag,jed+jstag
     do i=isd,ied+istag
        bc%north_t1(i,j,k) = 0.
        bc%north_t0(i,j,k) = 0.
     enddo
     enddo
     enddo
  else
     allocate(bc%north_t1(1,1,npz))
     allocate(bc%north_t0(1,1,npz))
  end if

  bc%allocated = .true.

end subroutine allocate_fv_nest_bc_type_3d

subroutine deallocate_fv_nest_bc_type_3d(BC)

  type(fv_nest_BC_type_3d) :: BC

  if (.not. bc%allocated) return

     deallocate(bc%north_t1)
     deallocate(bc%south_t1)
     deallocate(bc%west_t1)
     deallocate(bc%east_t1)

  if (allocated(bc%north_t0)) then
     deallocate(bc%north_t0)
     deallocate(bc%south_t0)
     deallocate(bc%west_t0)
     deallocate(bc%east_t0)
  endif

  bc%allocated = .false.

end subroutine deallocate_fv_nest_bc_type_3d


end module fv_arrays_nlm_mod