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


! <table>
! <tr>
!     <th>Module Name</th>
!     <th>Functions Included</th>
!   </tr>
!   <tr>
!     <td>boundary_mod</td>
!     <td>fill_nested_grid, nested_grid_BC, update_coarse_grid</td>
!   </tr>
!   <tr>
!     <td>constants_mod</td>
!     <td>kappa, pi=>pi_8, omega, rdgas, grav, rvgas, cp_air, radius</td>
!   </tr>
!   <tr>
!     <td>external_ic_mod</td>
!     <td>get_external_ic, get_cubed_sphere_terrain</td>
!   </tr>
!   <tr>
!     <td>field_manager_mod</td>
!     <td>MODEL_ATMOS</td>
!   </tr>
!   <tr>
!     <td>fms_mod</td>
!     <td>file_exist</td>
!   </tr>
!   <tr>
!     <td>fv_arrays_mod</td>
!     <td>fv_atmos_type, fv_nest_type, fv_grid_bounds_type, R_GRID</td>
!   </tr>
!   <tr>
!     <td>fv_control_mod</td>
!     <td>fv_init, fv_end, ngrids</td>
!   </tr>
!   <tr>
!     <td>fv_diagnostics_mod</td>
!     <td>prt_maxmin</td>
!   </tr>
!   <tr>
!     <td>fv_eta_mod</td>
!     <td>compute_dz_var, compute_dz_L32, set_hybrid_z</td>
!   </tr>
!   <tr>
!     <td>fv_grid_utils_mod</td>
!     <td>ptop_min, fill_ghost, g_sum, 
!         make_eta_level, cubed_to_latlon, great_circle_dist</td>
!   </tr>
!   <tr>
!     <td>fv_io_mod</td>
!     <td>fv_io_init, fv_io_read_restart, fv_io_write_restart, 
!        remap_restart, fv_io_register_restart, fv_io_register_nudge_restart, 
!        fv_io_register_restart_BCs, fv_io_register_restart_BCs_NH, fv_io_write_BCs, 
!        fv_io_read_BCs</td>
!   </tr>
!   <tr>
!     <td>fv_mp_mod</td>
!     <td>is_master, switch_current_Atm, mp_reduce_min, mp_reduce_max</td>
!   </tr>
!   <tr>
!     <td>fv_surf_map_mod</td>
!     <td>sgh_g, oro_g,del2_cubed_sphere, del4_cubed_sphere </td>
!   </tr>
!   <tr>
!     <td>fv_treat_da_inc_mod</td>
!     <td>read_da_inc</td>
!   </tr>
!   <tr>
!     <td>fv_timing_mod</td>
!     <td>timing_on, timing_off</td>
!   </tr>
!   <tr>
!     <td>fv_update_phys_mod</td>
!     <td>fv_update_phys</td>
!   </tr>
!   <tr>
!     <td>init_hydro_mod</td>
!     <td>p_var</td>
!   </tr>
!   <tr>
!     <td>mpp_mod</td>
!     <td>mpp_chksum, stdout, mpp_error, FATAL, NOTE, get_unit, mpp_sum,
!         mpp_get_current_pelist, mpp_set_current_pelist, mpp_send, mpp_recv, 
!         mpp_sync_self, mpp_npes, mpp_pe, mpp_sync</td>
!   </tr>
!   <tr>
!     <td>mpp_domains_mod</td>
!     <td>mpp_get_compute_domain, mpp_get_data_domain, mpp_get_global_domain, 
!         mpp_update_domains, domain2d, DGRID_NE, CENTER, CORNER, NORTH, EAST,
!         mpp_get_C2F_index, WEST, SOUTH, mpp_global_field</td>
!   </tr>
!   <tr>
!     <td>mpp_parameter_mod</td>
!     <td>EUPDATE, WUPDATE, SUPDATE, NUPDATE</td>
!   </tr>
!   <tr>
!     <td>time_manager_mod</td>
!     <td>time_type, get_time, set_time, operator(+), operator(-)</td>
!   </tr>
!   <tr>
!     <td>tracer_manager_mod</td>
!     <td>get_tracer_index, get_tracer_names</td>
!   </tr>
!   <tr>
!     <td>test_cases_mod</td>
!     <td>test_case, alpha, init_case, init_double_periodic, init_latlon</td>
!   </tr>
! </table>


  use constants_mod,       only: kappa, pi=>pi_8, omega, rdgas, grav, rvgas, cp_air, radius
  use fv_arrays_mod,       only: fv_atmos_type, fv_nest_type, fv_grid_bounds_type, r_grid
  use fv_io_mod,           only: fv_io_init, fv_io_read_restart, fv_io_write_restart, &
                                 remap_restart, fv_io_register_restart, fv_io_register_nudge_restart, &
                                 fv_io_register_restart_bcs, fv_io_write_bcs, fv_io_read_bcs
  use fv_grid_utils_mod,   only: ptop_min, fill_ghost, g_sum, &
                                 make_eta_level, cubed_to_latlon, great_circle_dist
  use fv_diagnostics_mod,  only: prt_maxmin
  use init_hydro_mod,      only: p_var
  use mpp_domains_mod,     only: mpp_update_domains, domain2d, dgrid_ne
  use mpp_mod,             only: mpp_chksum, stdout, mpp_error, fatal, note
  use mpp_mod,             only: get_unit, mpp_sum, mpp_broadcast
  use mpp_mod,             only: mpp_get_current_pelist, mpp_npes, mpp_set_current_pelist
  use test_cases_mod,      only: alpha, init_case, init_double_periodic!, init_latlon
  use fv_mp_mod,           only: is_master, mp_reduce_min, mp_reduce_max, corners_ydir => ydir, fill_corners, tile_fine, global_nest_domain
  use fv_surf_map_mod,     only: sgh_g, oro_g
  use tracer_manager_mod,  only: get_tracer_names
  use field_manager_mod,   only: model_atmos
  use external_ic_mod,     only: get_external_ic
  use fv_eta_mod,          only: compute_dz_var, compute_dz_l32, set_hybrid_z
  use fv_surf_map_mod,     only: del2_cubed_sphere, del4_cubed_sphere
  use boundary_mod,        only: fill_nested_grid, nested_grid_bc, update_coarse_grid
  use tracer_manager_mod,  only: get_tracer_index
  use field_manager_mod,   only: model_atmos
  use fv_timing_mod,       only: timing_on, timing_off
  use mpp_domains_mod,     only: mpp_get_compute_domain, mpp_get_data_domain, mpp_get_global_domain
  use mpp_mod,             only: mpp_send, mpp_recv, mpp_sync_self, mpp_set_current_pelist, mpp_get_current_pelist, mpp_npes, mpp_pe, mpp_sync
  use mpp_domains_mod,     only: center, corner, north, east,  mpp_get_c2f_index, west, south
  use mpp_domains_mod,     only: mpp_global_field
  use fms_mod,             only: file_exist
  use fv_treat_da_inc_mod, only: read_da_inc
  use fv_regional_mod,     only: write_full_fields
#ifdef MULTI_GASES
  use multi_gases_mod,  only:  virq
#endif

  implicit none
  private

  public :: fv_restart_init, fv_restart_end, fv_restart, fv_write_restart

  real(kind=R_GRID), parameter :: cnst_0p20=0.20d0
  !--- private data type
  logical                       :: module_is_initialized = .false.

contains 


  subroutine fv_restart_init()
    call fv_io_init()
    module_is_initialized = .true.
  end subroutine fv_restart_init

  subroutine fv_restart(fv_domain, Atm, dt_atmos, seconds, days, cold_start, grid_type, this_grid)
    type(domain2d),      intent(inout) :: fv_domain
    type(fv_atmos_type), intent(inout) :: Atm(:)
    real,                intent(in)    :: dt_atmos
    integer,             intent(out)   :: seconds
    integer,             intent(out)   :: days
    logical,             intent(inout)    :: cold_start
    integer,             intent(in)    :: grid_type, this_grid

    integer :: i, j, k, n, ntileMe, nt, iq
    integer :: isc, iec, jsc, jec, ncnst, ntprog, ntdiag
    integer :: isd, ied, jsd, jed, npz
    integer isd_p, ied_p, jsd_p, jed_p, isc_p, iec_p, jsc_p, jec_p, isg, ieg, jsg,jeg, npx_p, npy_p
    real, allocatable :: g_dat(:,:,:)

    integer :: unit
    real, allocatable :: dz1(:)
    real rgrav, f00, ztop, pertn, ph
    logical :: hybrid
    character(len=128):: tname, errstring, fname, tracer_name
    character(len=120):: fname_ne, fname_sw
    character(len=3) :: gn

    integer :: npts, sphum
    integer, allocatable :: pelist(:), smoothed_topo(:)
    real    :: sumpertn
    real    :: zvir

    integer :: i_butterfly, j_butterfly
    logical :: do_read_restart = .false.
    logical :: do_read_restart_bc = .false.
    integer, allocatable :: ideal_test_case(:), new_nest_topo(:)

    rgrav = 1. / grav

    if(.not.module_is_initialized) call mpp_error(fatal, 'You must call fv_restart_init.')

    ntileme = size(atm(:))
    allocate(smoothed_topo(ntileme))
    smoothed_topo(:) = 0
    allocate(ideal_test_case(ntileme))
    ideal_test_case(:) = 0
    allocate(new_nest_topo(ntileme))
    new_nest_topo(:) = 0

    do n = 1, ntileme

       isd = atm(n)%bd%isd
       ied = atm(n)%bd%ied
       jsd = atm(n)%bd%jsd
       jed = atm(n)%bd%jed
       isc = atm(n)%bd%isc
       iec = atm(n)%bd%iec
       jsc = atm(n)%bd%jsc
       jec = atm(n)%bd%jec
       ncnst = atm(n)%ncnst
       if( is_master() ) write(*,*) 'in fv_restart ncnst=', ncnst
       npz = atm(n)%npz
       ntprog = size(atm(n)%q,4)
       ntdiag = size(atm(n)%qdiag,4)

