fv_arrays.F90

Go to the documentation of this file.

!***********************************************************************
!*                   GNU Lesser General Public License                 
!*
!* This file is part of the FV3 dynamical core.
!*
!* The FV3 dynamical core is free software: you can redistribute it 
!* and/or modify it under the terms of the
!* GNU Lesser General Public License as published by the
!* Free Software Foundation, either version 3 of the License, or 
!* (at your option) any later version.
!*
!* The FV3 dynamical core is distributed in the hope that it will be 
!* useful, but WITHOUT ANYWARRANTY; without even the implied warranty 
!* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
!* See the GNU General Public License for more details.
!*
!* You should have received a copy of the GNU Lesser General Public
!* License along with the FV3 dynamical core.  
!* If not, see <http://www.gnu.org/licenses/>.
!***********************************************************************

module fv_arrays_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_mod,               only: mpp_broadcast
  use platform_mod,          only: r8_kind
  public

  integer, public, parameter :: r_grid = r8_kind

  !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, id_dbztop, id_dbz_m10c, &
           id_ctz, id_w1km, id_wmaxup, id_wmaxdn, id_cape, id_cin, id_diss

! Selected p-level fields from 3D variables:
 integer :: id_vort200, id_vort500, id_w500, id_w700
 integer :: id_vort850, id_w850, id_x850, id_srh25, &
            id_uh03, id_uh25, id_theta_e,  &
            id_w200, id_s200, id_sl12, id_sl13, id_w5km, id_rain5km, id_w2500m
 integer :: id_srh1, id_srh3, id_ustm, id_vstm
! 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_dp10,  id_dp50,  id_dp100, id_dp200,  id_dp250, id_dp300, &
            id_dp500, id_dp700, id_dp850, id_dp925,  id_dp1000

 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_hght3d, id_any_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(:)

     integer :: id_prer, id_prei, id_pres, id_preg
     integer :: id_qv_dt_gfdlmp, id_t_dt_gfdlmp, id_ql_dt_gfdlmp, id_qi_dt_gfdlmp
     integer :: id_u_dt_gfdlmp, id_v_dt_gfdlmp
     integer :: id_t_dt_phys, id_qv_dt_phys, id_ql_dt_phys, id_qi_dt_phys, id_u_dt_phys, id_v_dt_phys
     integer :: id_intqv, id_intql, id_intqi, id_intqr, id_intqs, id_intqg

     integer :: id_uw, id_vw, id_hw, id_qvw, id_qlw, id_qiw, id_o3w

     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


  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. 
     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
     
     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.

     integer, pointer :: grid_type

     logical, pointer :: nested
     logical, pointer :: regional
     logical :: bounded_domain

  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
!  -> moved to grid_tools

! Momentum (or KE) options:
   integer :: hord_mt = 9

   integer :: kord_mt = 8

   integer :: kord_wz = 8

   integer :: hord_vt = 9


! Heat & air mass (delp) transport options:
   integer :: hord_tm = 9

   integer :: hord_dp = 9

   integer :: kord_tm = -8
   integer :: hord_tr = 12

   integer :: kord_tr = 8
 
   real :: scale_z = 0.
   real :: w_max = 75.
   real :: z_min = 0.05
   real :: lim_fac = 1.0

   integer :: nord = 1

   integer :: nord_tr = 0

   real :: dddmp = 0.0

   real :: d2_bg = 0.0

   real :: d4_bg = 0.16

   real :: vtdm4 = 0.0

   real :: trdm2 = 0.0
   real :: d2_bg_k1 = 4.

   real :: d2_bg_k2 = 2.

   real :: d2_divg_max_k1 = 0.15
   real :: d2_divg_max_k2 = 0.08
   real :: damp_k_k1 = 0.2
   real :: damp_k_k2 = 0.12

   integer :: n_zs_filter = 0

   integer :: nord_zs_filter = 4

   logical :: full_zs_filter = .false.   

   logical :: rf_fast = .false.   

   logical :: consv_am  = .false.   


   logical :: do_sat_adj= .false.   ! 
   logical :: do_f3d = .false.   ! 
   logical :: no_dycore = .false.   

   logical :: convert_ke = .false.   
 
   logical :: do_vort_damp = .false.   

   logical :: use_old_omega = .true. 
