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_register_restart_bcs_nh, 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_get_compute_domain, mpp_get_data_domain, mpp_get_global_domain
  use mpp_domains_mod,     only: mpp_update_domains, domain2d, dgrid_ne
  use mpp_domains_mod,     only: center, corner, north, east,  mpp_get_c2f_index, west, south
  use mpp_domains_mod,     only: mpp_global_field
  use mpp_mod,             only: mpp_chksum, stdout, mpp_error, fatal, note
  use mpp_mod,             only: get_unit, mpp_sum
  use mpp_mod,             only: mpp_get_current_pelist, mpp_set_current_pelist
  use mpp_mod,             only: mpp_send, mpp_recv, mpp_sync_self, mpp_npes, mpp_pe, mpp_sync
  use test_cases_mod,      only: test_case, alpha, init_case, init_double_periodic, init_latlon
  use fv_mp_mod,           only: is_master, switch_current_atm, mp_reduce_min, mp_reduce_max
  use fv_surf_map_mod,     only: sgh_g, oro_g
  use fv_surf_map_mod,     only: del2_cubed_sphere, del4_cubed_sphere
  use tracer_manager_mod,  only: get_tracer_names
  use tracer_manager_mod,  only: get_tracer_index
  use field_manager_mod,   only: model_atmos
  use external_ic_mod,     only: get_external_ic, get_cubed_sphere_terrain
  use fv_eta_mod,          only: compute_dz_var, compute_dz_l32, set_hybrid_z
  use boundary_mod,        only: fill_nested_grid, nested_grid_bc, update_coarse_grid 
  use field_manager_mod,   only: model_atmos
  use fv_timing_mod,       only: timing_on, timing_off
  use fms_mod,             only: file_exist
  use fv_treat_da_inc_mod, only: read_da_inc
#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, setup_nested_boundary_halo
  public :: d2c_setup, d2a_setup

  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, grids_on_this_pe)
    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
    logical, intent(INOUT) :: grids_on_this_pe(:)


    integer :: i, j, k, n, ntileme, nt, iq
    integer :: isc, iec, jsc, jec, npz, npz_rst, ncnst, ntprog, ntdiag
    integer :: isd, ied, jsd, jed
    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
    logical :: hybrid
    logical :: cold_start_grids(size(atm))
    character(len=128):: tname, errstring, fname, tracer_name
    character(len=120):: fname_ne, fname_sw
    character(len=3) :: gn

    integer :: npts
    real    :: sumpertn

    rgrav = 1. / grav

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

    ntileme = size(atm(:))

    cold_start_grids(:) = cold_start
    do n = 1, ntileme

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

       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 (atm(n)%flagstruct%external_ic) then
             if (is_master()) print*, 'External IC set on grid', atm(n)%grid_number, ', re-initializing grid'
             cold_start_grids(n) = .true.
             atm(n)%flagstruct%warm_start = .false. !resetting warm_start flag to avoid FATAL error below
          else
             if (is_master()) print*, 'Searching for nested grid restart file ', trim(fname)
             cold_start_grids(n) = .not. file_exist(fname, atm(n)%domain)
             atm(n)%flagstruct%warm_start = file_exist(fname, atm(n)%domain)!resetting warm_start flag to avoid FATAL error below
          endif
       endif

       if (.not. grids_on_this_pe(n)) then
          
          !Even if this grid is not on this PE, if it has child grids we must send
          !along the data that is needed. 
          !This is a VERY complicated bit of code that attempts to follow the entire decision tree
          ! of the initialization without doing anything. This could very much be cleaned up.

          if (atm(n)%neststruct%nested) then
             if (cold_start_grids(n)) then
                if (atm(n)%parent_grid%flagstruct%n_zs_filter > 0) call fill_nested_grid_topo_halo(atm(n), .false.)
                if (atm(n)%flagstruct%nggps_ic) then
                   call fill_nested_grid_topo(atm(n), .false.)
                   call fill_nested_grid_topo_halo(atm(n), .false.)
                   call nested_grid_bc(atm(n)%ps, atm(n)%parent_grid%ps, atm(n)%neststruct%nest_domain, &
                        atm(n)%neststruct%ind_h, atm(n)%neststruct%wt_h, 0, 0, &
                        atm(n)%npx, atm(n)%npy,atm(n)%bd, isg, ieg, jsg, jeg, proc_in=.false.)         
                   call setup_nested_boundary_halo(atm(n),.false.) 
                else
                   call fill_nested_grid_topo(atm(n), .false.)
                   call setup_nested_boundary_halo(atm(n),.false.) 
                   if ( atm(n)%flagstruct%external_ic .and. grid_type < 4 ) call fill_nested_grid_data(atm(n:n), .false.)
                endif
             else
                if (is_master()) print*, 'Searching for nested grid BC files ', trim(fname_ne), ' ', trim(fname_sw)

                !!!! PROBLEM: file_exist doesn't know to look for fv_BC_ne.res.nest02.nc instead of fv_BC_ne.res.nc on coarse grid
                if (file_exist(fname_ne, atm(n)%domain) .and. file_exist(fname_sw, atm(n)%domain)) then
                else
                   if ( is_master() ) write(*,*) 'BC files not found, re-generating nested grid boundary conditions'
                   call fill_nested_grid_topo_halo(atm(n), .false.)
                   call setup_nested_boundary_halo(atm(n), .false.)
                   atm(n)%neststruct%first_step = .true.                   
                endif
             end if

             if (.not. atm(n)%flagstruct%hydrostatic .and. atm(n)%flagstruct%make_nh .and. &
                  (.not. atm(n)%flagstruct%nggps_ic .and. .not. atm(n)%flagstruct%ecmwf_ic) ) then
                call nested_grid_bc(atm(n)%delz, atm(n)%parent_grid%delz, atm(n)%neststruct%nest_domain, &
                     atm(n)%neststruct%ind_h, atm(n)%neststruct%wt_h, 0, 0, &
                     atm(n)%npx, atm(n)%npy, atm(n)%npz, atm(n)%bd, isg, ieg, jsg, jeg, proc_in=.false.)
                call nested_grid_bc(atm(n)%w, atm(n)%parent_grid%w, atm(n)%neststruct%nest_domain, &
                     atm(n)%neststruct%ind_h, atm(n)%neststruct%wt_h, 0, 0, &
                     atm(n)%npx, atm(n)%npy, atm(n)%npz, atm(n)%bd, isg, ieg, jsg, jeg, proc_in=.false.)
             endif

          endif

          cycle

       endif
       !This call still appears to be necessary to get isd, etc. correct
       call switch_current_atm(atm(n))

    !--- call fv_io_register_restart to register restart field to be written out in fv_io_write_restart
    call fv_io_register_restart(atm(n)%domain,atm(n:n))
    if (atm(n)%neststruct%nested) call fv_io_register_restart_bcs(atm(n))
    if( .not.cold_start_grids(n) .and. (.not. atm(n)%flagstruct%external_ic) ) 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:n), atm(n)%domain)
                if( is_master() ) write(*,*) 'Back from read_da_inc',atm(n)%pt(i,j,k)
             endif
