module fv_restart_mod

  use fms_mod,    only: file_exist, error_mesg,  lowercase, FATAL, &
       NOTE, read_data, uppercase, stdout, mpp_pe, mpp_root_pe, mpp_chksum
  use fms_io_mod,   only: open_restart_file, close_file, write_data, get_instance_filename
  use mpp_io_mod, only: mpp_open, mpp_close, axistype, fieldtype,  &
       mpp_read_meta, mpp_get_info, mpp_get_fields, mpp_read, &
       MPP_NETCDF, MPP_SINGLE, MPP_MULTI, MPP_SEQUENTIAL, MPP_RDONLY, &
       MPP_WRONLY, MPP_NATIVE, MPP_OVERWR, mpp_write_meta, mpp_write, &
       atttype, mpp_get_atts, mpp_get_att_name, mpp_get_att_real_scalar, &
       mpp_get_field_name
  use fv_pack,    only: ak, bk, ks, ptop, beglat, endlat,       &
       master, nlon, mlat, nt_phys, nqrst, ng_d, ng_s, &
       ighost, mountain, coslon, sinlon, lon, lat, lonb, latb,     &
       adjust_dry_mass, use_set_eta, p_var, d2a3d,          &
       age_tracer, age_time, get_eta_level, &
       do_fms_tracer_manager, fms_tracers_file, fv_domain, &
       restart_format, use_tendency, cold_start, gid
#ifdef MARS_GCM
  use fv_pack,    only: p_ref
#endif MARS_GCM

  use tracer_manager_mod, only : tr_get_tracer_names=>get_tracer_names, &
       get_tracer_names, get_number_tracers, &
       set_tracer_profile, &
       get_tracer_index, NO_TRACER
  use field_manager_mod, only  : MODEL_ATMOS      
  use         platform_mod, only: I8_KIND

!  use shr_kind_mod, only : r8 => shr_kind_r8
#ifdef SPMD
  use mod_comm, only: mp_bcst_int, mp_bcst_r2d, mp_scatter3d, mp_bcst_real, &
       mp_send4d_ns, mp_recv4d_ns, mp_gather3d
#endif

  use fv_arrays_mod, only: fv_array_check, fv_stack_push, fv_array_sync, &
       fv_print_chksums, ksp, kep
  use mpp_mod, only: mpp_sync, mpp_root_pe, mpp_broadcast
  use mpp_domains_mod, only: mpp_get_data_domain, mpp_get_compute_domain
  use constants_mod,    only: grav, kappa, rdgas

  use pmaxmin_mod, only: pmaxming, pmaxmin

  implicit none
  private

  public :: fv_restart, read_fv_rst, write_fv_rst, add_tracers

contains

  subroutine read_fv_rst(im, jm, km, nq, dflnm, days, seconds, oform)

#include "fv_arrays.h"
    character*(*), intent(in)       :: dflnm                ! dynamics filename
    integer, intent(in) :: nq
    integer, intent(out):: im, jm, km
    integer, intent(out):: days
    integer, intent(out):: seconds
    character(len=*), intent(in) :: oform

! Local variables
    integer iuic                       ! Unit number
    integer i, j, k, n, jn, it, ntracers
    logical is_tracer_in_restart_file
    character(len=128)    :: tname
    real zsurf(nlon, beglat: endlat)
    real ginv, qmax, qmin
    real, allocatable :: w2d(:,:), w3d(:,:,:)

    type(fieldtype), allocatable :: tracer_fields(:), fields(:)
    type(atttype), allocatable :: attributes(:)
    integer :: ndim,natt,nvar,ntime
    real p0

    character(len=*), parameter :: fms_tracers_restart_file = &
         'INPUT/'//fms_tracers_file

    character(len=128) :: directory, restart_name, tr_name, dflnm_mod
    logical :: use_native_format
    integer :: ids, ide, jds, jde, ics,ice,jcs,jce

#ifdef MARS_GCM
!      ------Dont require specific water vapor fields in netcdf version 
!         of tracer restart file:   Hence this code is not necessary 
#else
    integer :: nsphum, nliq_wat, nice_wat, ncld_amt
    character(len=20) :: fv_tr_name(4)
    data fv_tr_name / "sphum","liq_wat","ice_wat","cld_amt" /
#endif


#ifdef use_shared_pointers
    pointer( p_zsurf, zsurf )
#include "fv_point.inc"
    call fv_stack_push( p_zsurf, nlon*(endlat-beglat+1) )
#endif

    ginv = 1./ GRAV

    !Modify the filename if it has an ensemble instance appendix
    !The following call does not modify the filename unless there 
    !is an ensemble instance appendix and fms_io knows about it.
    call get_instance_filename(dflnm, dflnm_mod)