! PG off centering:
   real :: beta = 0.0

#ifdef SW_DYNAMICS
   integer :: n_sponge = 0

   real :: d_ext = 0.
   
   integer :: nwat  = 0

   logical :: warm_start = .false.   

   logical :: inline_q = .true.    

   logical :: adiabatic = .true.   
#else
   integer :: n_sponge = 1

   real    :: d_ext = 0.02

   integer :: nwat  = 3

   logical :: warm_start = .true.   

   logical :: inline_q = .false.   
   logical :: adiabatic = .false.     
#endif
!-----------------------------------------------------------
! Grid shifting, rotation, and cube transformations:
!-----------------------------------------------------------
   real :: shift_fac = 18.
 
! Defaults for Schmidt transformation:
   logical :: do_schmidt = .false.   
   logical :: do_cube_transform = .false. 
   real(kind=R_GRID) :: stretch_fac = 1.

   real(kind=R_GRID) :: target_lat = -90.

   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

   integer :: n_split = 0

   real    :: fac_n_spl = 1.0
   real    :: fhouri    = 0.0

   integer :: m_split = 0
   integer :: k_split = 1

   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

   integer :: print_freq = 0

   logical :: write_3d_diags = .true.   
!------------------------------------------
! Model Domain parameters
!------------------------------------------
   integer :: npx
          
   integer :: npy

   integer :: npz
#ifdef USE_GFSL63
   character(24) :: npz_type = 'gfs'
#else
   character(24) :: npz_type = ''
#endif

   integer :: npz_rst = 0
                      
   integer :: ncnst = 0

   integer :: pnats = 0
  
   integer :: dnats = 0
   integer :: dnrts = -1

   integer :: ntiles = 1

   integer :: ndims = 2     ! Lat-Lon Dims for Grid in Radians
   integer :: nf_omega  = 1

   integer :: fv_sg_adj = -1

   real    :: sg_cutoff = -1

   integer :: na_init = 0
   

   logical :: nudge_dz = .false.    

   real :: p_ref = 1.e5   

   real :: dry_mass = 98290.

   integer :: nt_prog = 0
   integer :: nt_phys = 0
   real :: tau_h2o = 0.


   real :: delt_max = 1.

   real :: d_con = 0.

   real :: ke_bg = 0.

   real :: consv_te = 0.
 
   real :: tau = 0.

   real :: rf_cutoff = 30.e2   

   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 :: fill_gfs = .true. ! default behavior
   logical :: check_negative = .false.   
   logical :: non_ortho = .true.
   logical :: moist_phys = .true.     
   logical :: do_held_suarez = .false.   
   logical :: do_reed_physics = .false.
   logical :: reed_cond_only = .false.
   logical :: reproduce_sum = .true.   

   logical :: adjust_dry_mass = .false.    
 
   logical :: fv_debug = .false.   
   logical :: srf_init = .false.
   logical :: mountain = .true.   
   logical :: old_divg_damp = .false. 

   logical :: remap_t = .true.   

   logical :: z_tracer = .false.   

   logical :: fv_land = .false.   
!--------------------------------------------------------------------------------------
! The following options are useful for NWP experiments using datasets on the lat-lon grid
!--------------------------------------------------------------------------------------
   logical :: nudge = .false.   

   logical :: nudge_ic = .false.   

   logical :: ncep_ic = .false.   

   logical :: nggps_ic = .false.   

   logical :: ecmwf_ic = .false.   

   logical :: gfs_phil = .false.   

   logical :: agrid_vel_rst = .false.   

   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.   

   logical :: external_eta = .false.   

   logical :: read_increment = .false.   
! following are namelist parameters for Stochastic Energy Baskscatter dissipation estimate
   logical :: do_skeb  = .false.         
   integer :: skeb_npass  = 11