!            ====== end PJP added DA functionailty======
        endif
    endif

!---------------------------------------------------------------------------------------------
! Read, interpolate (latlon to cubed), then remap vertically with terrain adjustment if needed
!---------------------------------------------------------------------------------------------
    if (atm(n)%neststruct%nested) then
          if (cold_start_grids(n)) call fill_nested_grid_topo(atm(n), .true.)
          !if (cold_start_grids(n) .and. .not. Atm(n)%flagstruct%nggps_ic) call fill_nested_grid_topo(Atm(n), .true.)
       if (cold_start_grids(n)) then
          if (atm(n)%parent_grid%flagstruct%n_zs_filter > 0 .or. atm(n)%flagstruct%nggps_ic) call fill_nested_grid_topo_halo(atm(n), .true.)
       end if
       if (atm(n)%flagstruct%external_ic .and. atm(n)%flagstruct%nggps_ic) then
          !Fill nested grid halo with ps
          call nested_grid_bc(atm(n)%ps, atm(n)%parent_grid%ps, atm(n)%neststruct%nest_domain, &
               atm(n)%neststruct%ind_h, atm(n)%neststruct%wt_h, 0, 0, &
               atm(n)%npx, atm(n)%npy,atm(n)%bd, isg, ieg, jsg, jeg, proc_in=.true.)         
       endif
    endif
    if ( atm(n)%flagstruct%external_ic ) then
         if( is_master() ) write(*,*) 'Calling get_external_ic'
         call get_external_ic(atm(n:n), atm(n)%domain, cold_start_grids(n)) 
         if( is_master() ) write(*,*) 'IC generated from the specified external source'
    endif

    seconds = 0; days = 0   ! Restart needs to be modified to record seconds and days.

! Notes by Jeff D.
  ! This logic doesn't work very well.
  ! Shouldn't have read for all tiles then loop over tiles

       isd = atm(n)%bd%isd
       ied = atm(n)%bd%ied
       jsd = atm(n)%bd%jsd
       jed = atm(n)%bd%jed
       ncnst = atm(n)%ncnst
       if( is_master() ) write(*,*) 'in fv_restart ncnst=', ncnst
       isc = atm(n)%bd%isc; iec = atm(n)%bd%iec; jsc = atm(n)%bd%jsc; jec = atm(n)%bd%jec

    ! Init model data
       if(.not.cold_start_grids(n))then
          atm(n)%neststruct%first_step = .false.
          if (atm(n)%neststruct%nested) then
             if ( atm(n)%flagstruct%npz_rst /= 0 .and. atm(n)%flagstruct%npz_rst /= atm(n)%npz ) then
                call setup_nested_boundary_halo(atm(n)) 
             else
                !If BC file is found, then read them in. Otherwise we need to initialize the BCs.
                if (is_master()) print*, 'Searching for nested grid BC files ', trim(fname_ne), ' ', trim(fname_sw)
                if (file_exist(fname_ne, atm(n)%domain) .and. file_exist(fname_sw, atm(n)%domain)) then
                   call fv_io_read_bcs(atm(n))
                else
                   if ( is_master() ) write(*,*) 'BC files not found, re-generating nested grid boundary conditions'
                   call fill_nested_grid_topo_halo(atm(n), .true.)
                   call setup_nested_boundary_halo(atm(n), .true.)
                   atm(n)%neststruct%first_step = .true.
                endif
                !Following line to make sure u and v are consistent across processor subdomains
                call mpp_update_domains(atm(n)%u, atm(n)%v, atm(n)%domain, gridtype=dgrid_ne, complete=.true.)
             endif
          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
                   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)%neststruct%nested, atm(n)%domain, atm(n)%bd, atm(n)%flagstruct%regional)
                   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)%neststruct%nested, &
                                          atm(n)%domain, atm(n)%bd, atm(n)%flagstruct%regional)
                 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