!!$       if (is_master()) then
!!$          print*, 'FV_RESTART: ', n, cold_start_grids(n)
!!$       endif

       !1. sort out restart, external_ic, and cold-start (idealized)
       if (atm(n)%neststruct%nested) then
          write(fname,   '(A, I2.2, A)') 'INPUT/fv_core.res.nest', atm(n)%grid_number, '.nc'
          write(fname_ne,'(A, I2.2, A)') 'INPUT/fv_BC_ne.res.nest', atm(n)%grid_number, '.nc'
          write(fname_sw,'(A, I2.2, A)') 'INPUT/fv_BC_sw.res.nest', atm(n)%grid_number, '.nc'
          if (is_master()) print*, 'Searching for nested grid BC files ', trim(fname_ne), ' ', trim(fname_sw)
          do_read_restart = file_exist(fname, atm(n)%domain)
          do_read_restart_bc = file_exist(fname_ne, atm(n)%domain) .and. file_exist(fname_sw, atm(n)%domain)
          if (is_master()) then
             print*, 'FV_RESTART: ', n, do_read_restart, do_read_restart_bc
             if (.not. do_read_restart_bc) write(*,*) 'BC files not found, re-generating nested grid boundary conditions'
          endif
          atm(n)%neststruct%first_step = .not. do_read_restart_bc
       else
          fname='INPUT/fv_core.res.nc'
          do_read_restart = file_exist('INPUT/fv_core.res.nc') .or. file_exist('INPUT/fv_core.res.tile1.nc') 
          if (is_master()) print*, 'FV_RESTART: ', n, do_read_restart, do_read_restart_bc
       endif

       !2. Register restarts
       !--- call fv_io_register_restart to register restart field to be written out in fv_io_write_restart
       if ( n==this_grid ) call fv_io_register_restart(atm(n)%domain,atm(n:n))
       !if (Atm(n)%neststruct%nested) call fv_io_register_restart_BCs(Atm(n)) !TODO put into fv_io_register_restart


       !3preN. Topography BCs for nest, including setup for blending
       
       if (atm(n)%neststruct%nested) then
          if (.not. allocated(pelist)) then
             allocate(pelist(0:mpp_npes()-1))
             call mpp_get_current_pelist(pelist)      
          endif
          call mpp_set_current_pelist() !global
          call mpp_broadcast(atm(n)%flagstruct%external_ic,atm(n)%pelist(1))
          call mpp_sync()
          call mpp_set_current_pelist(pelist)
          if ( ( smoothed_topo(atm(n)%parent_grid%grid_number) > 0 .or. &
                 .not. do_read_restart_bc .or. &
                 atm(n)%flagstruct%external_ic  ) ) then
             new_nest_topo(n) = 1
             if (n==this_grid) then
                
                call fill_nested_grid_topo(atm(n), n==this_grid)
                call fill_nested_grid_topo_halo(atm(n), n==this_grid) !TODO can we combine these?
                call nested_grid_bc(atm(n)%ps, atm(n)%parent_grid%ps, global_nest_domain, &
                     atm(n)%neststruct%ind_h, atm(n)%neststruct%wt_h, 0, 0, &
                     atm(n)%npx, atm(n)%npy, atm(n)%bd, 1, atm(n)%npx-1, 1, atm(n)%npy-1)
                
             elseif (this_grid==atm(n)%parent_grid%grid_number) then !this_grid is grid n's parent
                
                call fill_nested_grid_topo(atm(n), n==this_grid)
                call fill_nested_grid_topo_halo(atm(n), n==this_grid) !TODO can we combine these?
                !call mpp_get_data_domain( Atm(n)%parent_grid%domain, isd, ied, jsd, jed)
                call nested_grid_bc(atm(n)%parent_grid%ps, global_nest_domain, 0, 0, n-1)
                !Atm(n)%ps, Atm(n)%parent_grid%ps, global_nest_domain, &
                !Atm(n)%neststruct%ind_h, Atm(n)%neststruct%wt_h, 0, 0, &
                !Atm(n)%npx, Atm(n)%npy, Atm(n)%bd, isd, ied, jsd, jed, proc_in=n==this_grid)         

             endif

          endif
       endif

       !This call still appears to be necessary to get isd, etc. correct
       !call switch_current_Atm(Atm(n)) !TODO should NOT be necessary now that we manually set isd, etc.

       !--- call fv_io_register_restart to register restart field to be written out in fv_io_write_restart
       !if (n==this_grid) call fv_io_register_restart(Atm(n)%domain,Atm(n:n))
       !if (Atm(n)%neststruct%nested) call fv_io_register_restart_BCs(Atm(n)) !TODO put into fv_io_register_restart

       if (n==this_grid) then
       
          !3. External_ic
          if (atm(n)%flagstruct%external_ic) then
             if( is_master() ) write(*,*) 'Calling get_external_ic'
             call get_external_ic(atm(n), atm(n)%domain, .not. do_read_restart, dt_atmos) 
             if( is_master() ) write(*,*) 'IC generated from the specified external source'

             !4. Restart
          elseif (do_read_restart) then

             if ( atm(n)%flagstruct%npz_rst /= 0 .and. atm(n)%flagstruct%npz_rst /= atm(n)%npz ) then
                !Remap vertically the prognostic variables for the chosen vertical resolution
                if( is_master() ) then
                   write(*,*) ' '
                   write(*,*) '***** Important Note from FV core ********************'
                   write(*,*) 'Remapping dynamic IC from', atm(n)%flagstruct%npz_rst, 'levels to ', atm(n)%npz,'levels'
                   write(*,*) '***** End Note from FV core **************************'
                   write(*,*) ' '
                endif
                call remap_restart( atm(n)%domain, atm(n:n) )
                if( is_master() ) write(*,*) 'Done remapping dynamical IC'
             else
                if( is_master() ) write(*,*) 'Warm starting, calling fv_io_restart'
                call fv_io_read_restart(atm(n)%domain,atm(n:n))
                !====== PJP added DA functionality ======
                if (atm(n)%flagstruct%read_increment) then
                   ! print point in middle of domain for a sanity check
                   i = (atm(n)%bd%isc + atm(n)%bd%iec)/2
                   j = (atm(n)%bd%jsc + atm(n)%bd%jec)/2
                   k = atm(n)%npz/2
                   if( is_master() ) write(*,*) 'Calling read_da_inc',atm(n)%pt(i,j,k)
                   call read_da_inc(atm(n), atm(n)%domain, atm(n)%bd, atm(n)%npz, atm(n)%ncnst, &
                        atm(n)%u, atm(n)%v, atm(n)%q, atm(n)%delp, atm(n)%pt, atm(n)%delz, isd, jsd, ied, jed, &
                        isc, jsc, iec, jec )
                   if( is_master() ) write(*,*) 'Back from read_da_inc',atm(n)%pt(i,j,k)
                endif
                !====== end PJP added DA functionailty======
             endif
             
             seconds = 0; days = 0   ! Restart needs to be modified to record seconds and days.

             if (atm(n)%neststruct%nested) then
                if ( atm(n)%flagstruct%npz_rst /= 0 .and. atm(n)%flagstruct%npz_rst /= npz ) then
                   call mpp_error(fatal, "Remap-restart not implemented for nests.")
                endif
                if (do_read_restart_bc) call fv_io_read_bcs(atm(n))
                call mpp_update_domains(atm(n)%u, atm(n)%v, atm(n)%domain, gridtype=dgrid_ne, complete=.true.)
             endif

             if ( atm(n)%flagstruct%mountain ) then
                ! !!! Additional terrain filter -- should not be called repeatedly !!!
                if ( atm(n)%flagstruct%n_zs_filter > 0 ) then
                   if ( atm(n)%flagstruct%nord_zs_filter == 2 ) then
                      !!! TODO: move this block into its own routine or CLEAN UP these subroutine calls
                      call del2_cubed_sphere(atm(n)%npx, atm(n)%npy, atm(n)%phis, &
                           atm(n)%gridstruct%area_64, atm(n)%gridstruct%dx, atm(n)%gridstruct%dy,   &
                           atm(n)%gridstruct%dxc, atm(n)%gridstruct%dyc, atm(n)%gridstruct%sin_sg, &
                           atm(n)%flagstruct%n_zs_filter, cnst_0p20*atm(n)%gridstruct%da_min, &
                           .false., oro_g, atm(n)%gridstruct%bounded_domain, atm(n)%domain, atm(n)%bd)
                      if ( is_master() ) write(*,*) 'Warning !!! del-2 terrain filter has been applied ', &
                           atm(n)%flagstruct%n_zs_filter, ' times'
                   else if( atm(n)%flagstruct%nord_zs_filter == 4 ) then
                      call del4_cubed_sphere(atm(n)%npx, atm(n)%npy, atm(n)%phis, atm(n)%gridstruct%area_64, &
                           atm(n)%gridstruct%dx, atm(n)%gridstruct%dy,   &
                           atm(n)%gridstruct%dxc, atm(n)%gridstruct%dyc, atm(n)%gridstruct%sin_sg, &
                           atm(n)%flagstruct%n_zs_filter, .false., oro_g, atm(n)%gridstruct%bounded_domain, &
                           atm(n)%domain, atm(n)%bd)
                      if ( is_master() ) write(*,*) 'Warning !!! del-4 terrain filter has been applied ', &
                           atm(n)%flagstruct%n_zs_filter, ' times'
                   endif
                endif
                call mpp_update_domains( atm(n)%phis, atm(n)%domain, complete=.true. )
             else
                atm(n)%phis = 0.
                if( is_master() ) write(*,*) 'phis set to zero'
             endif !mountain



             !5. Idealized test case
          else
             
             ideal_test_case(n) = 1

             if ( atm(n)%flagstruct%make_hybrid_z ) then
                hybrid = .false.
             else
                hybrid = atm(n)%flagstruct%hybrid_z
             endif
             if (grid_type < 4) then
                if ( .not. atm(n)%flagstruct%external_ic ) then
                   call init_case(atm(n)%u,atm(n)%v,atm(n)%w,atm(n)%pt,atm(n)%delp,atm(n)%q, &
                        atm(n)%phis, atm(n)%ps,atm(n)%pe, atm(n)%peln,atm(n)%pk,atm(n)%pkz, &
                        atm(n)%uc,atm(n)%vc, atm(n)%ua,atm(n)%va,        & 
                        atm(n)%ak, atm(n)%bk, atm(n)%gridstruct, atm(n)%flagstruct,&
                        atm(n)%npx, atm(n)%npy, npz, atm(n)%ng, &
                        ncnst, atm(n)%flagstruct%nwat,  &
                        atm(n)%flagstruct%ndims, atm(n)%flagstruct%ntiles, &
                        atm(n)%flagstruct%dry_mass, &
                        atm(n)%flagstruct%mountain,       &
                        atm(n)%flagstruct%moist_phys, atm(n)%flagstruct%hydrostatic, &
                        hybrid, atm(n)%delz, atm(n)%ze0, &
                        atm(n)%flagstruct%adiabatic, atm(n)%ks, atm(n)%neststruct%npx_global, &
                        atm(n)%ptop, atm(n)%domain, atm(n)%tile_of_mosaic, atm(n)%bd)
                endif
             elseif (grid_type == 4) then
                call init_double_periodic(atm(n)%u,atm(n)%v,atm(n)%w,atm(n)%pt, &
                     atm(n)%delp,atm(n)%q,atm(n)%phis, atm(n)%ps,atm(n)%pe, &
                     atm(n)%peln,atm(n)%pk,atm(n)%pkz, &
                     atm(n)%uc,atm(n)%vc, atm(n)%ua,atm(n)%va,        & 
                     atm(n)%ak, atm(n)%bk, &
                     atm(n)%gridstruct, atm(n)%flagstruct, &
                     atm(n)%npx, atm(n)%npy, npz, atm(n)%ng, &
                     ncnst, atm(n)%flagstruct%nwat,  &
                     atm(n)%flagstruct%ndims, atm(n)%flagstruct%ntiles, &
                     atm(n)%flagstruct%dry_mass, atm(n)%flagstruct%mountain, &
                     atm(n)%flagstruct%moist_phys, atm(n)%flagstruct%hydrostatic, &
                     hybrid, atm(n)%delz, atm(n)%ze0, atm(n)%ks, atm(n)%ptop, &
                     atm(n)%domain, atm(n)%tile_of_mosaic, atm(n)%bd)
                if( is_master() ) write(*,*) 'Doubly Periodic IC generated'
             elseif (grid_type == 5 .or. grid_type == 6) then
                call mpp_error(fatal, "Idealized test cases for grid_type == 5,6 (global lat-lon) grid not supported")
             endif

             !Turn this off on the nested grid if you are just interpolating topography from the coarse grid!
             !These parameters are needed in LM3/LM4, and are communicated through restart files 
             if ( atm(n)%flagstruct%fv_land ) then
                do j=jsc,jec
                   do i=isc,iec
                      atm(n)%sgh(i,j) = sgh_g(i,j)
                      atm(n)%oro(i,j) = oro_g(i,j)
                   enddo
                enddo
             endif

          endif !external_ic vs. restart vs. idealized


       endif !n==this_grid


          !!!! NOT NEEDED??
          !Currently even though we do fill in the nested-grid IC from
          ! init_case or external_ic we appear to overwrite it using
          !  coarse-grid data