! 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.   

   logical :: use_hydro_pressure = .false.   

   logical :: do_uni_zfull = .false.   

   logical :: hybrid_z = .false.   
  
   logical :: make_nh = .false.   

   logical :: make_hybrid_z = .false.   

   logical :: nudge_qv = .false.   

   real :: add_noise = -1.

   logical :: butterfly_effect = .false.   

   integer :: a2b_ord = 4

   integer :: c2l_ord = 4

   integer :: nrows_blend = 0

   real(kind=R_GRID) :: dx_const = 1000.

   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., &  
                        deglat_start = -30., deglat_stop = 30.

   logical :: regional = .false.       

   integer :: bc_update_interval = 3

   logical :: regional_bcs_from_gsi = .false.   

   logical :: write_restart_with_bcs = .false.   

   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 nest_level_type 
     !Interpolation arrays for grid nesting
     logical                                :: on_level ! indicate if current processor on this level.
     logical                                :: do_remap_bc
     integer, allocatable, dimension(:,:,:) :: ind_h, ind_u, ind_v, ind_b ! I don't think these are necessary since BC interpolation is done locally
     real, allocatable, dimension(:,:,:) :: wt_h, wt_u, wt_v, wt_b
  end type nest_level_type

  type fv_nest_type

!nested grid flags:

     integer :: refinement = 3


     integer :: parent_tile = 1

     logical :: nested = .false.   
                                         
     integer :: nestbctype = 1
     integer :: nsponge = 0
     integer :: nestupdate = 7
       
     logical :: twowaynest = .false.   
     integer :: ioffset, joffset
     integer :: nlevel = 0 ! levels down from top-most domain

     integer :: nest_timestep = 0
     integer :: tracer_nest_timestep = 0
     real    :: s_weight = 1.e-6 
     logical :: first_step = .true.
     integer :: refinement_of_global = 1
     integer :: npx_global
     integer :: upoff = 1
     integer :: isu = -999, ieu = -1000, jsu = -999, jeu = -1000 
     real    :: update_blend = 1. ! option for controlling how much "blending" is done during two-way update
     logical, allocatable :: do_remap_bc(:)

     !nest_domain now a global structure defined in fv_mp_mod
     !type(nest_domain_type) :: nest_domain !Structure holding link from this grid to its parent
     !type(nest_domain_type), allocatable :: nest_domain_all(:)
     integer                :: num_nest_level ! number of nest levels.
     type(nest_level_type), allocatable :: nest(:) ! store data for each level.

     !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

     !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

     !points to same parent grid as does Atm%parent_grid
     type(fv_atmos_type), pointer :: parent_grid => null()


     !These are for tracer flux BCs
     logical :: do_flux_bcs, do_2way_flux_bcs
     type(restart_file_type) :: bcfile_ne, bcfile_sw

  end type fv_nest_type

  type phys_diag_type

     real, _allocatable :: phys_t_dt(:,:,:)
     real, _allocatable :: phys_qv_dt(:,:,:)
     real, _allocatable :: phys_ql_dt(:,:,:)
     real, _allocatable :: phys_qi_dt(:,:,:)
     real, _allocatable :: phys_u_dt(:,:,:)
     real, _allocatable :: phys_v_dt(:,:,:)

  end type phys_diag_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 = 3 !default

  end type fv_grid_bounds_type

  type fv_regional_bc_bounds_type

     integer :: is_north ,ie_north ,js_north ,je_north &
               ,is_south ,ie_south ,js_south ,je_south &
               ,is_east  ,ie_east  ,js_east  ,je_east  &
               ,is_west  ,ie_west  ,js_west  ,je_west

     integer :: is_north_uvs ,ie_north_uvs ,js_north_uvs ,je_north_uvs &
               ,is_south_uvs ,ie_south_uvs ,js_south_uvs ,je_south_uvs &
               ,is_east_uvs  ,ie_east_uvs  ,js_east_uvs  ,je_east_uvs  &
               ,is_west_uvs  ,ie_west_uvs  ,js_west_uvs  ,je_west_uvs   

     integer :: is_north_uvw ,ie_north_uvw ,js_north_uvw ,je_north_uvw &
               ,is_south_uvw ,ie_south_uvw ,js_south_uvw ,je_south_uvw &
               ,is_east_uvw  ,ie_east_uvw  ,js_east_uvw  ,je_east_uvw  &
               ,is_west_uvw  ,ie_west_uvw  ,js_west_uvw  ,je_west_uvw

  end type fv_regional_bc_bounds_type

  type fv_atmos_type

     logical :: allocated = .false.
     logical :: dummy = .false. ! same as grids_on_this_pe(n)
     integer :: grid_number = 1
     character(len=32) :: nml_filename = "input.nml"

     !Timestep-related variables.

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