#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(atm(n)%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(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
     else
       if ( atm(n)%flagstruct%warm_start ) then
         call mpp_error(fatal, 'FV restart files not found; set warm_start = .F. if cold_start is desired.')
      endif
! Cold start
       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, atm(n)%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, 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, atm(n)%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, atm(n)%bd)
            if( is_master() ) write(*,*) 'Doubly Periodic IC generated'
         elseif (grid_type == 5 .or. grid_type == 6) then
            call init_latlon(atm(n)%u,atm(n)%v,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)%npx, atm(n)%npy, atm(n)%npz, atm(n)%ng, ncnst, &
                             atm(n)%flagstruct%ndims, atm(n)%flagstruct%ntiles, &
                             atm(n)%flagstruct%dry_mass, &
                             atm(n)%flagstruct%mountain,       &
                             atm(n)%flagstruct%moist_phys, hybrid, atm(n)%delz, &
                             atm(n)%ze0, atm(n)%domain, atm(n)%tile)
         endif

         !Turn this off on the nested grid if you are just interpolating topography from the coarse grid!
        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


          !Set up nested grids
          !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
          ! Only fill nested-grid data if external_ic is called for the cubed-sphere grid
          if (atm(n)%neststruct%nested) then
             call setup_nested_boundary_halo(atm(n), .true.) 
             if (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

       if ( (.not.atm(n)%flagstruct%hydrostatic) .and. atm(n)%flagstruct%make_nh .and. atm(n)%neststruct%nested) then
          call nested_grid_bc(atm(n)%delz, atm(n)%parent_grid%delz, atm(n)%neststruct%nest_domain, &
               atm(n)%neststruct%ind_h, atm(n)%neststruct%wt_h, 0, 0, &
               atm(n)%npx, atm(n)%npy, npz, atm(n)%bd, isg, ieg, jsg, jeg, proc_in=.true.)
          call nested_grid_bc(atm(n)%w, atm(n)%parent_grid%w, atm(n)%neststruct%nest_domain, &
               atm(n)%neststruct%ind_h, atm(n)%neststruct%wt_h, 0, 0, &
               atm(n)%npx, atm(n)%npy, npz, atm(n)%bd, isg, ieg, jsg, jeg, proc_in=.true.)
          call fv_io_register_restart_bcs_nh(atm(n)) !needed to register nested-grid BCs not registered earlier
       endif

  end do


    do n = ntileme,1,-1
       if (atm(n)%neststruct%nested .and. atm(n)%flagstruct%external_ic .and. &
            atm(n)%flagstruct%grid_type < 4 .and. cold_start_grids(n)) then
          call fill_nested_grid_data_end(atm(n), grids_on_this_pe(n))
       endif
    end do

    do n = 1, ntileme
       if (.not. grids_on_this_pe(n)) cycle

       isd = atm(n)%bd%isd
       ied = atm(n)%bd%ied
       jsd = atm(n)%bd%jsd
       jed = atm(n)%bd%jed
       ncnst = atm(n)%ncnst
       ntprog = size(atm(n)%q,4)
       ntdiag = size(atm(n)%qdiag,4)
       isc = atm(n)%bd%isc; iec = atm(n)%bd%iec; jsc = atm(n)%bd%jsc; jec = atm(n)%bd%jec

!---------------------------------------------------------------------------------------------
! 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(atm(n)%npz) )
          if( atm(n)%npz==32 ) then
              call compute_dz_l32(atm(n)%npz, ztop, dz1)
          else
              ztop = 45.e3
              call compute_dz_var(atm(n)%npz, ztop, dz1)
          endif
          call set_hybrid_z(isc, iec, jsc, jec, atm(n)%ng, atm(n)%npz, ztop, dz1, rgrav,  &
                            atm(n)%phis, atm(n)%ze0)
          deallocate ( dz1 )
!         call prt_maxmin('ZE0', Atm(n)%ze0,  isc, iec, jsc, jec, 0, Atm(n)%npz, 1.E-3)
!         call prt_maxmin('DZ0', Atm(n)%delz, isc, iec, jsc, jec, 0, Atm(n)%npz, 1.   )
       endif
!      call make_eta_level(Atm(n)%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,atm(n)%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)%grid_number > 1) then
         write(gn,'(A2, I1)') " g", atm(n)%grid_number
      else
         gn = ''
      end if

      unit = stdout()
      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, atm(n)%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:atm(n)%npz,i:i), isc, iec, jsc, jec, atm(n)%npz, &
                        1., atm(n)%gridstruct%area_64, atm(n)%domain)
      enddo