!!$          if (Atm(n)%neststruct%nested) then
!!$             if (.not. Atm(n)%flagstruct%external_ic .and.  .not. Atm(n)%flagstruct%nggps_ic .and. grid_type < 4 ) then
!!$                call fill_nested_grid_data(Atm(n:n))
!!$             endif
!!$          end if

!       endif  !end cold_start check

       !5n. Nesting setup (part I)

       !Broadcast data for nesting
       if (ntileme > 1) then
          if (.not. allocated(pelist)) then
             allocate(pelist(0:mpp_npes()-1))
             call mpp_get_current_pelist(pelist)      
          endif

          call mpp_set_current_pelist()!global
          !for remap BCs
          call mpp_broadcast(atm(n)%ptop,atm(n)%pelist(1))
          call mpp_broadcast(atm(n)%ak,atm(n)%npz+1,atm(n)%pelist(1))
          call mpp_broadcast(atm(n)%bk,atm(n)%npz+1,atm(n)%pelist(1))
          !smoothed_topo
          call mpp_broadcast(smoothed_topo(n),atm(n)%pelist(1))

          call mpp_sync()
          call mpp_set_current_pelist(pelist)


          if (atm(n)%neststruct%nested) then
             atm(n)%neststruct%do_remap_BC(ntileme) = .false.

             if (atm(n)%npz /= atm(n)%parent_grid%npz) then
                atm(n)%neststruct%do_remap_BC(n) = .true.
             else
                do k=1,atm(n)%npz+1
                   if (atm(n)%ak(k) /= atm(n)%parent_grid%ak(k)) then
                      atm(n)%neststruct%do_remap_BC(n) = .true.
                      exit
                   endif
                   if (atm(n)%bk(k) /= atm(n)%parent_grid%bk(k)) then
                      atm(n)%neststruct%do_remap_BC(n) = .true.
                      exit
                   endif
                enddo
             endif

             atm(n)%parent_grid%neststruct%do_remap_BC(n) = atm(n)%neststruct%do_remap_BC(n)
             if (is_master() .and. n==this_grid) then
                if (atm(n)%neststruct%do_remap_BC(n)) then
                   print*, ' Remapping BCs ENABLED on grid', n
                else
                   print*, ' Remapping BCs DISABLED (not necessary) on grid', n             
                endif
                write(*,'(A, I3, A, F8.2, A)') ' Nested grid ', n, ',  ptop = ', atm(n)%ak(1), ' Pa'
                write(*,'(A, I3, A, F8.2, A)') ' Parent grid ', n, ',  ptop = ', atm(n)%parent_grid%ak(1), ' Pa'
                if (atm(n)%ak(1) < atm(n)%parent_Grid%ak(1)) then
                   print*, ' WARNING nested grid top above parent grid top. May have problems with remapping BCs.'
                endif
             endif
          endif

       endif

    end do !break cycling loop to finish nesting setup


    do n = ntileme,1,-1
       if (new_nest_topo(n) > 0 ) then
          call twoway_topo_update(atm(n), n==this_grid)
       endif
    end do

    !6. Data Setup
    do n = 1, ntileme

       if (n/=this_grid) cycle

       isd = atm(n)%bd%isd
       ied = atm(n)%bd%ied
       jsd = atm(n)%bd%jsd
       jed = atm(n)%bd%jed
       isc = atm(n)%bd%isc
       iec = atm(n)%bd%iec
       jsc = atm(n)%bd%jsc
       jec = atm(n)%bd%jec
       ncnst = atm(n)%ncnst
       if( is_master() ) write(*,*) 'in fv_restart ncnst=', ncnst
       npz = atm(n)%npz
       ntprog = size(atm(n)%q,4)
       ntdiag = size(atm(n)%qdiag,4)


       if ( ideal_test_case(n) == 0 ) then