! output format
    if (trim(uppercase(oform)) == 'NATIVE') then
        use_native_format = .true.
        restart_name = trim(dflnm_mod)
    else if (trim(uppercase(oform)) == 'NETCDF') then
        use_native_format = .false.
        if (trim(dflnm_mod(LEN_TRIM(dflnm_mod)-2:LEN_TRIM(dflnm_mod))) /= ".nc" ) &
             restart_name = trim(dflnm_mod)//".nc"
        if (file_exist(dflnm_mod) .and. .not. file_exist(trim(dflnm_mod)//".nc")) then
            call error_mesg( 'read_fv_rst', &
                 'NetCDF file does not exist. Switching to Native format.', &
                 NOTE)
            restart_name = trim(dflnm_mod)
            use_native_format = .true.
        endif

    else
        call error_mesg ('read_fv_rst', &
             'invalid value for argument oform='//trim(oform), FATAL)
    endif
    if (.not. file_exist(dflnm_mod) .and. .not. file_exist(trim(dflnm_mod)//".nc")) &
         call error_mesg ('read_fv_rst', &
         'Neither a Native or NetCDF restart file is present.', FATAL)

    if( master ) then
        if (use_native_format) then
            if (file_exist(dflnm_mod)) then
!!$         call mpp_open (iuic, dflnm_mod, action=MPP_RDONLY, form=MPP_NATIVE, access=MPP_SEQUENTIAL)
                iuic = open_restart_file(dflnm_mod, 'read')
                read(iuic) im, jm, km, nqrst
            endif
        else
            if (file_exist(trim(dflnm_mod)//'.nc')) then
                call mpp_open(iuic, dflnm_mod, &
                     action=MPP_RDONLY, form=MPP_NETCDF, threading=MPP_SINGLE )
                call mpp_read_meta(iuic)
                call mpp_get_info(iuic,ndim,nvar,natt,ntime)

                allocate(attributes(natt))
                call mpp_get_atts(iuic,attributes)

                do n= 1,natt
!            select case (attributes(n)%name)
                   select case ( mpp_get_att_name(attributes(n)))
                   case ('nlon')
                       im= int(mpp_get_att_real_scalar(attributes(n)))
                   case ('mlat')
                       jm= int(mpp_get_att_real_scalar(attributes(n)))
                   case ('nlev')
                       km= int(mpp_get_att_real_scalar(attributes(n)))
                   case ('nt_phys')
                       nqrst= int(mpp_get_att_real_scalar(attributes(n)))
                   case ('days')
                       days= int(mpp_get_att_real_scalar(attributes(n)))
                   case ('seconds')
                       seconds= int(mpp_get_att_real_scalar(attributes(n)))
                   case ('ks')
                       ks = int(mpp_get_att_real_scalar(attributes(n)))
                   end select
                enddo
                deallocate(attributes)


                allocate(fields(nvar))
                call mpp_get_fields(iuic, fields)
                do n=1,nvar
!            select case (fields(n)%name)
                   select case ( mpp_get_field_name(fields(n)))
                   case ('ak')
                       call mpp_read(iuic,fields(n),ak)
                   case ('bk')
                       call mpp_read(iuic,fields(n),bk)
                   end select
                enddo
                call mpp_close(iuic)
                deallocate(fields)
            endif
        endif
        if ( im /= nlon .or. jm /= mlat .or. km /= nlev ) then
            write(6,*) 'Dataset resolution:',im, jm, km, nqrst
            call error_mesg('read_fv_rst:','resolution inconsistent', FATAL)
        endif

        if ( nqrst > nt_phys .or. nqrst > nq) then
            write(6,*) ' '
            write(6,*) '---- Warning -------'
            write(6,*) 'Total number of tracers in restart file is greater than either'
            write(6,*) 'NT_PHYS or NQ.  The number NT_PHYS from the namelist is used.'
            write(6,*) 'NQRST=', nqrst, 'NT_PHYS=', nt_phys, ' NQ=', nq
            nqrst = nt_phys
        endif
        if (use_native_format) then
            read(iuic) days, seconds
            read(iuic) ak, bk, ks
        endif
        write(6,*) 'Restart days=', days,' sec=',seconds
        write(6,*) '----------------------------------'
        write(6,*) 'Model Top Pressure (pa)= ', ak(1)
        write(6,*) '----------------------------------'
        write(6,*) 'Checking Initial condition ...'
    endif

#ifdef SPMD           
!    call mp_bcst_int(days)
!    call mp_bcst_int(seconds)
!    call mp_bcst_int(im)
!    call mp_bcst_int(jm)
!    call mp_bcst_int(km)
!    call mp_bcst_int(nqrst)
!    call mp_bcst_int(ks)
!
!    call mp_bcst_r2d(km+1, 1, 1, 1, ak, 0)
!    call mp_bcst_r2d(km+1, 1, 1, 1, bk, 0)
    call mpp_broadcast( days, mpp_root_pe() )
    call mpp_broadcast( seconds, mpp_root_pe() )
    call mpp_broadcast( im, mpp_root_pe() )
    call mpp_broadcast( jm, mpp_root_pe() )
    call mpp_broadcast( km, mpp_root_pe() )
    call mpp_broadcast( nqrst, mpp_root_pe() )
    call mpp_broadcast( ks, mpp_root_pe() )

    call mpp_broadcast( ak, km+1, mpp_root_pe() )
    call mpp_broadcast( bk, km+1, mpp_root_pe() )
#endif

    ptop = ak(1)

    if (use_native_format) then   
!---------------------------
! read surface geopotential
!---------------------------
        allocate ( w2d(im,jm) )
#ifdef SPMD
        if ( master ) read(iuic) w2d
        call mp_scatter3d(w2d, phis, im,  jm, 1, beglat, endlat,  &
             1, 1, ighost, ighost, 0)
        call fv_array_sync
#else
        read(iuic) w2d
        do j=beglat,endlat
           do i=1,im
              phis(i,j) = w2d(i,j)
           enddo
        enddo
#endif
!       call pmaxming('PHIS', phis, im, jm, 1, beglat, endlat,   &
!                     ighost, ighost, ginv)

        allocate ( w3d(im, jm, km) )

!---------------------------
! Read u-wind
!---------------------------

#ifdef SPMD
        if ( master ) read(iuic) w3d
        call mp_scatter3d(w3d, u, im,  jm, km, beglat, endlat,  &
             1, km, ng_d, ng_s, 0)
        call fv_array_sync
#else
        read(iuic) u
#endif

!---------------------------
! Read v-wind
!---------------------------

#ifdef SPMD
        if ( master ) read(iuic) w3d
        call mp_scatter3d(w3d, v, im,  jm, km, beglat, endlat,  &
             1, km, ng_d, ng_d, 0)
        call fv_array_sync
#else
        read(iuic) v
#endif

!-----------------------------------
! Read virtual potential temperature
!-----------------------------------

#ifdef SPMD
        if ( master ) read(iuic) w3d
        call mp_scatter3d(w3d, pt, im,  jm, km, beglat, endlat,  &
             1, km, ng_d, ng_d, 0)
        call fv_array_sync
#else
        read(iuic) pt
#endif

!-----------------------------------
! Read pressure thickness (pascal)
!-----------------------------------

#ifdef SPMD
        if ( master ) read(iuic) w3d
        call mp_scatter3d(w3d, delp, im,  jm, km, beglat, endlat,  &
             1, km, 0, 0, 0)
        call fv_array_sync
#else
        read(iuic) delp
#endif

    else  
!this code may not work? all these fields do not have the same data domain
        call mpp_get_data_domain( fv_domain, ids, ide, jds, jde )
        call mpp_get_compute_domain( fv_domain, ics, ice, jcs, jce )
        call read_data ( restart_name, 'Surface_geopotential', phis(ids:ide,jcs:jce), domain = fv_domain)
        call read_data ( restart_name, 'U', u(ids:ide,jds:jde,:) , domain = fv_domain)
        call read_data ( restart_name, 'V', v(ids:ide,jds:jde,:), domain = fv_domain)
        call read_data ( restart_name, 'T', pt(ids:ide,jds:jde,:), domain = fv_domain)
        call read_data ( restart_name, 'DELP', delp(ids:ide,jcs:jce,:), domain = fv_domain) !RASF Bug fix for indices
    endif ! not use_native_format

#ifdef SHAVE_P
    call set_eta(km, ks, ptop, ak, bk)
    do j=beglat,endlat
       do i=1,im
          delp(i,j,1) = ak(2) - ak(1)
       enddo
    enddo
#endif

!------------------------------------
! Read tracers one at a time
!------------------------------------

    if (use_native_format) then 
        do n=1, nqrst

#ifdef SPMD
           if ( master ) read(iuic) w3d
           call mp_scatter3d(w3d, q(1,beglat-ng_d,1,n), im,  jm, km, beglat, endlat,  &
                1, km, ng_d, ng_d , 0)
           call fv_array_sync
#else
           read(iuic) w3d
           do k=1,nlev
              do j=beglat,endlat
                 do i=1,im
                    q(i,j,k,n) = w3d(i,j,k)
                 enddo
              enddo
           enddo
#endif
        enddo
    else        ! ----------read "dynamics" tracers  from tracer restart file 
#ifdef MARS_GCM
!rjw      ------No longer require specific water vapor fields in netcdf version of tracer restart file
!rjw          It now makes more sense to read all the tracers from the tracer file
!rjw           in a single loop, that appears below in the processing of field manager tracers
#else
        do n =1,nqrst
           call read_data(fms_tracers_restart_file,fv_tr_name(n),q(ids:ide,jds:jde,:,n) , domain = fv_domain)
        enddo
#endif MARS_GCM 
    endif ! --------------- Finished reading dynamcics tracers ------------


    if( age_tracer ) call error_mesg('read_fv_rst', 'age_tracer == .T. not yet working', FATAL)

    if (use_native_format) then 
        if ( age_tracer .and. nqrst == nq ) then  ! this check will not work properly (pletzer)
!           q = 0. ; age_time = 0.
            if(master) read(iuic) age_time
#ifdef SPMD
!            call mp_bcst_real(age_time)
            call mpp_broadcast( age_time, mpp_root_pe() )
#endif
        endif

        if ( master ) call close_file(iuic)
    endif

!-----------------------------------------
! Read additional tracers from netCDF file
!-----------------------------------------

#ifdef MARS_GCM
    if (do_fms_tracer_manager ) then !{
#else
    if (do_fms_tracer_manager .and. nq > nt_phys) then !{
#endif MARS_GCM

       !Modify the filename if it has an ensemble instance appendix
       call get_instance_filename(fms_tracers_restart_file,restart_name )
       restart_name=trim(restart_name)
        if( .not. file_exist(restart_name//'.nc') ) then !{

            if(master) then
                write(*,*)'---- Warning -------'
                write(*,*)'do_fms_tracer_manager=',do_fms_tracer_manager, &
                     ' no of tracers nq=', nq
                write(*,*)'but netcdf restart file ', restart_name//'.nc', ' does not exist'
                write(*,*)'--------------------'               
            endif

        else !}{

            if(master) then
                write(*,*)'Additional tracers present in field_table. '
                write(*,*)'These will be read from netcdf file ',restart_name
            endif

            call mpp_open(iuic, restart_name, &
                 action=MPP_RDONLY, form=MPP_NETCDF, threading=MPP_SINGLE )

            call mpp_read_meta(iuic)
            call mpp_get_info(iuic,ndim,nvar,natt,ntime)
            if ( nt_phys+nvar /= nq ) then
                if (master) then 
                    write(6,*) ' '
                    write(6,*) '---- Warning -------'
                    write(6,*) 'Total number of tracers nq from netcdf and binary restart files does not match'
                    write(6,*) 'nt_phys=', nt_phys,' nvar=', nvar,' nq=', nq, ' nt_phys+nvar /= nq'
                    write(*,*)'---------------------'        
                endif
            endif
            allocate(tracer_fields(nvar))
!if(master) 
            call mpp_get_fields(iuic, tracer_fields)

! read only the tracers that are registered in the tracer manager
            call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers)
            n = 0 
            call mpp_get_data_domain( fv_domain, ids, ide, jds, jde )
            do it = 1, ntracers !{
! loop over tracers
               call get_tracer_names(MODEL_ATMOS, it, tname)
#ifdef MARS_GCM

#else
               if (use_native_format) then
! exclude tracers read from native format file 
                   if(tname == "sphum"  ) cycle
                   if(tname == "liq_wat") cycle
                   if(tname == "ice_wat") cycle
                   if(tname == "cld_amt") cycle
               endif
#endif MARS_GCM
! this tracer will be added
               n = n + 1
               is_tracer_in_restart_file = .FALSE.
               do jn = 1, nvar !{
! loop over netcdf file tracers                  
!                  if(  lowercase( trim(tracer_fields(jn)%name) ) == lowercase( trim(tname) )  ) then !{
                  if(  lowercase( trim(mpp_get_field_name(tracer_fields(jn))) ) == lowercase( trim(tname) )  ) then !{
                      if(master) print *,'tracer ', trim(tname), ' found!'
! ok found the tracer in the restart file
                      is_tracer_in_restart_file = .TRUE.
! read data and store in array w3d
!                     call mpp_read(iuic, tracer_fields(jn), w3d)
                      call read_data( restart_name, tname, &
                           q(ids:ide,jds:jde,:,it), domain = fv_domain)
                  endif !}
                  if(is_tracer_in_restart_file) exit
               enddo !}
               if(.NOT. is_tracer_in_restart_file) then !{
! missing tracer from restart file
                   if(master) then
                       write(6,*) '---- Warning -------'
                       write(6,*) 'Tracer ', tname, ' was defined in field table but could not be found in restart file.'
                       write(*,*)'---------------------'        
                   endif
! set default profile
                   call set_tracer_profile (MODEL_ATMOS, it, &
                        q(ids:ide,jds:jde,:,it))
               endif !}
            enddo !}


            call mpp_close(iuic)
            deallocate(tracer_fields)

        endif !}

    endif !}

    if (use_native_format) deallocate ( w3d )
!..............................................................................................

    do j=beglat,endlat
       do i=1,im
          zsurf(i,j) = phis(i,j)
       enddo
    enddo

    call pmaxmin('ZS', zsurf, qmin, qmax, im*(endlat-beglat+1), 1, ginv)
    if ( qmax > 1. ) then
        mountain = .true.
    else
        mountain = .false.
    endif

    call pmaxming('U ',       u, im, jm, km, beglat, endlat, ng_d, ng_s, 1.)
    call pmaxming('V ',       v, im, jm, km, beglat, endlat, ng_d, ng_d, 1.)
!    call pmaxming('DELP ', delp, im, jm, km, beglat, endlat, 0,       0, 1.)
    call pmaxming('T ',      pt, im, jm, km, beglat, endlat, ng_d, ng_d, 1.)
    do n=1, nqrst
       call pmaxming('Q ', q(1,beglat-ng_d,1,n), im, jm, km, beglat, endlat, ng_d, ng_d, 1.)
    enddo

!-----------------------
! Add additional tracers:
!-----------------------

#ifdef SET_Q_TOP
! Set water vapor at top layer to climatology value (for stratosphere)
    do j=beglat, endlat
       do i=1,im
          q(i,j,1,1) = 3.E-6
       enddo
    enddo
#endif

    if( use_set_eta ) call set_eta(km, ks, ptop, ak,  bk)

    call p_var(nlon, mlat, nlev, beglat, endlat, ptop, delp, ps,           &
         pe, peln,  pk,  pkz,  kappa, q, ng_d, ncnst, adjust_dry_mass )

!---------------------------------------------------
! Read surface pressure & winds for coupling purpose
!---------------------------------------------------

    !Modify the filename if it has an ensemble instance appendix
    call get_instance_filename('INPUT/fv_srf_wnd.res',restart_name )
    restart_name=trim(restart_name)

    if (file_exist(restart_name) .or. file_exist(trim(restart_name)//'.nc')) then
        if (use_native_format) then 
            if( master ) then
!!$         call mpp_open (iuic, restart_name, action=MPP_RDONLY,      &
!!$                        form=MPP_NATIVE, access=MPP_SEQUENTIAL)
                iuic = open_restart_file(restart_name, 'read')
                read(iuic) im, jm
                if ( im /= nlon .or. jm /= mlat ) then
                    write(6,*) 'Dataset resolution:',im, jm
                    call error_mesg('read_fv_rst:','resolution inconsistent in fv_srf_wnd', FATAL)
                endif
            endif

! Read ps_bp (surface pressure before physics update):
            if(master) read(iuic) w2d
#ifdef SPMD
            call mp_scatter3d(w2d, ps_bp, im,  jm, 1, beglat, endlat,  &
                 1, 1, 0, 0, 0)
            call fv_array_sync
#endif
! Read u-comp
            if(master) read(iuic) w2d
#ifdef SPMD
            call mp_scatter3d(w2d, u_srf, im,  jm, 1, beglat, endlat,  &
                 1, 1, 0, 0, 0)
            call fv_array_sync
#endif
! Read v-comp
            if(master)  read(iuic) w2d
#ifdef SPMD
            call mp_scatter3d(w2d, v_srf, im,  jm, 1, beglat, endlat,  &
                 1, 1, 0, 0, 0)
            call fv_array_sync
#endif
            if(master)  call close_file(iuic)


        else ! not native format

            call mpp_open(iuic, restart_name, &
                 action=MPP_RDONLY, form=MPP_NETCDF, threading=MPP_SINGLE )
            if (master ) then
                call mpp_read_meta(iuic)
                call mpp_get_info(iuic,ndim,nvar,natt,ntime)

                allocate(attributes(natt))
                call mpp_get_atts(iuic,attributes)

                do n= 1,natt
!            select case (attributes(n)%name)
                   select case ( mpp_get_att_name(attributes(n)))
                   case ('nlon')
                       im= int(mpp_get_att_real_scalar(attributes(n)))
                   case ('mlat')
                       jm= int(mpp_get_att_real_scalar(attributes(n)))
                   end select
                enddo
                deallocate(attributes)
            endif

            call mpp_close(iuic)

            if ( im /= nlon .or. jm /= mlat ) then
                write(6,*) 'Dataset resolution:',im, jm
                call error_mesg('read_fv_rst:','resolution inconsistent in fv_srf_wnd', FATAL)
            endif
            call read_data ( restart_name, 'ps_bp', ps_bp(ids:ide,jcs:jce), domain = fv_domain)
            call read_data ( restart_name, 'U_srf', u_srf(ids:ide,jcs:jce), domain = fv_domain)
            call read_data ( restart_name, 'V_srf', v_srf(ids:ide,jcs:jce), domain = fv_domain)
        endif

    else

        if(master) then
            write(*,*) ' '
            write(*,*) '--------------------------------------------'
            write(*,*) 'File not found: ',restart_name
            write(*,*) 'Surface winds will be generated from D grid'
            write(*,*) '--------------------------------------------'
        endif


        do j=beglat,endlat
           do i=1,im
              ps_bp(i,j) = ps(i,j)
           enddo
        enddo

        call d2a3d(u(1,beglat-ng_d,km), v(1,beglat-ng_d,km), u_srf, v_srf, &
             nlon, mlat, 1, beglat, endlat, ng_d, ng_s, coslon, sinlon)
    endif

    call pmaxming('u_srf ', u_srf, im, jm, 1, beglat, endlat, 0,  0, 1.)
    call pmaxming('v_srf ', v_srf, im, jm, 1, beglat, endlat, 0,  0, 1.)

    if (use_native_format) deallocate ( w2d )


! check sum input data
    call print_check_sum ('Check sums for FV input data:')
    if( master ) then 
        call error_mesg ('read_fv_rst', 'FV netcdf restart file read.', NOTE)
    endif


#ifdef PLET
! Do Pure Lagrangian Error Test
    if ( nq >= 7 ) then
        p0 = (1.E5) ** kappa
!$omp  parallel do private(i,j,k)
        do k=ksp,kep
           do j=beglat, endlat
              do i=1,im
                 q(i,j,k,nq-2) = pt(i,j,k)                   
                 q(i,j,k,nq-1) = pt(i,j,k)/pkz(i,j,k) * p0
                 q(i,j,k,nq) = pkz(i,j,k) * p0
              enddo
           enddo
        enddo
        call fv_array_sync()
    endif
#endif

  end subroutine read_fv_rst


  subroutine add_tracers( q, im, jm, km, nq, nqrst, beglat, endlat, ng,  &
       age_tracer, age_time  )

    integer, intent(in):: im,jm,km,nq
    integer, intent(in):: nqrst
    integer, intent(in):: beglat, endlat
    integer, intent(in):: ng
    logical, intent(in):: age_tracer

    real, intent(inout):: age_time
    real, intent(inout):: q(im, beglat-ng:endlat+ng, km, nq)
! Local:
    integer i,j,k,iq

    call fv_array_check( LOC(q) )
    do iq=nqrst+1, nq
!$omp  parallel do private(i,j,k)
       do k=ksp,kep
          do j=beglat, endlat
             do i=1,im
                q(i,j,k,iq) = 0.
             enddo
          enddo
       enddo
    enddo
    call fv_array_sync()

    if ( age_tracer ) age_time = 0.

  end subroutine add_tracers

  subroutine fv_restart(days, seconds)

#include "fv_arrays.h"

    integer, intent(out):: days
    integer, intent(out):: seconds

    real, allocatable :: g_phis(:,:)
    real, allocatable :: qtmp(:,:,:)
    integer :: iuhs          ! Unit to read in terrain data for cold start
    integer :: i,j,k
    integer :: im, jm, km
    character(len=128) :: restart_name

#ifdef MARS_GCM 
    character (len=128) :: filename, fieldname
#endif MARS_GCM

#include "fv_point.inc"

    call fv_print_chksums( 'Entering  fv_restart' )
    if (cold_start) then

#ifdef SW_DYN

        if(nlev /= 1)  call error_mesg('FV_init:','Set nlev=1 for Shallow Water Dynamics', FATAL)

!-----------------------
! Shallow water dynamics
!-----------------------
        ptop  = 0.
        ak(1) = 0.
        ak(2) = 0.
        bk(1) = 0.
        bk(2) = 1.
        KS    = 0
#ifndef USE_LIMA
        n_spong = 0
#endif
        call init_sw_ic(icase, nlon, mlat, nlev, ncnst, beglat, endlat, ng_d, ng_s,  &
             phis, u, v, pt, delp, q, grav, radius, omega, f_d, ighost )
#else

! Start from a "dry & cold" isothermal atmosphere at rest

#ifdef MARS_GCM
        filename= 'INPUT/surf.res.nc' 
        fieldname= 'topo'

        if (file_exist( trim(filename) )) then
            mountain = .true.

           if(master) then
              allocate ( g_phis(nlon,mlat) )
              call read_data( trim(filename), trim(fieldname), g_phis, no_domain=.true. )
           endif

           call mp_scatter3d(g_phis, phis,  nlon, mlat,  1,  &
                             beglat, endlat, 1, 1, ighost, ighost, 0)

           call fv_array_sync
           if(master) deallocate ( g_phis )

            phis = phis * grav
        else
            mountain = .false.
            if(master) write(6,*) 'No surface dataset found; setting phis = 0'
            phis = 0.
        endif
#else

        !Modify the filename if it has an ensemble instance appendix
        call get_instance_filename('INPUT/surf.res',restart_name )
        restart_name=trim(restart_name)

        if (file_exist(restart_name)) then
            if(master) write(6,*) 'Reading surface dataset into phis from ', restart_name
            mountain = .true.
            !if(master) &
                 call mpp_open(iuhs, restart_name, action=MPP_RDONLY,        &
                 form=MPP_NATIVE, access=MPP_SEQUENTIAL)
#ifdef SPMD
            !if(master) then
                allocate ( g_phis(nlon,mlat) )
                read(iuhs) g_phis
            !endif
            call mp_scatter3d(g_phis, phis,  nlon, mlat,  1,  &
                 beglat, endlat, 1, 1, ighost, ighost, 0)
            call fv_array_sync
            !if(master) &
                 deallocate ( g_phis )
#else
            read(iuhs) phis
#endif
            if(master) close (iuhs)
            phis = phis * grav
        else
            mountain = .false.
            if(master) write(6,*) 'No surface dataset found; setting phis = 0'
            phis = 0.
        endif

#endif MARS_GCM

        call set_eta(nlev, ks, ptop, ak, bk)
        call init_dry_atm(mountain, kappa, grav, rdgas)

        if ( age_tracer ) then
            age_time = 0.
!$omp parallel do private (i, j, k)
!         do k=1,nlev
            do k=ksp,kep
               do j=beglat, endlat
                  do i=1,nlon
                     q(i,j,k,ncnst) = 0.
                  enddo
               enddo
            enddo
            call fv_array_sync()
        endif

#endif
    else
!---------------------------
! Read from the restart file
!---------------------------
        restart_name ='INPUT/fv_rst.res'
        !The call get_instance_filename() is inside the read_fv_rst subroutine 
        call read_fv_rst(im, jm, km, ncnst, restart_name, days, seconds, restart_format)

!   call p_var(nlon, mlat, nlev, beglat, endlat, ptop, delp, ps,           &
!              pe, peln,  pk,  pkz,  kappa, q, ng_d, ncnst, adjust_dry_mass )
    endif

#ifdef SPMD
! Ghost phis boundary N/S
    if ( ighost /= 0 ) then
        call mp_send4d_ns(nlon, mlat, 1, 1, beglat, endlat,       &
             1, 1, ighost, ighost, phis)
        call mp_recv4d_ns(nlon, mlat, 1, 1, beglat, endlat,       &
             1, 1, ighost, ighost, phis)
        call fv_array_sync()
    endif
#endif

!------------------------------------------------
! Safety check of the (a,b) of the eta coordinate
!------------------------------------------------

    call check_eta(nlev, ks, ptop, ak, bk)

! Generate winds on the mass grid for the physics/coupler:
    if ( use_tendency ) then
        call d2a3d(u, v, ua, va, nlon, mlat, nlev, beglat, endlat,  &
             ng_d, ng_s, coslon, sinlon)
    endif
    call fv_print_chksums( 'Exiting  fv_restart' )
  end subroutine fv_restart
! $Id: fv_restart.F90,v 17.0 2009/07/21 02:53:25 fms Exp $

  subroutine write_fv_rst( dflnm, days, seconds, grav, oform)

#include "fv_arrays.h"

! !INPUT PARAMETERS:
    integer,  intent(in):: days
    integer,  intent(in):: seconds
    real, intent(in):: grav
    character*(*), intent(in):: dflnm           ! dynamics restart filename
    character(len=*), intent(in) :: oform

! local
    integer :: iuic                 ! Unit number
    real, allocatable :: w2d(:,:), w3d(:,:,:), latu(:), phis_alloc(:,:), delp_alloc(:,:,:)
    real :: ginv, qmax, qmin
    real :: pfull(nlev), phalf(nlev+1)
    integer i, j, k, n

    character(len=128)           :: tracer_name, tracer_longname, tracer_units
    type(axistype)               :: axis_lon, axis_lat, axis_lev, axis_levh, &
                                    axis_lonv, axis_latu, axis_time
    type(fieldtype), allocatable :: tracer_fields(:)
    type(fieldtype)              :: u_field, v_field, t_field, phis_field, delp_field

    type(fieldtype)              :: ak_field, bk_field

    character(len=*), parameter :: fms_tracers_restart_file = &
         'RESTART/'//fms_tracers_file

    character(len=128) :: directory, restart_name, tr_name
    logical :: use_native_format
    integer :: ids,ide, jds,jde, ics,ice,jcs,jce
    integer :: nsphum, nliq_wat, nice_wat, ncld_amt, ntracers
    real :: dummy(1) =1.0

#ifdef MARS_GCM
!rjw     fv_tr_name is not used in this subroutine      
#else
    character(len=20) :: fv_tr_name(4)
    data fv_tr_name / "sphum","liq_wat","ice_wat","cld_amt" /
#endif MARS_GCM

#include "fv_point.inc"

! output format
    if (trim(uppercase(oform)) == 'NATIVE') then
        use_native_format = .true.
    else if (trim(uppercase(oform)) == 'NETCDF') then
        use_native_format = .false.
    else
        call error_mesg ('write_fv_rst', &
             'invalid value for argument oform='//trim(oform), FATAL)
    endif

! check sum output data
  call print_check_sum ('Check sums for FV output data:')

    if (use_native_format) then
        if ( master ) then
!!$         if (file_exist(dflnm)) then
!!$             call mpp_open ( iuic, dflnm, action=MPP_OVERWR,        &
!!$                             form=MPP_NATIVE, access=MPP_SEQUENTIAL )
!!$         else
!!$             call mpp_open ( iuic, dflnm, action=MPP_WRONLY,        &
!!$                             form=MPP_NATIVE, access=MPP_SEQUENTIAL )
!!$            
!!$         endif
           !Modify the filename if it has an ensemble instance appendix
            call get_instance_filename(dflnm,restart_name )
            restart_name=trim(restart_name)
            iuic = open_restart_file(restart_name, 'write')
            write(iuic) nlon, mlat, nlev, nt_phys
            write(iuic) days, seconds
            write(iuic) ak, bk, ks
            write(6,*) ' '
            write(6,*) 'FV restart file header written for days=', days, 'sec=',seconds
        endif

        call pmaxmin('PS', ps(1,beglat), qmin, qmax, nlon, endlat-beglat+1, 0.01)

!---------------------------
! Write surface geopotential
!---------------------------
#ifdef SPMD
        allocate ( w2d(nlon,mlat) )
        call mp_gather3d(phis, w2d, nlon, mlat, 1, beglat, endlat, 1, 1,  &
             ighost, ighost, 0)
        if ( master ) then
            write(iuic) w2d
            ginv = 1./ grav
            call maxmin_global('Surface height', w2d, nlon, 1, mlat, ginv)
        endif
#else
        write(iuic) phis
#endif

        allocate ( w3d(nlon,mlat,nlev) )

!---------------------------
! Write u-wind
!---------------------------

#ifdef SPMD
        call mp_gather3d(u, w3d, nlon, mlat, nlev, beglat, endlat, 1, nlev, ng_d, ng_s, 0)
        if ( master ) then
            write(iuic) w3d
            call maxmin_global('U', w3d, nlon, mlat, nlev, 1.)
        endif
#else
        write(iuic) u
#endif

!---------------------------
! Write v-wind
!---------------------------

#ifdef SPMD
        call mp_gather3d(v, w3d, nlon, mlat, nlev, beglat, endlat, 1, nlev, ng_d, ng_d, 0)
        if ( master ) then
            write(iuic) w3d
            call maxmin_global('V', w3d, nlon, mlat, nlev, 1.)
        endif
#else
        write(iuic) v
#endif

!------------------------------------
! Write virtual potential temperature
!------------------------------------

#ifdef SPMD
        call mp_gather3d(pt, w3d, nlon, mlat, nlev, beglat, endlat, 1, nlev, ng_d, ng_d, 0)
        if ( master ) then
            write(iuic) w3d
            call maxmin_global('T', w3d, nlon, mlat, nlev, 1.)
        endif
#else
        write(iuic) pt
#endif

!------------------------------------
! Write pressure thickness (pascal)
!------------------------------------

#ifdef SPMD
        call mp_gather3d(delp, w3d, nlon, mlat, nlev, beglat, endlat, 1, nlev, 0, 0, 0)
        if ( master ) then
            write(iuic) w3d
!            call maxmin_global('DELP', w3d, nlon, mlat, nlev, 1.)
        endif
#else
        write(iuic) delp
#endif


        if(ncnst /=  0) then 
!------------------------------------
! Write tracers one at a time
!------------------------------------

            do n=1,nt_phys 

#ifdef SPMD
               call mp_gather3d(q(1,beglat-ng_d,1,n), w3d, nlon, mlat, nlev,    &
                    beglat, endlat, 1, nlev, ng_d, ng_d, 0)
#else
               do k=1,nlev
                  do j=beglat,endlat
                     do i=1,nlon
                        w3d(i,j,k) = q(i,j,k,n)
                     enddo
                  enddo
               enddo
#endif
               if ( master ) then
                   write(iuic) w3d
                   call maxmin_global('Q', w3d, nlon, mlat, nlev, 1.)
               endif
            enddo
            if ( master .and. age_tracer ) write(iuic) age_time
        endif

        if( master ) then 
            call close_file(iuic)
            write(*,*) 'FV restart file written and closed'
        endif

!--------------------------------------------------------------
! Write surface winds on A grid for coupler to restart properly
!--------------------------------------------------------------

        if ( master ) then
           !Modify the filename if it has an ensemble instance appendix
            call get_instance_filename('RESTART/fv_srf_wnd.res',restart_name )
            restart_name=trim(restart_name)
            iuic = open_restart_file(restart_name, 'write')
            write(iuic) nlon, mlat
        endif

#ifdef SPMD
        call mp_gather3d(ps_bp, w2d, nlon, mlat, 1, beglat, endlat, 1, 1,  &
             0, 0, 0)
#endif
        if (master) write(iuic) w2d

#ifdef SPMD
        call mp_gather3d(u_srf, w2d, nlon, mlat, 1, beglat, endlat, 1, 1,  &
             0, 0, 0)
#endif
        if (master) write(iuic) w2d

#ifdef SPMD
        call mp_gather3d(v_srf, w2d, nlon, mlat, 1, beglat, endlat, 1, 1,  &
             0, 0, 0)
#endif
        if (master) then
            write(iuic) w2d
            call close_file(iuic)
        endif

!----------------------------------------
! Write additional tracers to netCDF file
!----------------------------------------

        if (do_fms_tracer_manager .and. ncnst > nt_phys) then !{

            if(master) &
                 & write(*,*)'Additional tracers present in field_table. These will be saved in a netcdf file.'

#ifdef MARS_GCM
            call get_eta_level(nlev, p_ref, pfull, phalf, 0.01)
#else
            call get_eta_level(nlev, 1.E5, pfull, phalf, 0.01)
#endif MARS_GCM

            !Modify the filename if it has an ensemble instance appendix
            call get_instance_filename(fms_tracers_restart_file,restart_name )
            restart_name=trim(restart_name)
            call mpp_open(iuic, restart_name, &
                 action=MPP_OVERWR, form=MPP_NETCDF, threading=MPP_SINGLE )

! define axes
            call mpp_write_meta(iuic, axis_lon, 'lon', 'degree_E', 'longitude', &
                 data=lon )
            call mpp_write_meta(iuic, axis_lat, 'lat', 'degree_N', 'latitude', &
                 data=lat )
            call mpp_write_meta(iuic, axis_lev, 'pfull', 'mb', 'ref full pressure level', &
                 data=pfull )

! define fields
            allocate(tracer_fields(ncnst))

            do n = nt_phys+1, ncnst
               call tr_get_tracer_names(MODEL_ATMOS, n, &
                    tracer_name, tracer_longname, tracer_units)
               call mpp_write_meta(iuic, tracer_fields(n), &
                    (/axis_lon, axis_lat, axis_lev/), &
                    tracer_name, tracer_units, tracer_longname, pack=1) ! want doubles
            enddo

! write axes
            call mpp_write(iuic, axis_lon)
            call mpp_write(iuic, axis_lat)
            call mpp_write(iuic, axis_lev)

! write tracers
            do n = nt_phys+1, ncnst
#ifdef SPMD
               call mp_gather3d(q(1,beglat-ng_d,1,n), w3d, nlon, mlat, nlev,    &
                    beglat, endlat, 1, nlev, ng_d, ng_d, 0)
#else
               do k=1,nlev
                  do j=beglat,endlat
                     do i=1,nlon
                        w3d(i,j,k) = q(i,j,k,n)
                     enddo
                  enddo
               enddo
#endif
               if(master) call mpp_write(iuic, tracer_fields(n), data=w3d)

            enddo

            call mpp_close(iuic)
            deallocate(tracer_fields)

        endif !}


        deallocate ( w3d )
        deallocate ( w2d )

    else ! ! ------------------NetCDF restart files---------
        call mpp_get_data_domain   ( fv_domain, ids, ide, jds, jde )
        call mpp_get_compute_domain( fv_domain, ics, ice, jcs, jce )
        restart_name = trim(dflnm) ! add netcdf file suffix

        !Modify the filename if it has an ensemble instance appendix
        call get_instance_filename(restart_name,restart_name )
        restart_name=trim(restart_name)
            
        call mpp_open(iuic, restart_name, &
             action=MPP_OVERWR, form=MPP_NETCDF, threading=MPP_SINGLE )

! define axes
#ifdef MARS_GCM
        call get_eta_level(nlev, p_ref, pfull, phalf, 0.01)
#else
        call get_eta_level(nlev, 1.E5, pfull, phalf, 0.01)
#endif MARS_GCM
!        if(master) write(*,*) "phalf values are ",phalf    
! define longitude axis for u and tracer fields
        call mpp_write_meta(iuic, axis_lon, 'lon', 'degree_E', 'longitude', &
             data=lon )

! define longitude axis for v field
        call mpp_write_meta(iuic, axis_lonv, 'lonv', 'degree_E', 'longitude', &
             data=lonb(1:size(lonb)-1) )


! define latitude axis for v and tracer fields
        call mpp_write_meta(iuic, axis_lat, 'lat', 'degree_N', 'latitude', &
             data=lat )

! define latitude axis for u field
        allocate(latu(size(lat)))
        ! The first latitude corresponds to meaningless data for the South Pole. 
        ! For interpolation reasons, it has been decided to set this latitude to 
        ! the second latitude with a small offset for Ferret compatibility.
        latu(1)  = latb(2) - 1e-5
        latu(2:) = latb(2:size(lat))
        call mpp_write_meta(iuic, axis_latu, 'latu', 'degree_N', 'latitude', &
             data=latu )

! define pressure axis
        call mpp_write_meta(iuic, axis_lev, 'pfull', 'mb', 'ref full pressure level', &
             data=pfull )
        call mpp_write_meta(iuic, axis_levh, 'phalf', 'mb', 'ref half pressure level', &
             data=phalf )

        call mpp_write_meta(iuic, axis_time, 'Time', 'time level', 'Time', &
                            cartesian='T', data=dummy)

        call mpp_write_meta(iuic, ak_field, (/axis_levh/), &
             'ak', 'Pa', 'ak-field', pack=1) ! want doubles
        call mpp_write_meta(iuic, bk_field, (/axis_levh/), &
             'bk', ' ', 'bk-field', pack=1)

        call mpp_write_meta(iuic, 0, 'nlon', ival=nlon)     
        call mpp_write_meta(iuic, 0, 'mlat', ival=mlat)     
        call mpp_write_meta(iuic, 0, 'nlev', ival=nlev)     

        call mpp_write_meta(iuic, 0, 'nt_phys', nt_phys)     
        call mpp_write_meta(iuic, 0, 'days', days)     
        call mpp_write_meta(iuic, 0, 'seconds', seconds)     
        call mpp_write_meta(iuic, 0, 'ks', ks)     

        call mpp_write_meta(iuic, u_field, (/axis_lon , axis_latu, axis_lev, axis_time/), &
             'U', 'm/s', 'u-field', pack=1)
        call mpp_write_meta(iuic, v_field, (/axis_lonv, axis_lat , axis_lev, axis_time/), &
             'V', 'm/s', 'v-field', pack=1)
        call mpp_write_meta(iuic, t_field, (/axis_lon , axis_lat , axis_lev, axis_time/), &
             'T', 'K', 'Temperature', pack=1)
        call mpp_write_meta(iuic, phis_field, (/axis_lon, axis_lat, axis_time/), &
             'Surface_geopotential', 'm**2/s**2', 'Surface Height', pack=1)
        call mpp_write_meta(iuic, delp_field, (/axis_lon, axis_lat, axis_lev, axis_time/), &
             'DELP', 'pa', 'delp-field', pack=1)

! write axes
        call mpp_write(iuic, axis_time)
        call mpp_write(iuic, axis_lon)
        call mpp_write(iuic, axis_lonv)
        call mpp_write(iuic, axis_lat)
        call mpp_write(iuic, axis_latu)
        call mpp_write(iuic, axis_lev)

        call mpp_write(iuic, bk_field, bk)     
        call mpp_write(iuic, ak_field, ak)     

        call mpp_write(iuic, u_field, fv_domain, u(ids:ide,jds:jde,:))
        call mpp_write(iuic, v_field, fv_domain, v(ids:ide,jds:jde,:))
        call mpp_write(iuic, t_field, fv_domain, pt(ids:ide,jds:jde,:))
! The mpp_io code uses Cray pointers which doesn't work well with the PSET 
! 2D arrays of the FV code. Therefore a copy of the 2D array is needed before 
! passing the array to the mpp_write routine.
        allocate (phis_alloc(ids:ide,jds:jde))
        phis_alloc = -999.9
        phis_alloc(ids:ide,jcs:jce) = phis(ids:ide,jcs:jce)
        call mpp_write(iuic, phis_field, fv_domain, phis_alloc)
        allocate (delp_alloc(ids:ide,jds:jde,nlev))
        delp_alloc = -999.9
        delp_alloc(ids:ide,jcs:jce,:) = delp(ids:ide,jcs:jce,:)   !RASF Bug fix for array bounds
        call mpp_write(iuic, delp_field, fv_domain, delp_alloc)   !RASF Bug fix for array bounds
        deallocate(delp_alloc)


        call mpp_close(iuic)


        call error_mesg ('write_fv_rst', 'Writing NetCDF formatted restart file: '//trim(restart_name), NOTE)

        call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers)


        restart_name = trim(fms_tracers_restart_file)!//'1.nc' ! add netcdf file suffix
        !Modify the filename if it has an ensemble instance appendix
        call get_instance_filename(restart_name,restart_name )
        restart_name=trim(restart_name)

        call mpp_open(iuic, restart_name, &
             action=MPP_OVERWR, form=MPP_NETCDF, threading=MPP_SINGLE )

! define axes

        call mpp_write_meta(iuic, axis_time, 'Time', 'time level', &
                            'Time',cartesian='T', data=dummy)

! define longitude axis for tracer fields
        call mpp_write_meta(iuic, axis_lon, 'lon', 'degree_E', 'longitude', &
             data=lon )

! define latitude axis for tracer fields
        call mpp_write_meta(iuic, axis_lat, 'lat', 'degree_N', 'latitude', &
             data=lat )

! define pressure axis
        call mpp_write_meta(iuic, axis_lev, 'pfull', 'mb', 'ref full pressure level', &
             data=pfull )

! define fields
        allocate(tracer_fields(ntracers))

        do n = 1,ntracers
           call tr_get_tracer_names(MODEL_ATMOS, n, &
                tracer_name, tracer_longname, tracer_units)
           call mpp_write_meta(iuic, tracer_fields(n), &
                (/axis_lon, axis_lat, axis_lev, axis_time/), &
                tracer_name, tracer_units, tracer_longname, pack=1)
        enddo

! write axes
        call mpp_write(iuic, axis_lon)
        call mpp_write(iuic, axis_lat)
        call mpp_write(iuic, axis_lev)
        call mpp_write(iuic, axis_time)

        do n = 1, ntracers
           call mpp_write(iuic, tracer_fields(n), fv_domain, q(ids:ide,jds:jde,:,n))
        enddo

        deallocate(tracer_fields)
        call mpp_close(iuic)


        !Modify the filename if it has an ensemble instance appendix
        call get_instance_filename('RESTART/fv_srf_wnd.res',restart_name )
        restart_name=trim(restart_name)
        call mpp_open(iuic, restart_name, &
             action=MPP_OVERWR, form=MPP_NETCDF, threading=MPP_SINGLE )

! define axes

! define longitude axis for u and tracer fields
        call mpp_write_meta(iuic, axis_lon, 'lon', 'degree_E', 'longitude', &
             data=lon )

! define longitude axis for v field
        call mpp_write_meta(iuic, axis_lonv, 'lonv', 'degree_E', 'longitude', &
             data=lonb(1:size(lonb)-1) )


! define latitude axis for v and tracer fields
        call mpp_write_meta(iuic, axis_lat, 'lat', 'degree_N', 'latitude', &
             data=lat )

! define latitude axis for u field
        latu(1)  = lat(1)
        latu(2:) = latb(2:size(lat))
        call mpp_write_meta(iuic, axis_latu, 'latu', 'degree_N', 'latitude', &
             data=latu )

        call mpp_write_meta(iuic, 'nlon', ival=nlon)     
        call mpp_write_meta(iuic, 'mlat', ival=mlat)     

        call mpp_write_meta(iuic, u_field, (/axis_lon , axis_latu/), &
             'U_srf', 'm/s', 'surface u-field', pack=1)
        call mpp_write_meta(iuic, v_field, (/axis_lonv, axis_lat /), &
             'V_srf', 'm/s', 'surface v-field', pack=1)
        call mpp_write_meta(iuic, phis_field, (/axis_lon, axis_lat/), &
             'ps_bp', 'pa', 'surface pressure', pack=1)

! write axes
        call mpp_write(iuic, axis_lon)
        call mpp_write(iuic, axis_lonv)
        call mpp_write(iuic, axis_lat)
        call mpp_write(iuic, axis_latu)
! The mpp_io code uses Cray pointers which doesn't work well with the PSET 
! 2D arrays of the FV code. Therefore a copy of the 2D array is needed before 
! passing the array to the mpp_write routine.
        phis_alloc = -999.9
        phis_alloc(ids:ide,jcs:jce) = ps_bp(ids:ide,jcs:jce)
        call mpp_write(iuic, phis_field, fv_domain, phis_alloc)
        phis_alloc(ids:ide,jcs:jce) = u_srf(ids:ide,jcs:jce)
        call mpp_write(iuic, u_field, fv_domain, phis_alloc)
        phis_alloc(ids:ide,jcs:jce) = v_srf(ids:ide,jcs:jce)
        call mpp_write(iuic, v_field, fv_domain, phis_alloc)




    endif ! end IF NETCDF
    call fv_array_sync()

  end subroutine write_fv_rst
!---------------------------------------------------------------------

  subroutine maxmin_global(qname, a, im, jm, km, fac)

    implicit none

    character*(*)  qname
    integer im, jm, km
    real a(im,jm, km)
    real fac
    real pmax(km), pmin(km)
    real qmax, qmin
    integer i, j, k

!$omp parallel do private(i, j, k)
!Balaji: not parallelized, called with unshared arrays
    do k=1, km
       pmax(k) = a(1,1,k)
       pmin(k) = a(1,1,k)
       do j=1,jm
          do i=1, im
             pmax(k) = max (pmax(k), a(i,j,k))
             pmin(k) = min (pmin(k), a(i,j,k))
          enddo
       enddo
    enddo

    qmax = pmax(1)
    qmin = pmin(1)
    do k=2,km
       qmax = max (qmax, pmax(k))
       qmin = min (qmin, pmin(k))
    enddo

    write(6,*) qname, ' max = ', qmax * fac, ' min = ', qmin * fac
  end subroutine maxmin_global


subroutine print_check_sum (label)

#include "fv_arrays.h"

 character(len=*),      intent(in) :: label
 integer :: n, chksum_unit
 integer :: isl, iel, jsl, jel
 integer(I8_KIND) :: zsum
#ifdef use_shared_pointers
#include "fv_point.inc"
#endif

  call mpp_get_compute_domain( fv_domain, isl, iel, jsl, jel )
  chksum_unit = stdout()

  if (mpp_pe() == mpp_root_pe()) write (chksum_unit,'(/a)') trim(label)
! sum all restart variables (except eta,peta)
  zsum = mpp_chksum(pt(isl:iel,jsl:jel,:))
  if (mpp_pe() == mpp_root_pe()) write (chksum_unit,10) '(t)',zsum
  zsum = mpp_chksum(u(isl:iel,jsl:jel,:))
  if (mpp_pe() == mpp_root_pe()) write (chksum_unit,10) '(u)',zsum
  zsum = mpp_chksum(v(isl:iel,jsl:jel,:))
  if (mpp_pe() == mpp_root_pe()) write (chksum_unit,10) '(v)',zsum
  zsum = 0
! sum of all tracers
  do n = 1, size(q,4)
    zsum = zsum + mpp_chksum(q(isl:iel,jsl:jel,:,n))
  enddo
  if (mpp_pe() == mpp_root_pe()) write (chksum_unit,10) '(r)',zsum
10 format ('chksum',a7,' = ',z16)

end subroutine print_check_sum

end module fv_restart_mod