#endif
      call prt_maxmin('U ', atm(n)%u(isc:iec,jsc:jec,1:atm(n)%npz), isc, iec, jsc, jec, 0, atm(n)%npz, 1.)
      call prt_maxmin('V ', atm(n)%v(isc:iec,jsc:jec,1:atm(n)%npz), isc, iec, jsc, jec, 0, atm(n)%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.
         if ( .not.atm(n)%flagstruct%hybrid_z ) then
            call mpp_error(note, "  Initializing delz from hydrostatic state")
             do k=1,atm(n)%npz
                do j=jsc,jec
                   do i=isc,iec
                      atm(n)%delz(i,j,k) = (rdgas*rgrav)*atm(n)%pt(i,j,k)*(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, atm(n)%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, atm(n)%npz, 1, &              
              atm(n)%gridstruct%grid_type, atm(n)%domain, &
              atm(n)%gridstruct%nested, 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,atm(n)%npz)
               atm(n)%v_srf(i,j) = atm(n)%va(i,j,atm(n)%npz)
            enddo
         enddo
         atm(n)%flagstruct%srf_init = .true.
    endif

    end do   ! n_tile

  end subroutine fv_restart

  subroutine setup_nested_boundary_halo(Atm, proc_in)

    !This routine is now taking the "easy way out" with regards
    ! to pt (virtual potential temperature), q_con, and cappa;
    ! their halo values are now set up when the BCs are set up
    ! in fv_dynamics

    type(fv_atmos_type), intent(INOUT) :: atm
    logical, INTENT(IN), OPTIONAL :: proc_in
    real, allocatable :: g_dat(:,:,:), g_dat2(:,:,:)
    real, allocatable :: pt_coarse(:,:,:)
    integer i,j,k,nq, sphum, ncnst, istart, iend, npz, nwat
    integer isc, iec, jsc, jec, isd, ied, jsd, jed, is, ie, js, je
    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 zvir
    logical process
    integer :: liq_wat, ice_wat, rainwat, snowwat, graupel
    real :: qv, dp1, q_liq, q_sol, q_con, cvm, cappa, dp, pt, dz, pkz, rdg

    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
    is    = atm%bd%is ; ie  = atm%bd%ie ; js  = atm%bd%js ; je  = atm%bd%je
    npz   = atm%npz    
    nwat  = atm%flagstruct%nwat

   if (nwat >= 3) then
      liq_wat = get_tracer_index(model_atmos, 'liq_wat')
      ice_wat = get_tracer_index(model_atmos, 'ice_wat')
   endif
   if ( nwat== 5) then
      rainwat = get_tracer_index(model_atmos, 'rainwat')
      snowwat = get_tracer_index(model_atmos, 'snowwat')
   elseif (nwat == 6) then
      rainwat = get_tracer_index(model_atmos, 'rainwat')
      snowwat = get_tracer_index(model_atmos, 'snowwat')
      graupel = get_tracer_index(model_atmos, 'graupel')
   endif

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

    call nested_grid_bc(atm%delp, atm%parent_grid%delp, atm%neststruct%nest_domain, &
         atm%neststruct%ind_h, atm%neststruct%wt_h, 0, 0, &
         atm%npx, atm%npy, npz, atm%bd, isg, ieg, jsg, jeg, proc_in=process)
    do nq=1,ncnst
       call nested_grid_bc(atm%q(:,:,:,nq), &
            atm%parent_grid%q(:,:,:,nq), atm%neststruct%nest_domain, &
            atm%neststruct%ind_h, atm%neststruct%wt_h, 0, 0, &
            atm%npx, atm%npy, npz, atm%bd, isg, ieg, jsg, jeg, proc_in=process)
    end do

    if (process) then
       if (is_master()) print*, 'FILLING NESTED GRID HALO'
    else
       if (is_master()) print*, 'SENDING DATA TO FILL NESTED GRID HALO'
    endif


    !Filling phis?
    !In idealized test cases, where the topography is EXACTLY known (ex case 13),
    !interpolating the topography yields a much worse result. In comparison in
    !real topography cases little difference is seen.

    !This is probably because the halo phis, which is used to compute
    !geopotential height (gz, gh), only affects the interior by being
    !used to compute corner gz in a2b_ord[24]. We might suppose this
    !computation would be more accurate when using values of phis which
    !are more consistent with those on the interior (ie the exactly-known
    !values) than the crude values given through linear interpolation.

    !For real topography cases, or in cases in which the coarse-grid topo
    ! is smoothed, we fill the boundary halo with the coarse-grid topo.

#ifndef SW_DYNAMICS
    !pt --- actually temperature

    call nested_grid_bc(atm%pt, atm%parent_grid%pt, atm%neststruct%nest_domain, &
         atm%neststruct%ind_h, atm%neststruct%wt_h, 0, 0, &
         atm%npx, atm%npy, npz, atm%bd, isg, ieg, jsg, jeg, proc_in=process)    

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

       !w
       call nested_grid_bc(atm%w(:,:,:), &
            atm%parent_grid%w(:,:,:), &
            atm%neststruct%nest_domain, atm%neststruct%ind_h, atm%neststruct%wt_h, 0, 0, &
            atm%npx, atm%npy, npz, atm%bd, isg, ieg, jsg, jeg, proc_in=process)


       !delz
       call nested_grid_bc(atm%delz(:,:,:), &
            atm%parent_grid%delz(:,:,:), &
            atm%neststruct%nest_domain, atm%neststruct%ind_h, atm%neststruct%wt_h, 0, 0, &
            atm%npx, atm%npy, npz, atm%bd, isg, ieg, jsg, jeg, proc_in=process)

    end if

#endif

    if (atm%neststruct%child_proc) then
       call nested_grid_bc(atm%u, atm%parent_grid%u(:,:,:), &
            atm%neststruct%nest_domain, atm%neststruct%ind_u, atm%neststruct%wt_u, 0, 1, &
            atm%npx, atm%npy, npz, atm%bd, isg, ieg, jsg, jeg, proc_in=process)
       call nested_grid_bc(atm%v, atm%parent_grid%v(:,:,:), &
            atm%neststruct%nest_domain, atm%neststruct%ind_v, atm%neststruct%wt_v, 1, 0, &
            atm%npx, atm%npy, npz, atm%bd, isg, ieg, jsg, jeg, proc_in=process)
    else
       call nested_grid_bc(atm%parent_grid%u(:,:,:), &
            atm%neststruct%nest_domain, 0, 1)
       call nested_grid_bc(atm%parent_grid%v(:,:,:), &
            atm%neststruct%nest_domain, 1, 0)
    endif


    if (process) then
!!$#ifdef SW_DYNAMICS
!!$    !ps: first level only
!!$    !This is only valid for shallow-water simulations
!!$    do j=jsd,jed
!!$       do i=isd,ied
!!$
!!$          Atm%ps(i,j) = Atm%delp(i,j,1)/grav
!!$
!!$       end do
!!$    end do
!!$#endif
       call mpp_update_domains(atm%u, atm%v, atm%domain, gridtype=dgrid_ne)
       call mpp_update_domains(atm%w, atm%domain, complete=.true.) ! needs an update-domain for rayleigh damping
      endif

      call mpp_sync_self()

  end subroutine setup_nested_boundary_halo

  subroutine fill_nested_grid_topo_halo(Atm, proc_in)

    type(fv_atmos_type), intent(INOUT) :: Atm
    logical, intent(IN), OPTIONAL :: proc_in
    integer :: isg, ieg, jsg, jeg

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

    call mpp_get_global_domain( atm%parent_grid%domain, &
         isg, ieg, jsg, jeg)

    if (is_master()) print*, '  FILLING NESTED GRID HALO WITH INTERPOLATED TERRAIN'
    call nested_grid_bc(atm%phis, atm%parent_grid%phis, atm%neststruct%nest_domain, &
         atm%neststruct%ind_h, atm%neststruct%wt_h, 0, 0, &
         atm%npx, atm%npy, atm%bd, isg, ieg, jsg, jeg, proc_in=proc_in)
    
  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-1)*atm%parent_grid%npes_per_tile
    if (atm%neststruct%parent_proc .and. atm%neststruct%parent_tile == atm%parent_grid%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

  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

    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%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%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%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%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%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%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%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%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 fill_nested_grid_data_end(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, atm%neststruct%nest_domain, &
               atm%neststruct%ind_update_h(isd_p:ied_p+1,jsd_p:jed_p+1,:), &
               atm%gridstruct%dx, atm%gridstruct%dy, atm%gridstruct%area, &
               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%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
    !Sets up flow to be initially hydrostatic (shouldn't be the case for all ICs?)
    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)
#endif

 

  end subroutine fill_nested_grid_data_end

  subroutine fv_write_restart(Atm, grids_on_this_pe, timestamp)
    type(fv_atmos_type), intent(inout) :: atm(:)
    character(len=*),    intent(in)    :: timestamp
    logical, intent(IN) :: grids_on_this_pe(:)
    integer n

    call fv_io_write_restart(atm, grids_on_this_pe, timestamp)
    do n=1,size(atm)
       if (atm(n)%neststruct%nested .and. grids_on_this_pe(n)) then
          call fv_io_write_bcs(atm(n))
       endif
    enddo

  end subroutine fv_write_restart

  subroutine fv_restart_end(Atm, grids_on_this_pe)
    type(fv_atmos_type), intent(inout) :: atm(:)
    logical, intent(INOUT) :: grids_on_this_pe(:)

    integer :: isc, iec, jsc, jec
    integer :: iq, n, ntileme, 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


    ntileme = size(atm(:))

    do n = 1, ntileme

      if (.not. grids_on_this_pe(n)) then
         cycle
      endif

      call mpp_set_current_pelist(atm(n)%pelist)

      isc = atm(n)%bd%isc; iec = atm(n)%bd%iec; jsc = atm(n)%bd%jsc; jec = atm(n)%bd%jec

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

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

      unit = stdout()
      write(unit,*)
      write(unit,*) 'fv_restart_end u   ', trim(gn),' = ', mpp_chksum(atm(n)%u(isc:iec,jsc:jec,:))
      write(unit,*) 'fv_restart_end v   ', trim(gn),' = ', mpp_chksum(atm(n)%v(isc:iec,jsc:jec,:))
      if ( .not. atm(n)%flagstruct%hydrostatic )    &
         write(unit,*) 'fv_restart_end w   ', trim(gn),' = ', mpp_chksum(atm(n)%w(isc:iec,jsc:jec,:))
      write(unit,*) 'fv_restart_end delp', trim(gn),' = ', mpp_chksum(atm(n)%delp(isc:iec,jsc:jec,:))
      write(unit,*) 'fv_restart_end phis', trim(gn),' = ', mpp_chksum(atm(n)%phis(isc:iec,jsc:jec))
#ifndef SW_DYNAMICS
      write(unit,*) 'fv_restart_end pt  ', trim(gn),' = ', mpp_chksum(atm(n)%pt(isc:iec,jsc:jec,:))
      if (ntprog>0) &
           write(unit,*) 'fv_restart_end q(prog) nq  ', trim(gn),' =',ntprog, mpp_chksum(atm(n)%q(isc:iec,jsc:jec,:,:))
      if (ntdiag>0) &
           write(unit,*) 'fv_restart_end 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_end '//trim(tracer_name)// trim(gn),' = ', mpp_chksum(atm(n)%q(isc:iec,jsc:jec,:,iq))
      enddo

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

   enddo

   call fv_io_write_restart(atm, grids_on_this_pe)
   do n=1,ntileme
      if (atm(n)%neststruct%nested .and. grids_on_this_pe(n)) call fv_io_write_bcs(atm(n))
    end do

    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 d2c_setup(u, v, &
      ua, va, &
      uc, vc, dord4, &
      isd,ied,jsd,jed, is,ie,js,je, npx,npy, &
      grid_type, nested, &
      se_corner, sw_corner, ne_corner, nw_corner, &
      rsin_u,rsin_v,cosa_s,rsin2,regional )

  logical, intent(in):: dord4
  real, intent(in) ::  u(isd:ied,jsd:jed+1)
  real, intent(in) ::  v(isd:ied+1,jsd:jed)
  real, intent(out), dimension(isd:ied  ,jsd:jed  ):: ua
  real, intent(out), dimension(isd:ied  ,jsd:jed  ):: va
  real, intent(out), dimension(isd:ied+1,jsd:jed  ):: uc
  real, intent(out), dimension(isd:ied  ,jsd:jed+1):: vc
  integer, intent(in) :: isd,ied,jsd,jed, is,ie,js,je, npx,npy,grid_type
  logical, intent(in) :: nested, se_corner, sw_corner, ne_corner, nw_corner, regional
  real, intent(in) :: rsin_u(isd:ied+1,jsd:jed)
  real, intent(in) :: rsin_v(isd:ied,jsd:jed+1)
  real, intent(in) :: cosa_s(isd:ied,jsd:jed)
  real, intent(in) :: rsin2(isd:ied,jsd:jed)

! Local 
  real, dimension(isd:ied,jsd:jed):: utmp, vtmp
  real, parameter:: t11=27./28., t12=-13./28., t13=3./7., t14=6./7., t15=3./28.
  real, parameter:: a1 =  0.5625
  real, parameter:: a2 = -0.0625
  real, parameter:: c1 = -2./14.
  real, parameter:: c2 = 11./14.
  real, parameter:: c3 =  5./14.
  integer npt, i, j, ifirst, ilast, id

  if ( dord4) then
       id = 1
  else
       id = 0
  endif


  if (grid_type < 3 .and. .not. (nested .or. regional)) then
     npt = 4
  else
     npt = -2
  endif

  if ( nested) then  

     do j=jsd+1,jed-1
        do i=isd,ied
           utmp(i,j) = a2*(u(i,j-1)+u(i,j+2)) + a1*(u(i,j)+u(i,j+1))
        enddo
     enddo
     do i=isd,ied
        j = jsd
        utmp(i,j) = 0.5*(u(i,j)+u(i,j+1))
        j = jed
        utmp(i,j) = 0.5*(u(i,j)+u(i,j+1))
     end do

     do j=jsd,jed
        do i=isd+1,ied-1
           vtmp(i,j) = a2*(v(i-1,j)+v(i+2,j)) + a1*(v(i,j)+v(i+1,j))
        enddo
        i = isd
        vtmp(i,j) = 0.5*(v(i,j)+v(i+1,j)) 
        i = ied
        vtmp(i,j) = 0.5*(v(i,j)+v(i+1,j))
     enddo

     do j=jsd,jed
        do i=isd,ied
           ua(i,j) = (utmp(i,j)-vtmp(i,j)*cosa_s(i,j)) * rsin2(i,j)
           va(i,j) = (vtmp(i,j)-utmp(i,j)*cosa_s(i,j)) * rsin2(i,j)
        enddo
     enddo

  else

     !----------
     ! Interior:
     !----------
     utmp = 0.
     vtmp = 0.


     do j=max(npt,js-1),min(npy-npt,je+1)
        do i=max(npt,isd),min(npx-npt,ied)
           utmp(i,j) = a2*(u(i,j-1)+u(i,j+2)) + a1*(u(i,j)+u(i,j+1))
        enddo
     enddo
     do j=max(npt,jsd),min(npy-npt,jed)
        do i=max(npt,is-1),min(npx-npt,ie+1)
           vtmp(i,j) = a2*(v(i-1,j)+v(i+2,j)) + a1*(v(i,j)+v(i+1,j))
        enddo
     enddo

     !----------
     ! edges:
     !----------
     if (grid_type < 3) then

        if ( js==1 .or. jsd<npt) then
           do j=jsd,npt-1
              do i=isd,ied
                 utmp(i,j) = 0.5*(u(i,j) + u(i,j+1))
                 vtmp(i,j) = 0.5*(v(i,j) + v(i+1,j))
              enddo
           enddo
        endif

        if ( (je+1)==npy .or. jed>=(npy-npt)) then
           do j=npy-npt+1,jed
              do i=isd,ied
                 utmp(i,j) = 0.5*(u(i,j) + u(i,j+1))
                 vtmp(i,j) = 0.5*(v(i,j) + v(i+1,j))
              enddo
           enddo
        endif

        if ( is==1 .or. isd<npt ) then
           do j=max(npt,jsd),min(npy-npt,jed)
              do i=isd,npt-1
                 utmp(i,j) = 0.5*(u(i,j) + u(i,j+1))
                 vtmp(i,j) = 0.5*(v(i,j) + v(i+1,j))
              enddo
           enddo
        endif

        if ( (ie+1)==npx .or. ied>=(npx-npt)) then
           do j=max(npt,jsd),min(npy-npt,jed)
              do i=npx-npt+1,ied
                 utmp(i,j) = 0.5*(u(i,j) + u(i,j+1))
                 vtmp(i,j) = 0.5*(v(i,j) + v(i+1,j))
              enddo
           enddo
        endif

     endif
     do j=js-1-id,je+1+id
        do i=is-1-id,ie+1+id
           ua(i,j) = (utmp(i,j)-vtmp(i,j)*cosa_s(i,j)) * rsin2(i,j)
           va(i,j) = (vtmp(i,j)-utmp(i,j)*cosa_s(i,j)) * rsin2(i,j)
        enddo
     enddo

  end if

! A -> C
!--------------
! Fix the edges
!--------------
! Xdir:
     if( sw_corner ) then
         do i=-2,0
            utmp(i,0) = -vtmp(0,1-i)
         enddo
     endif
     if( se_corner ) then
         do i=0,2
            utmp(npx+i,0) = vtmp(npx,i+1)
         enddo
     endif
     if( ne_corner ) then
         do i=0,2
            utmp(npx+i,npy) = -vtmp(npx,je-i)
         enddo
     endif
     if( nw_corner ) then
         do i=-2,0
            utmp(i,npy) = vtmp(0,je+i)
         enddo
     endif

  if (grid_type < 3 .and. .not. (nested .or. regional)) then
     ifirst = max(3,    is-1)
     ilast  = min(npx-2,ie+2)
  else
     ifirst = is-1
     ilast  = ie+2
  endif
!---------------------------------------------
! 4th order interpolation for interior points:
!---------------------------------------------
     do j=js-1,je+1
        do i=ifirst,ilast
           uc(i,j) = a1*(utmp(i-1,j)+utmp(i,j))+a2*(utmp(i-2,j)+utmp(i+1,j))
        enddo
     enddo

     if (grid_type < 3) then
! Xdir:
     if( is==1 .and. .not. (nested .or. regional) ) then
        do j=js-1,je+1
           uc(0,j) = c1*utmp(-2,j) + c2*utmp(-1,j) + c3*utmp(0,j) 
           uc(1,j) = ( t14*(utmp( 0,j)+utmp(1,j))    &
                     + t12*(utmp(-1,j)+utmp(2,j))    &
                     + t15*(utmp(-2,j)+utmp(3,j)) )*rsin_u(1,j)
           uc(2,j) = c1*utmp(3,j) + c2*utmp(2,j) + c3*utmp(1,j)
        enddo
     endif

     if( (ie+1)==npx .and. .not. (nested .or. regional) ) then
        do j=js-1,je+1
           uc(npx-1,j) = c1*utmp(npx-3,j)+c2*utmp(npx-2,j)+c3*utmp(npx-1,j) 
           uc(npx,j) = (t14*(utmp(npx-1,j)+utmp(npx,j))+      &
                        t12*(utmp(npx-2,j)+utmp(npx+1,j))     &
                      + t15*(utmp(npx-3,j)+utmp(npx+2,j)))*rsin_u(npx,j)
           uc(npx+1,j) = c3*utmp(npx,j)+c2*utmp(npx+1,j)+c1*utmp(npx+2,j) 
        enddo
     endif

     endif

!------
! Ydir:
!------
     if( sw_corner ) then
         do j=-2,0
            vtmp(0,j) = -utmp(1-j,0)
         enddo
     endif
     if( nw_corner ) then
         do j=0,2
            vtmp(0,npy+j) = utmp(j+1,npy)
         enddo
     endif
     if( se_corner ) then
         do j=-2,0
            vtmp(npx,j) = utmp(ie+j,0)
         enddo
     endif
     if( ne_corner ) then
         do j=0,2
            vtmp(npx,npy+j) = -utmp(ie-j,npy)
         enddo
     endif

     if (grid_type < 3) then

     do j=js-1,je+2
      if ( j==1  .and. .not. (nested .or. regional)) then
        do i=is-1,ie+1
           vc(i,1) = (t14*(vtmp(i, 0)+vtmp(i,1))    &
                    + t12*(vtmp(i,-1)+vtmp(i,2))    &
                    + t15*(vtmp(i,-2)+vtmp(i,3)))*rsin_v(i,1)
        enddo
      elseif ( (j==0 .or. j==(npy-1))  .and. .not. (nested .or. regional)) then
        do i=is-1,ie+1
           vc(i,j) = c1*vtmp(i,j-2) + c2*vtmp(i,j-1) + c3*vtmp(i,j)
        enddo
      elseif ( (j==2 .or. j==(npy+1))  .and. .not. (nested .or. regional)) then
        do i=is-1,ie+1
           vc(i,j) = c1*vtmp(i,j+1) + c2*vtmp(i,j) + c3*vtmp(i,j-1)
        enddo
      elseif ( j==npy  .and. .not. (nested .or. regional)) then
        do i=is-1,ie+1
           vc(i,npy) = (t14*(vtmp(i,npy-1)+vtmp(i,npy))    &
                      + t12*(vtmp(i,npy-2)+vtmp(i,npy+1))  &
                      + t15*(vtmp(i,npy-3)+vtmp(i,npy+2)))*rsin_v(i,npy)
        enddo
      else
! 4th order interpolation for interior points:
        do i=is-1,ie+1
           vc(i,j) = a2*(vtmp(i,j-2)+vtmp(i,j+1))+a1*(vtmp(i,j-1)+vtmp(i,j))
        enddo
     endif
     enddo
    else
! 4th order interpolation:
       do j=js-1,je+2
          do i=is-1,ie+1
             vc(i,j) = a2*(vtmp(i,j-2)+vtmp(i,j+1))+a1*(vtmp(i,j-1)+vtmp(i,j))
          enddo
       enddo
    endif

  end subroutine d2c_setup

 subroutine d2a_setup(u, v, ua, va, dord4, &
      isd,ied,jsd,jed, is,ie,js,je, npx,npy, &
      grid_type, nested, &
      cosa_s,rsin2,regional )

  logical, intent(in):: dord4
  real, intent(in) ::  u(isd:ied,jsd:jed+1)
  real, intent(in) ::  v(isd:ied+1,jsd:jed)
  real, intent(out), dimension(isd:ied  ,jsd:jed  ):: ua
  real, intent(out), dimension(isd:ied  ,jsd:jed  ):: va
  integer, intent(in) :: isd,ied,jsd,jed, is,ie,js,je, npx,npy,grid_type
  real, intent(in) :: cosa_s(isd:ied,jsd:jed)
  real, intent(in) :: rsin2(isd:ied,jsd:jed)
  logical, intent(in) :: nested, regional

! Local 
  real, dimension(isd:ied,jsd:jed):: utmp, vtmp
  real, parameter:: t11=27./28., t12=-13./28., t13=3./7., t14=6./7., t15=3./28.
  real, parameter:: a1 =  0.5625
  real, parameter:: a2 = -0.0625
  real, parameter:: c1 = -2./14.
  real, parameter:: c2 = 11./14.
  real, parameter:: c3 =  5./14.
  integer npt, i, j, ifirst, ilast, id

  if ( dord4) then
       id = 1
  else
       id = 0
  endif


  if (grid_type < 3 .and. .not. (nested .or. regional)) then
     npt = 4
  else
     npt = -2
  endif

  if ( nested) then  

     do j=jsd+1,jed-1
        do i=isd,ied
           utmp(i,j) = a2*(u(i,j-1)+u(i,j+2)) + a1*(u(i,j)+u(i,j+1))
        enddo
     enddo
     do i=isd,ied
        j = jsd
        utmp(i,j) = 0.5*(u(i,j)+u(i,j+1))
        j = jed
        utmp(i,j) = 0.5*(u(i,j)+u(i,j+1))
     end do

     do j=jsd,jed
        do i=isd+1,ied-1
           vtmp(i,j) = a2*(v(i-1,j)+v(i+2,j)) + a1*(v(i,j)+v(i+1,j))
        enddo
        i = isd
        vtmp(i,j) = 0.5*(v(i,j)+v(i+1,j)) 
        i = ied
        vtmp(i,j) = 0.5*(v(i,j)+v(i+1,j))
     enddo

  else

     !----------
     ! Interior:
     !----------

     do j=max(npt,js-1),min(npy-npt,je+1)
        do i=max(npt,isd),min(npx-npt,ied)
           utmp(i,j) = a2*(u(i,j-1)+u(i,j+2)) + a1*(u(i,j)+u(i,j+1))
        enddo
     enddo
     do j=max(npt,jsd),min(npy-npt,jed)
        do i=max(npt,is-1),min(npx-npt,ie+1)
           vtmp(i,j) = a2*(v(i-1,j)+v(i+2,j)) + a1*(v(i,j)+v(i+1,j))
        enddo
     enddo

     !----------
     ! edges:
     !----------
     if (grid_type < 3) then

        if ( js==1 .or. jsd<npt) then
           do j=jsd,npt-1
              do i=isd,ied
                 utmp(i,j) = 0.5*(u(i,j) + u(i,j+1))
                 vtmp(i,j) = 0.5*(v(i,j) + v(i+1,j))
              enddo
           enddo
        endif

        if ( (je+1)==npy .or. jed>=(npy-npt)) then
           do j=npy-npt+1,jed
              do i=isd,ied
                 utmp(i,j) = 0.5*(u(i,j) + u(i,j+1))
                 vtmp(i,j) = 0.5*(v(i,j) + v(i+1,j))
              enddo
           enddo
        endif

        if ( is==1 .or. isd<npt ) then
           do j=max(npt,jsd),min(npy-npt,jed)
              do i=isd,npt-1
                 utmp(i,j) = 0.5*(u(i,j) + u(i,j+1))
                 vtmp(i,j) = 0.5*(v(i,j) + v(i+1,j))
              enddo
           enddo
        endif

        if ( (ie+1)==npx .or. ied>=(npx-npt)) then
           do j=max(npt,jsd),min(npy-npt,jed)
              do i=npx-npt+1,ied
                 utmp(i,j) = 0.5*(u(i,j) + u(i,j+1))
                 vtmp(i,j) = 0.5*(v(i,j) + v(i+1,j))
              enddo
           enddo
        endif

     endif

  end if



  do j=js-1-id,je+1+id
     do i=is-1-id,ie+1+id
        ua(i,j) = (utmp(i,j)-vtmp(i,j)*cosa_s(i,j)) * rsin2(i,j)
        va(i,j) = (vtmp(i,j)-utmp(i,j)*cosa_s(i,j)) * rsin2(i,j)
     enddo
  enddo

end subroutine d2a_setup

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