#ifdef SW_DYNAMICS
          atm(n)%pt(:,:,:)=1.
#else
          if ( .not.atm(n)%flagstruct%hybrid_z ) then
             if(atm(n)%ptop/=atm(n)%ak(1)) call mpp_error(fatal,'FV restart: ptop not equal Atm(n)%ak(1)')
          else
             atm(n)%ptop = atm(n)%ak(1);  atm(n)%ks = 0
          endif
          call p_var(npz,         isc,         iec,       jsc,     jec,   atm(n)%ptop,     ptop_min,  &
               atm(n)%delp, atm(n)%delz, atm(n)%pt, atm(n)%ps, atm(n)%pe, atm(n)%peln,   &
               atm(n)%pk,   atm(n)%pkz, kappa, atm(n)%q, atm(n)%ng, &
               ncnst,  atm(n)%gridstruct%area_64, atm(n)%flagstruct%dry_mass,  &
               atm(n)%flagstruct%adjust_dry_mass,  atm(n)%flagstruct%mountain, &
               atm(n)%flagstruct%moist_phys,  atm(n)%flagstruct%hydrostatic, &
               atm(n)%flagstruct%nwat, atm(n)%domain, atm(1)%flagstruct%adiabatic, atm(n)%flagstruct%make_nh)
#endif
          if ( grid_type < 7 .and. grid_type /= 4 ) then
             ! Fill big values in the non-existing corner regions:
             !          call fill_ghost(Atm(n)%phis, Atm(n)%npx, Atm(n)%npy, big_number)
             do j=jsd,jed+1
                do i=isd,ied+1
                   atm(n)%gridstruct%fc(i,j) = 2.*omega*( -cos(atm(n)%gridstruct%grid(i,j,1))*cos(atm(n)%gridstruct%grid(i,j,2))*sin(alpha) + &
                        sin(atm(n)%gridstruct%grid(i,j,2))*cos(alpha) )
                enddo
             enddo
             do j=jsd,jed
                do i=isd,ied
                   atm(n)%gridstruct%f0(i,j) = 2.*omega*( -cos(atm(n)%gridstruct%agrid(i,j,1))*cos(atm(n)%gridstruct%agrid(i,j,2))*sin(alpha) + &
                        sin(atm(n)%gridstruct%agrid(i,j,2))*cos(alpha) )
                enddo
             enddo
          else
             f00 = 2.*omega*sin(atm(n)%flagstruct%deglat/180.*pi)
             do j=jsd,jed+1
                do i=isd,ied+1
                   atm(n)%gridstruct%fc(i,j) = f00
                enddo
             enddo
             do j=jsd,jed
                do i=isd,ied
                   atm(n)%gridstruct%f0(i,j) = f00
                enddo
             enddo
          endif
          call mpp_update_domains( atm(n)%gridstruct%f0, atm(n)%domain )
          if ( atm(n)%gridstruct%cubed_sphere .and. (.not. atm(n)%gridstruct%bounded_domain))then
             call fill_corners(atm(n)%gridstruct%f0, atm(n)%npx, atm(n)%npy, corners_ydir)
          endif
       endif


!---------------------------------------------------------------------------------------------
! Transform the (starting) Eulerian vertical coordinate from sigma-p to hybrid_z
     if ( atm(n)%flagstruct%hybrid_z ) then
       if ( atm(n)%flagstruct%make_hybrid_z ) then
          allocate ( dz1(npz) )
          if( npz==32 ) then
              call compute_dz_l32(npz, ztop, dz1)
          else
              ztop = 45.e3
              call compute_dz_var(npz, ztop, dz1)
          endif
          call set_hybrid_z(isc, iec, jsc, jec, atm(n)%ng, npz, ztop, dz1, rgrav,  &
                            atm(n)%phis, atm(n)%ze0)
          deallocate ( dz1 )
!         call prt_maxmin('ZE0', Atm(n)%ze0,  isc, iec, jsc, jec, 0, npz, 1.E-3)
!         call prt_maxmin('DZ0', Atm(n)%delz, isc, iec, jsc, jec, 0, npz, 1.   )
       endif
!      call make_eta_level(npz, Atm(n)%pe, area, Atm(n)%ks, Atm(n)%ak, Atm(n)%bk, Atm(n)%ptop)
     endif
!---------------------------------------------------------------------------------------------

     if (atm(n)%flagstruct%add_noise > 0.) then
        write(errstring,'(A, E16.9)') "Adding thermal noise of amplitude ", atm(n)%flagstruct%add_noise
        call mpp_error(note, errstring)
        call random_seed
        npts = 0
        sumpertn = 0.
        do k=1,npz
        do j=jsc,jec
        do i=isc,iec
           call random_number(pertn)
           atm(n)%pt(i,j,k) = atm(n)%pt(i,j,k) + pertn*atm(n)%flagstruct%add_noise
           npts = npts + 1
           sumpertn = sumpertn + pertn*atm(n)%flagstruct%add_noise ** 2
        enddo
        enddo
        enddo
        call mpp_update_domains(atm(n)%pt, atm(n)%domain)
        call mpp_sum(sumpertn)
        call mpp_sum(npts)
        write(errstring,'(A, E16.9)') "RMS added noise: ", sqrt(sumpertn/npts)
        call mpp_error(note, errstring)
     endif

     if (atm(n)%flagstruct%butterfly_effect) then
        if (n==1 .and. atm(n)%tile_of_mosaic == 1) then
           i_butterfly = atm(n)%npx / 2
           j_butterfly = atm(n)%npy / 2
           if (isc <= i_butterfly .and. i_butterfly <= iec) then
           if (jsc <= j_butterfly .and. j_butterfly <= jec) then

              write(*,'(A, I0, A, I0)') "Adding butterfly effect at (i,j) ", i_butterfly, ", ", j_butterfly
              write(*,'(A, E24.17)') "pt (before) :", atm(n)%pt(i_butterfly,j_butterfly,atm(n)%npz)

              atm(n)%pt(i_butterfly,j_butterfly,atm(n)%npz) = nearest(atm(n)%pt(i_butterfly,j_butterfly,atm(n)%npz), -1.0)

              write(*,'(A, E24.17)') "pt (after)  :", atm(n)%pt(i_butterfly,j_butterfly,atm(n)%npz)

           endif
           endif
        endif
      endif
      if (atm(n)%flagstruct%fv_sg_adj > 0 .and. atm(n)%flagstruct%sg_cutoff > 0) then
        !Choose n_sponge from first reference level above sg_cutoff
        do k=1,npz
           ph = atm(n)%ak(k+1) +  atm(n)%bk(k+1)*atm(n)%flagstruct%p_ref
           if (ph > atm(n)%flagstruct%sg_cutoff) exit
        enddo
        atm(n)%flagstruct%n_sponge = min(k,npz)
        write(errstring,'(A, I3, A)') ' Override n_sponge: applying 2dz filter to ', k , ' levels'
        call mpp_error(note, errstring)   
      endif

      if (atm(n)%grid_number > 1) then
         write(gn,'(A2, I1)') " g", atm(n)%grid_number
      else
         gn = ''
      end if

      unit = stdout()
      !!!NOTE: Checksums not yet working in stand-alone regional model!!
      write(unit,*)
      write(unit,*) 'fv_restart u   ', trim(gn),' = ', mpp_chksum(atm(n)%u(isc:iec,jsc:jec,:))
      write(unit,*) 'fv_restart v   ', trim(gn),' = ', mpp_chksum(atm(n)%v(isc:iec,jsc:jec,:))
      if ( .not.atm(n)%flagstruct%hydrostatic )   &
        write(unit,*) 'fv_restart w   ', trim(gn),' = ', mpp_chksum(atm(n)%w(isc:iec,jsc:jec,:))
      write(unit,*) 'fv_restart delp', trim(gn),' = ', mpp_chksum(atm(n)%delp(isc:iec,jsc:jec,:))
      write(unit,*) 'fv_restart phis', trim(gn),' = ', mpp_chksum(atm(n)%phis(isc:iec,jsc:jec))

#ifdef SW_DYNAMICS
      call prt_maxmin('H ', atm(n)%delp, isc, iec, jsc, jec, atm(n)%ng, 1, rgrav)