#ifdef GFS_PHYS
     !--- DUMMY for backwards-compatibility. Will be removed
     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
    real, _allocatable :: v(:,:,:)    _null
    real, _allocatable :: pt(:,:,:)   _null
    real, _allocatable :: delp(:,:,:) _null
    real, _allocatable :: q(:,:,:,:)  _null
    real, _allocatable :: qdiag(:,:,:,:)  _null

!----------------------
! non-hydrostatic state:
!----------------------------------------------------------------------
    real, _allocatable ::     w(:,:,:)  _null
    real, _allocatable ::  delz(:,:,:)  _null
    real, _allocatable ::   ze0(:,:,:)  _null
    real, _allocatable ::  q_con(:,:,:) _null

!-----------------------------------------------------------------------
! Auxilliary pressure arrays:
! The 5 vars below can be re-computed from delp and ptop.
!-----------------------------------------------------------------------
! dyn_aux:
    real, _allocatable :: ps (:,:)      _null
    real, _allocatable :: pe (:,:,: )   _null
    real, _allocatable :: pk  (:,:,:)   _null
    real, _allocatable :: peln(:,:,:)   _null
    real, _allocatable :: pkz (:,:,:)   _null

! For phys coupling:
    real, _allocatable :: u_srf(:,:)    _null
    real, _allocatable :: v_srf(:,:)    _null
    real, _allocatable :: sgh(:,:)      _null
    real, _allocatable :: oro(:,:)      _null
    real, _allocatable :: ts(:,:)       _null
! For stochastic kinetic energy backscatter (SKEB)
    real, _allocatable :: diss_est(:,:,:) _null
 
!-----------------------------------------------------------------------
! Others:
!-----------------------------------------------------------------------
    real, _allocatable :: phis(:,:)     _null
    real, _allocatable :: omga(:,:,:)   _null
    real, _allocatable :: ua(:,:,:)     _null
    real, _allocatable :: va(:,:,:)     _null
    real, _allocatable :: uc(:,:,:)     _null
    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(fv_regional_bc_bounds_type) :: regional_bc_bounds

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

    type(domain2d) :: domain_for_coupler

     integer :: num_contact, npes_per_tile, global_tile, tile_of_mosaic, 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)

     type(phys_diag_type) :: phys_diag

  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, 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.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
    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
    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
    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 
    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=   npz_in !for ak, bk, which are 1D arrays and thus OK to allocate
       ndims_2d=   1 
       ncnst_2d=   1 
       nq_2d=   1 
    endif

!This should be set up in fv_mp_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

    !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%diss_est(isd:ied  ,jsd:jed  ,npz) ) 
    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(is:is, js:js  ,1) )
       allocate (  atm%ze0(is:is, js:js  ,1) )
    else
       allocate (    atm%w(isd:ied, jsd:jed  ,npz  ) )
       allocate ( atm%delz(is:ie, js:je  ,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

! 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(is,ie,js,je,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) = 0.
              atm%vc(i,j,k) = real_big
           enddo
        enddo
        do j=jsd, jed
           do i=isd, ied+1
               atm%v(i,j,k) = 0.
              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
              enddo
           enddo
           do j=js, je
              do i=is, ie
                 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
     do j=js, je
        do i=is, ie
           atm%ts(i,j) = 300.
           atm%phis(i,j) = real_big
        enddo
     enddo

    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
    
    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

    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-atm%ng:npx_2d  +atm%ng,1-atm%ng:npy_2d  +atm%ng,2,1))
          else
             allocate(atm%grid_global(1-atm%ng:npx_2d  +atm%ng,1-atm%ng:npy_2d  +atm%ng,2,1:6))
          endif
       end if
    endif


    !!Convenience pointers
    atm%gridstruct%nested    => atm%neststruct%nested
    atm%gridstruct%grid_type => atm%flagstruct%grid_type
    atm%flagstruct%grid_number => atm%grid_number
    atm%gridstruct%regional  => atm%flagstruct%regional
    atm%gridstruct%bounded_domain = atm%flagstruct%regional .or. atm%neststruct%nested
    if (atm%neststruct%nested) atm%neststruct%parent_grid => atm%parent_grid

    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

    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_mod