#else
      write(unit,*) 'fv_restart pt  ', trim(gn),' = ', mpp_chksum(atm(n)%pt(isc:iec,jsc:jec,:))
      if (ntprog>0) &
           write(unit,*) 'fv_restart q(prog) nq  ', trim(gn),' =',ntprog, mpp_chksum(atm(n)%q(isc:iec,jsc:jec,:,:))
      if (ntdiag>0) &
           write(unit,*) 'fv_restart q(diag) nq  ', trim(gn),' =',ntdiag, mpp_chksum(atm(n)%qdiag(isc:iec,jsc:jec,:,:))
      do iq=1,min(17, ntprog)     ! Check up to 17 tracers
        call get_tracer_names(model_atmos, iq, tracer_name)
        write(unit,*) 'fv_restart '//trim(tracer_name)//' = ', mpp_chksum(atm(n)%q(isc:iec,jsc:jec,:,iq))
      enddo

!---------------
! Check Min/Max:
!---------------
      call pmaxmn_g('ZS', atm(n)%phis, isc, iec, jsc, jec, 1, rgrav, atm(n)%gridstruct%area_64, atm(n)%domain)
      call pmaxmn_g('PS', atm(n)%ps,   isc, iec, jsc, jec, 1, 0.01,  atm(n)%gridstruct%area_64, atm(n)%domain)
      call pmaxmn_g('T ', atm(n)%pt,   isc, iec, jsc, jec, npz, 1.,  atm(n)%gridstruct%area_64, atm(n)%domain)

! Check tracers:
      do i=1, ntprog
          call get_tracer_names ( model_atmos, i, tname )
          call pmaxmn_g(trim(tname), atm(n)%q(isd:ied,jsd:jed,1:npz,i:i), isc, iec, jsc, jec, npz, &
                        1., atm(n)%gridstruct%area_64, atm(n)%domain)
      enddo
#endif
      call prt_maxmin('U ', atm(n)%u(isc:iec,jsc:jec,1:npz), isc, iec, jsc, jec, 0, npz, 1.)
      call prt_maxmin('V ', atm(n)%v(isc:iec,jsc:jec,1:npz), isc, iec, jsc, jec, 0, npz, 1.)

      if ( (.not.atm(n)%flagstruct%hydrostatic) .and. atm(n)%flagstruct%make_nh ) then
         call mpp_error(note, "  Initializing w to 0")
         atm(n)%w = 0.
         sphum = get_tracer_index(model_atmos, 'sphum')
         if ( .not.atm(n)%flagstruct%hybrid_z ) then
            if (atm(n)%flagstruct%adiabatic .or. sphum < 0) then
             zvir = 0.
            else
             zvir = rvgas/rdgas - 1.
            endif
             do k=1,npz
                do j=jsc,jec
                   do i=isc,iec
                      atm(n)%delz(i,j,k) = (rdgas*rgrav)*atm(n)%pt(i,j,k)*(1.+zvir*atm(n)%q(i,j,k,sphum))*(atm(n)%peln(i,k,j)-atm(n)%peln(i,k+1,j))
                   enddo
                enddo
             enddo
         endif
      endif

      if ( .not.atm(n)%flagstruct%hydrostatic )   &
      call pmaxmn_g('W ', atm(n)%w, isc, iec, jsc, jec, npz, 1., atm(n)%gridstruct%area_64, atm(n)%domain)

      if (is_master()) write(unit,*)

!--------------------------------------------
! Initialize surface winds for flux coupler:
!--------------------------------------------
    if ( .not. atm(n)%flagstruct%srf_init ) then
         call cubed_to_latlon(atm(n)%u, atm(n)%v, atm(n)%ua, atm(n)%va, &
              atm(n)%gridstruct, &
              atm(n)%npx, atm(n)%npy, npz, 1, &              
              atm(n)%gridstruct%grid_type, atm(n)%domain, &
              atm(n)%gridstruct%bounded_domain, atm(n)%flagstruct%c2l_ord, atm(n)%bd)
         do j=jsc,jec
            do i=isc,iec
               atm(n)%u_srf(i,j) = atm(n)%ua(i,j,npz)
               atm(n)%v_srf(i,j) = atm(n)%va(i,j,npz)
            enddo
         enddo
         atm(n)%flagstruct%srf_init = .true.
    endif

    end do   ! n_tile

  end subroutine fv_restart
  ! </SUBROUTINE> NAME="fv_restart"


  subroutine fill_nested_grid_topo_halo(Atm, proc_in)

    type(fv_atmos_type), intent(INOUT) :: Atm
    logical, intent(IN), OPTIONAL :: proc_in
    integer :: isd, ied, jsd, jed

    if (.not. atm%neststruct%nested) return

    call mpp_get_data_domain( atm%parent_grid%domain, &
         isd, ied, jsd, jed)

    !This is 2D and doesn't need remapping
    if (is_master()) print*, '  FILLING NESTED GRID HALO WITH INTERPOLATED TERRAIN'
    call nested_grid_bc(atm%phis, atm%parent_grid%phis, global_nest_domain, &
         atm%neststruct%ind_h, atm%neststruct%wt_h, 0, 0, &
         atm%npx, atm%npy, atm%bd, isd, ied, jsd, jed, proc_in=proc_in, nest_level=atm%grid_number-1)
    
  end subroutine fill_nested_grid_topo_halo

  subroutine fill_nested_grid_topo(Atm, proc_in) 
    type(fv_atmos_type), intent(INOUT) :: Atm
    logical, intent(IN), OPTIONAL :: proc_in
    real, allocatable :: g_dat(:,:,:)
    integer :: p, sending_proc
    integer :: isd_p, ied_p, jsd_p, jed_p
    integer :: isg, ieg, jsg,jeg

    logical :: process

    process = .true.
    if (present(proc_in)) then
       process = proc_in
    else
       process = .true.
    endif

!!$    if (.not. Atm%neststruct%nested) return

    call mpp_get_global_domain( atm%parent_grid%domain, &
         isg, ieg, jsg, jeg)
    call mpp_get_data_domain( atm%parent_grid%domain, &
         isd_p,  ied_p,  jsd_p,  jed_p  )

    allocate(g_dat( isg:ieg, jsg:jeg, 1) )
    call timing_on('COMM_TOTAL')

    !!! FIXME: For whatever reason this code CRASHES if the lower-left corner
    !!!        of the nested grid lies within the first PE of a grid tile.

    if (is_master() .and. .not. atm%flagstruct%external_ic ) print*, ' FILLING NESTED GRID INTERIOR WITH INTERPOLATED TERRAIN'

    sending_proc = (atm%parent_grid%pelist(1)) + &
         (atm%neststruct%parent_tile-tile_fine(atm%parent_grid%grid_number)+atm%parent_grid%flagstruct%ntiles-1)*atm%parent_grid%npes_per_tile
    if (atm%neststruct%parent_tile == atm%parent_grid%global_tile) then
    !if (Atm%neststruct%parent_proc .and. Atm%neststruct%parent_tile == Atm%parent_grid%global_tile) then
       call mpp_global_field( &
            atm%parent_grid%domain, &
            atm%parent_grid%phis(isd_p:ied_p,jsd_p:jed_p), g_dat(isg:,jsg:,1), position=center)
       if (mpp_pe() == sending_proc) then 
          do p=1,size(atm%pelist)
             call mpp_send(g_dat,size(g_dat),atm%pelist(p))
          enddo
       endif
    endif

    if (any(atm%pelist == mpp_pe())) then
       call mpp_recv(g_dat, size(g_dat), sending_proc)
    endif

    call timing_off('COMM_TOTAL')
    if (process) call fill_nested_grid(atm%phis, g_dat(isg:,jsg:,1), &
         atm%neststruct%ind_h, atm%neststruct%wt_h, &
         0, 0,  isg, ieg, jsg, jeg, atm%bd)

    call mpp_sync_self

    deallocate(g_dat)


  end subroutine fill_nested_grid_topo

  !This will still probably be needed for moving nests
  !NOTE: this has NOT been maintained and so %global_tile is now meaningless if not referring to data on the current PE
  !      needs to be re-coded to follow method in fill_nested_grid_Topo
  subroutine fill_nested_grid_data(Atm, proc_in)

    type(fv_atmos_type), intent(INOUT) :: Atm(:) !Only intended to be one element; needed for cubed_sphere_terrain
    logical, intent(IN), OPTIONAL :: proc_in
    real, allocatable :: g_dat(:,:,:), pt_coarse(:,:,:)
    integer :: i,j,k,nq, sphum, ncnst, istart, iend, npz
    integer :: isc, iec, jsc, jec, isd, ied, jsd, jed
    integer :: isd_p, ied_p, jsd_p, jed_p, isc_p, iec_p, jsc_p, jec_p
    integer :: isg, ieg, jsg,jeg, npx_p, npy_p
    integer :: isg_n, ieg_n, jsg_n, jeg_n, npx_n, npy_n
    real zvir, gh0, p1(2), p2(2), r, r0

    integer :: p, sending_proc, gid, n
    logical process

    call mpp_error(fatal, " FILL_NESTED_GRID_DATA not yet updated for remap BCs")

    if (present(proc_in)) then
       process = proc_in
    else
       process = .true.
    endif

    isd = atm(1)%bd%isd
    ied = atm(1)%bd%ied
    jsd = atm(1)%bd%jsd
    jed = atm(1)%bd%jed
    ncnst = atm(1)%ncnst
    isc = atm(1)%bd%isc; iec = atm(1)%bd%iec; jsc = atm(1)%bd%jsc; jec = atm(1)%bd%jec
    npz     = atm(1)%npz    
    

    gid = mpp_pe()

    sending_proc = atm(1)%parent_grid%pelist(1) + (atm(1)%neststruct%parent_tile-1)*atm(1)%parent_grid%npes_per_tile

       call mpp_get_data_domain( atm(1)%parent_grid%domain, &
            isd_p,  ied_p,  jsd_p,  jed_p  )
       call mpp_get_compute_domain( atm(1)%parent_grid%domain, &
            isc_p,  iec_p,  jsc_p,  jec_p  )
    call mpp_get_global_domain( atm(1)%parent_grid%domain, &
         isg, ieg, jsg, jeg, xsize=npx_p, ysize=npy_p)

    if (process) then 
       
       call mpp_error(note, "FILLING NESTED GRID DATA")

    else

       call mpp_error(note, "SENDING TO FILL NESTED GRID DATA")

    endif

    !delp

    allocate(g_dat( isg:ieg, jsg:jeg, npz) )

    call timing_on('COMM_TOTAL')

    !Call mpp_global_field on the procs that have the required data.
       !Then broadcast from the head PE to the receiving PEs
       if (atm(1)%neststruct%parent_proc .and. atm(1)%neststruct%parent_tile == atm(1)%parent_grid%global_tile) then
          call mpp_global_field( &
               atm(1)%parent_grid%domain, &
               atm(1)%parent_grid%delp(isd_p:ied_p,jsd_p:jed_p,:), g_dat, position=center)
          if (gid == sending_proc) then !crazy logic but what we have for now
             do p=1,size(atm(1)%pelist)
                call mpp_send(g_dat,size(g_dat),atm(1)%pelist(p))
             enddo
          endif
       endif
       if (any(atm(1)%pelist == gid)) then
          call mpp_recv(g_dat, size(g_dat), sending_proc)
       endif

    call timing_off('COMM_TOTAL')
    if (process) call fill_nested_grid(atm(1)%delp, g_dat, &
         atm(1)%neststruct%ind_h, atm(1)%neststruct%wt_h, &
         0, 0,  isg, ieg, jsg, jeg, npz, atm(1)%bd)

    call mpp_sync_self

    !tracers
    do nq=1,ncnst

       call timing_on('COMM_TOTAL')

          if (any(atm(1)%parent_grid%pelist == gid) .and. atm(1)%neststruct%parent_tile == atm(1)%parent_grid%global_tile) then
             call mpp_global_field( &
               atm(1)%parent_grid%domain, &
               atm(1)%parent_grid%q(isd_p:ied_p,jsd_p:jed_p,:,nq), g_dat, position=center)
             if (gid == sending_proc) then
                do p=1,size(atm(1)%pelist)
                   call mpp_send(g_dat,size(g_dat),atm(1)%pelist(p))
                enddo
             endif
          endif
          if (any(atm(1)%pelist == gid)) then
             call mpp_recv(g_dat, size(g_dat), sending_proc)
          endif

       call timing_off('COMM_TOTAL')
       if (process) call fill_nested_grid(atm(1)%q(isd:ied,jsd:jed,:,nq), g_dat, &
            atm(1)%neststruct%ind_h, atm(1)%neststruct%wt_h, &
            0, 0,  isg, ieg, jsg, jeg, npz, atm(1)%bd)

    call mpp_sync_self

    end do

    !Note that we do NOT fill in phis (surface geopotential), which should 
    !be computed exactly instead of being interpolated.


#ifndef SW_DYNAMICS
    !pt --- actually temperature

    call timing_on('COMM_TOTAL')

       if (any(atm(1)%parent_grid%pelist == gid) .and. atm(1)%neststruct%parent_tile == atm(1)%parent_grid%global_tile) then
             call mpp_global_field( &
               atm(1)%parent_grid%domain, &
               atm(1)%parent_grid%pt(isd_p:ied_p,jsd_p:jed_p,:), g_dat, position=center)
          if (gid == sending_proc) then
             do p=1,size(atm(1)%pelist)
                call mpp_send(g_dat,size(g_dat),atm(1)%pelist(p))
             enddo
          endif
       endif
       if (any(atm(1)%pelist == gid)) then
          call mpp_recv(g_dat, size(g_dat), sending_proc)
       endif

    call mpp_sync_self

    call timing_off('COMM_TOTAL')
    if (process) call fill_nested_grid(atm(1)%pt, g_dat, &
         atm(1)%neststruct%ind_h, atm(1)%neststruct%wt_h, &
         0, 0,  isg, ieg, jsg, jeg, npz, atm(1)%bd)


    if ( atm(1)%flagstruct%nwat > 0 ) then
       sphum = get_tracer_index(model_atmos, 'sphum')
    else
       sphum = 1
    endif
    if ( atm(1)%parent_grid%flagstruct%adiabatic .or. atm(1)%parent_grid%flagstruct%do_Held_Suarez ) then
       zvir = 0.         ! no virtual effect
    else
       zvir = rvgas/rdgas - 1.
    endif

    call timing_on('COMM_TOTAL')

       if (any(atm(1)%parent_grid%pelist == gid) .and. atm(1)%neststruct%parent_tile == atm(1)%parent_grid%global_tile) then
             call mpp_global_field( &
               atm(1)%parent_grid%domain, &
               atm(1)%parent_grid%pkz(isc_p:iec_p,jsc_p:jec_p,:), g_dat, position=center)
          if (gid == sending_proc) then
             do p=1,size(atm(1)%pelist)
                call mpp_send(g_dat,size(g_dat),atm(1)%pelist(p))
             enddo
          endif
       endif
       if (any(atm(1)%pelist == gid)) then
          call mpp_recv(g_dat, size(g_dat), sending_proc)
       endif

    call mpp_sync_self

    call timing_off('COMM_TOTAL')
    if (process) then 
       allocate(pt_coarse(isd:ied,jsd:jed,npz))
       call fill_nested_grid(pt_coarse, g_dat, &
            atm(1)%neststruct%ind_h, atm(1)%neststruct%wt_h, &
            0, 0,  isg, ieg, jsg, jeg, npz, atm(1)%bd)

       if (atm(1)%bd%is == 1) then
          do k=1,npz
             do j=atm(1)%bd%jsd,atm(1)%bd%jed
                do i=atm(1)%bd%isd,0
#ifdef MULTI_GASES
                   atm(1)%pt(i,j,k) = cp_air*atm(1)%pt(i,j,k)/pt_coarse(i,j,k)*virq(atm(1)%q(i,j,k,:))
#else
                   atm(1)%pt(i,j,k) = cp_air*atm(1)%pt(i,j,k)/pt_coarse(i,j,k)*(1.+zvir*atm(1)%q(i,j,k,sphum))
#endif
                end do
             end do
          end do
       end if

       if (atm(1)%bd%js == 1) then
          if (atm(1)%bd%is == 1) then
             istart = atm(1)%bd%is
          else
             istart = atm(1)%bd%isd
          end if
          if (atm(1)%bd%ie == atm(1)%npx-1) then
             iend = atm(1)%bd%ie
          else
             iend = atm(1)%bd%ied
          end if

          do k=1,npz
             do j=atm(1)%bd%jsd,0
                do i=istart,iend
#ifdef MULTI_GASES
                   atm(1)%pt(i,j,k) = cp_air*atm(1)%pt(i,j,k)/pt_coarse(i,j,k)*virq(atm(1)%q(i,j,k,:))
#else
                   atm(1)%pt(i,j,k) = cp_air*atm(1)%pt(i,j,k)/pt_coarse(i,j,k)*(1.+zvir*atm(1)%q(i,j,k,sphum))
#endif
                end do
             end do
          end do
       end if

       if (atm(1)%bd%ie == atm(1)%npx-1) then
          do k=1,npz
             do j=atm(1)%bd%jsd,atm(1)%bd%jed
                do i=atm(1)%npx,atm(1)%bd%ied
#ifdef MULTI_GASES
                   atm(1)%pt(i,j,k) = cp_air*atm(1)%pt(i,j,k)/pt_coarse(i,j,k)*virq(atm(1)%q(i,j,k,:))
#else
                   atm(1)%pt(i,j,k) = cp_air*atm(1)%pt(i,j,k)/pt_coarse(i,j,k)*(1.+zvir*atm(1)%q(i,j,k,sphum))
#endif
                end do
             end do
          end do
       end if

       if (atm(1)%bd%je == atm(1)%npy-1) then
          if (atm(1)%bd%is == 1) then
             istart = atm(1)%bd%is
          else
             istart = atm(1)%bd%isd
          end if
          if (atm(1)%bd%ie == atm(1)%npx-1) then
             iend = atm(1)%bd%ie
          else
             iend = atm(1)%bd%ied
          end if

          do k=1,npz
             do j=atm(1)%npy,atm(1)%bd%jed
                do i=istart,iend
#ifdef MULTI_GASES
                   atm(1)%pt(i,j,k) = cp_air*atm(1)%pt(i,j,k)/pt_coarse(i,j,k)*virq(atm(1)%q(i,j,k,:))
#else
                   atm(1)%pt(i,j,k) = cp_air*atm(1)%pt(i,j,k)/pt_coarse(i,j,k)*(1.+zvir*atm(1)%q(i,j,k,sphum))
#endif
                end do
             end do
          end do
       end if

       deallocate(pt_coarse)

    end if

    if (.not. atm(1)%flagstruct%hydrostatic) then

       !delz
       call timing_on('COMM_TOTAL')

          if (any(atm(1)%parent_grid%pelist == gid) .and. atm(1)%neststruct%parent_tile == atm(1)%parent_grid%global_tile) then
             call mpp_global_field( &
               atm(1)%parent_grid%domain, &
               atm(1)%parent_grid%delz(isd_p:ied_p,jsd_p:jed_p,:), g_dat, position=center)
             if (gid == sending_proc) then
                do p=1,size(atm(1)%pelist)
                   call mpp_send(g_dat,size(g_dat),atm(1)%pelist(p))
                enddo
             endif
          endif
          if (any(atm(1)%pelist == gid)) then
             call mpp_recv(g_dat, size(g_dat), sending_proc)
          endif

    call mpp_sync_self

       call timing_off('COMM_TOTAL')
       if (process) call fill_nested_grid(atm(1)%delz, g_dat, &
            atm(1)%neststruct%ind_h, atm(1)%neststruct%wt_h, &
            0, 0,  isg, ieg, jsg, jeg, npz, atm(1)%bd)

       !w

       call timing_on('COMM_TOTAL')

          if (any(atm(1)%parent_grid%pelist == gid) .and. atm(1)%neststruct%parent_tile == atm(1)%parent_grid%global_tile) then
             call mpp_global_field( &
               atm(1)%parent_grid%domain, &
               atm(1)%parent_grid%w(isd_p:ied_p,jsd_p:jed_p,:), g_dat, position=center)
             if (gid == sending_proc) then
                do p=1,size(atm(1)%pelist)
                   call mpp_send(g_dat,size(g_dat),atm(1)%pelist(p))
                enddo
             endif
          endif
          if (any(atm(1)%pelist == gid)) then
             call mpp_recv(g_dat, size(g_dat), sending_proc)
          endif

    call mpp_sync_self

       call timing_off('COMM_TOTAL')
       if (process) call fill_nested_grid(atm(1)%w, g_dat, &
            atm(1)%neststruct%ind_h, atm(1)%neststruct%wt_h, &
            0, 0,  isg, ieg, jsg, jeg, npz, atm(1)%bd)
       !

    end if

#endif
    deallocate(g_dat) 

    !u

    allocate(g_dat( isg:ieg, jsg:jeg+1, npz) )
    g_dat = 1.e25

    call timing_on('COMM_TOTAL')

       if (any(atm(1)%parent_grid%pelist == gid) .and. atm(1)%neststruct%parent_tile == atm(1)%parent_grid%global_tile) then
             call mpp_global_field( &
               atm(1)%parent_grid%domain, &
               atm(1)%parent_grid%u(isd_p:ied_p,jsd_p:jed_p+1,:), g_dat, position=north)
          if (gid == sending_proc) then
             do p=1,size(atm(1)%pelist)
                call mpp_send(g_dat,size(g_dat),atm(1)%pelist(p))
             enddo
          endif
       endif
       if (any(atm(1)%pelist == gid)) then
          call mpp_recv(g_dat, size(g_dat), sending_proc)
       endif

    call mpp_sync_self

    call timing_off('COMM_TOTAL')
    call mpp_sync_self
    if (process) call fill_nested_grid(atm(1)%u, g_dat, &
         atm(1)%neststruct%ind_u, atm(1)%neststruct%wt_u, &
         0, 1,  isg, ieg, jsg, jeg, npz, atm(1)%bd)
    deallocate(g_dat)

    !v

    allocate(g_dat( isg:ieg+1, jsg:jeg, npz) )
    g_dat = 1.e25

    call timing_on('COMM_TOTAL')

       if (any(atm(1)%parent_grid%pelist == gid) .and. atm(1)%neststruct%parent_tile == atm(1)%parent_grid%global_tile) then
             call mpp_global_field( &
               atm(1)%parent_grid%domain, &
               atm(1)%parent_grid%v(isd_p:ied_p+1,jsd_p:jed_p,:), g_dat, position=east)
          if (gid == sending_proc) then
             do p=1,size(atm(1)%pelist)
                call mpp_send(g_dat,size(g_dat),atm(1)%pelist(p))
             enddo
          endif
       endif
       if (any(atm(1)%pelist == gid)) then
          call mpp_recv(g_dat, size(g_dat), sending_proc)
       endif

    call mpp_sync_self
                                      call timing_off('COMM_TOTAL')

    if (process) call fill_nested_grid(atm(1)%v, g_dat, &
         atm(1)%neststruct%ind_v, atm(1)%neststruct%wt_v, &
         1, 0,  isg, ieg, jsg, jeg, npz, atm(1)%bd)

    deallocate(g_dat)

  end subroutine fill_nested_grid_data

 subroutine twoway_topo_update(Atm, proc_in)
    type(fv_atmos_type), intent(INOUT) :: Atm  
    logical, intent(IN), OPTIONAL :: proc_in
    real, allocatable :: g_dat(:,:,:), pt_coarse(:,:,:)
    integer :: i,j,k,nq, sphum, ncnst, istart, iend, npz
    integer :: isc, iec, jsc, jec, isd, ied, jsd, jed
    integer :: isd_p, ied_p, jsd_p, jed_p, isc_p, iec_p, jsc_p, jec_p
    integer :: isg, ieg, jsg,jeg, npx_p, npy_p
    integer :: isg_n, ieg_n, jsg_n, jeg_n, npx_n, npy_n
    real zvir

    integer :: p , sending_proc
    logical :: process

    if (present(proc_in)) then
       process = proc_in
    else
       process = .true.
    endif

    isd = atm%bd%isd
    ied = atm%bd%ied
    jsd = atm%bd%jsd
    jed = atm%bd%jed
    ncnst = atm%ncnst
    isc = atm%bd%isc; iec = atm%bd%iec; jsc = atm%bd%jsc; jec = atm%bd%jec
    npz     = atm%npz    
    
    isd_p = atm%parent_grid%bd%isd
    ied_p = atm%parent_grid%bd%ied
    jsd_p = atm%parent_grid%bd%jsd
    jed_p = atm%parent_grid%bd%jed
    isc_p = atm%parent_grid%bd%isc
    iec_p = atm%parent_grid%bd%iec
    jsc_p = atm%parent_grid%bd%jsc
    jec_p = atm%parent_grid%bd%jec
    sending_proc = atm%parent_grid%pelist(1) + (atm%neststruct%parent_tile-1)*atm%parent_grid%npes_per_tile

    call mpp_get_global_domain( atm%parent_grid%domain, &
         isg, ieg, jsg, jeg, xsize=npx_p, ysize=npy_p)


    !NOW: what we do is to update the nested-grid terrain to the coarse grid,
    !to ensure consistency between the two grids.
    if ( process ) call mpp_update_domains(atm%phis, atm%domain, complete=.true.)
    if (atm%neststruct%twowaynest) then
       if (any(atm%parent_grid%pelist == mpp_pe()) .or. atm%neststruct%child_proc) then
          call update_coarse_grid(atm%parent_grid%phis, &
               atm%phis, global_nest_domain, &
               atm%gridstruct%dx, atm%gridstruct%dy, atm%gridstruct%area, &
               atm%bd, isd_p, ied_p, jsd_p, jed_p, isd, ied, jsd, jed, &
               atm%neststruct%isu, atm%neststruct%ieu, atm%neststruct%jsu, atm%neststruct%jeu, &
               atm%npx, atm%npy, 0, 0, &
               atm%neststruct%refinement, atm%neststruct%nestupdate, 0, 0, &
               atm%neststruct%parent_proc, atm%neststruct%child_proc, atm%parent_grid, atm%grid_number-1)
          atm%parent_grid%neststruct%parent_of_twoway = .true.
          !NOTE: mpp_update_nest_coarse (and by extension, update_coarse_grid) does **NOT** pass data
          !allowing a two-way update into the halo of the coarse grid. It only passes data so that the INTERIOR
          ! can have the two-way update. Thus, on the nest's cold start, if this update_domains call is not done,
          ! the coarse grid will have the wrong topography in the halo, which will CHANGE when a restart is done!!
          if (atm%neststruct%parent_proc) call mpp_update_domains(atm%parent_grid%phis, atm%parent_grid%domain)
       end if

    end if


#ifdef SW_DYNAMICS
!!$    !ps: first level only
!!$    !This is only valid for shallow-water simulations
!!$    if (process) then
!!$    do j=jsd,jed
!!$       do i=isd,ied
!!$
!!$          Atm%ps(i,j) = Atm%delp(i,j,1)/grav
!!$
!!$       end do
!!$    end do
!!$    endif
#else
    !Reset p_var after updating topography
    if (process) call p_var(npz, isc, iec, jsc, jec, atm%ptop, ptop_min, atm%delp, &
         atm%delz, atm%pt, atm%ps,   &
         atm%pe, atm%peln, atm%pk, atm%pkz, kappa, atm%q, &
         atm%ng, ncnst, atm%gridstruct%area_64, atm%flagstruct%dry_mass, .false., atm%flagstruct%mountain, &
         atm%flagstruct%moist_phys, .true., atm%flagstruct%nwat, atm%domain, atm%flagstruct%adiabatic)
#endif

 

  end subroutine twoway_topo_update

  subroutine fv_write_restart(Atm, timestamp)
    type(fv_atmos_type), intent(inout) :: Atm
    character(len=*),    intent(in)    :: timestamp

    call fv_io_write_restart(atm, timestamp)
    if (atm%neststruct%nested) then
       call fv_io_write_bcs(atm)
    endif

  end subroutine fv_write_restart

  subroutine fv_restart_end(Atm, restart_endfcst)
    type(fv_atmos_type), intent(inout) :: Atm
    logical, intent(in) :: restart_endfcst

    integer :: isc, iec, jsc, jec
    integer :: iq, ncnst, ntprog, ntdiag
    integer :: isd, ied, jsd, jed, npz
    integer :: unit
    integer :: file_unit
    integer, allocatable :: pelist(:)
    character(len=128):: tracer_name
    character(len=3):: gn

    call mpp_set_current_pelist(atm%pelist)

    isc = atm%bd%isc; iec = atm%bd%iec; jsc = atm%bd%jsc; jec = atm%bd%jec

    isd = atm%bd%isd
    ied = atm%bd%ied
    jsd = atm%bd%jsd
    jed = atm%bd%jed
    npz = atm%npz
    ncnst = atm%ncnst
    ntprog = size(atm%q,4)
    ntdiag = size(atm%qdiag,4)

    if (atm%grid_number > 1) then
       write(gn,'(A2, I1)') " g", atm%grid_number
    else
       gn = ''
    end if

    unit = stdout()
    write(unit,*)
    write(unit,*) 'fv_restart_end u   ', trim(gn),' = ', mpp_chksum(atm%u(isc:iec,jsc:jec,:))
    write(unit,*) 'fv_restart_end v   ', trim(gn),' = ', mpp_chksum(atm%v(isc:iec,jsc:jec,:))
    if ( .not. atm%flagstruct%hydrostatic )    &
         write(unit,*) 'fv_restart_end w   ', trim(gn),' = ', mpp_chksum(atm%w(isc:iec,jsc:jec,:))
    write(unit,*) 'fv_restart_end delp', trim(gn),' = ', mpp_chksum(atm%delp(isc:iec,jsc:jec,:))
    write(unit,*) 'fv_restart_end phis', trim(gn),' = ', mpp_chksum(atm%phis(isc:iec,jsc:jec))
#ifndef SW_DYNAMICS
    write(unit,*) 'fv_restart_end pt  ', trim(gn),' = ', mpp_chksum(atm%pt(isc:iec,jsc:jec,:))
    if (ntprog>0) &
         write(unit,*) 'fv_restart_end q(prog) nq  ', trim(gn),' =',ntprog, mpp_chksum(atm%q(isc:iec,jsc:jec,:,:))
    if (ntdiag>0) &
         write(unit,*) 'fv_restart_end q(diag) nq  ', trim(gn),' =',ntdiag, mpp_chksum(atm%qdiag(isc:iec,jsc:jec,:,:))
    do iq=1,min(17, ntprog)     ! Check up to 17 tracers
       call get_tracer_names(model_atmos, iq, tracer_name)
       write(unit,*) 'fv_restart_end '//trim(tracer_name)// trim(gn),' = ', mpp_chksum(atm%q(isc:iec,jsc:jec,:,iq))
    enddo

    !---------------
    ! Check Min/Max:
    !---------------
    !     call prt_maxmin('ZS', Atm%phis, isc, iec, jsc, jec, Atm%ng, 1, 1./grav)
    call pmaxmn_g('ZS', atm%phis, isc, iec, jsc, jec, 1, 1./grav, atm%gridstruct%area_64, atm%domain)
    call pmaxmn_g('PS ', atm%ps,   isc, iec, jsc, jec, 1, 0.01   , atm%gridstruct%area_64, atm%domain)
    call prt_maxmin('PS*', atm%ps, isc, iec, jsc, jec, atm%ng, 1, 0.01)
    call prt_maxmin('U ', atm%u(isd:ied,jsd:jed,1:npz), isc, iec, jsc, jec, atm%ng, npz, 1.)
    call prt_maxmin('V ', atm%v(isd:ied,jsd:jed,1:npz), isc, iec, jsc, jec, atm%ng, npz, 1.)
    if ( .not. atm%flagstruct%hydrostatic )    &
         call prt_maxmin('W ', atm%w , isc, iec, jsc, jec, atm%ng, npz, 1.)
    call prt_maxmin('T ', atm%pt, isc, iec, jsc, jec, atm%ng, npz, 1.)
    do iq=1, ntprog
       call get_tracer_names ( model_atmos, iq, tracer_name )
       call pmaxmn_g(trim(tracer_name), atm%q(isd:ied,jsd:jed,1:npz,iq:iq), isc, iec, jsc, jec, npz, &
            1., atm%gridstruct%area_64, atm%domain)
    enddo
    ! Write4 energy correction term
#endif

   if ( restart_endfcst ) then
    call fv_io_write_restart(atm)
    if (atm%neststruct%nested) call fv_io_write_bcs(atm)
   endif
   if(atm%flagstruct%write_restart_with_bcs)then
       call write_full_fields(atm)
   endif

 module_is_initialized = .false.

#ifdef EFLUX_OUT
 if( is_master() ) then
    write(*,*) steps, 'Mean equivalent Heat flux for this integration period=',atm(1)%idiag%efx_sum/real(max(1,Atm(1)%idiag%steps)), &
         'Mean nesting-related flux for this integration period=',Atm(1)%idiag%efx_sum_nest/real(max(1,Atm(1)%idiag%steps)), &
         'Mean mountain torque=',Atm(1)%idiag%mtq_sum/real(max(1,atm(1)%idiag%steps))
    file_unit = get_unit()
    open (unit=file_unit, file='e_flux.data', form='unformatted',status='unknown', access='sequential')
    do n=1,steps
       write(file_unit) atm(1)%idiag%efx(n)
       write(file_unit) atm(1)%idiag%mtq(n)    ! time series global mountain torque
       !write(file_unit) Atm(1)%idiag%efx_nest(n)  
    enddo
    close(unit=file_unit)
 endif
#endif

  end subroutine fv_restart_end


subroutine pmaxmn_g(qname, q, is, ie, js, je, km, fac, area, domain)
      character(len=*), intent(in)::  qname
      integer, intent(in):: is, ie, js, je
      integer, intent(in):: km
      real, intent(in)::    q(is-3:ie+3, js-3:je+3, km)
      real, intent(in)::    fac
      real(kind=R_GRID), intent(IN)::    area(is-3:ie+3, js-3:je+3)
      type(domain2d), intent(INOUT) :: domain
!
      real qmin, qmax, gmean
      integer i,j,k

      qmin = q(is,js,1)
      qmax = qmin

      do k=1,km
      do j=js,je
         do i=is,ie
            !if ( (q(i,j,k) >= 1e30) .eqv. (q(i,j,k) < 1e30) ) then !NAN checking
            !   print*, ' NAN found for ', qname, mpp_pe(), i,j,k
            !else
            if( q(i,j,k) < qmin) then
                qmin = q(i,j,k)
            elseif( q(i,j,k) > qmax ) then
                qmax = q(i,j,k)
            endif
          enddo
      enddo
      enddo

      call mp_reduce_min(qmin)
      call mp_reduce_max(qmax)

      gmean = g_sum(domain, q(is:ie,js:je,km), is, ie, js, je, 3, area, 1, .true.) 
      if(is_master()) write(6,*) qname, qmax*fac, qmin*fac, gmean*fac

end subroutine pmaxmn_g
end module fv_restart_mod