diff --git a/atmos_cubed_sphere b/atmos_cubed_sphere
index 3a4dfd8c6..9a290ee5f 160000
--- a/atmos_cubed_sphere
+++ b/atmos_cubed_sphere
@@ -1 +1 @@
-Subproject commit 3a4dfd8c6c4ceb8cec06397f25cb229ecd98065b
+Subproject commit 9a290ee5f3307a9d7b02762fe6c29a1f5ccd7f55
diff --git a/atmos_model.F90 b/atmos_model.F90
index 23e30e76c..9a73e0151 100644
--- a/atmos_model.F90
+++ b/atmos_model.F90
@@ -113,7 +113,8 @@ module atmos_model_mod
                               FV3GFS_diag_register, FV3GFS_diag_output,  &
                               DIAG_SIZE
 use fv_iau_mod,         only: iau_external_data_type,getiauforcing,iau_initialize
-use module_fv3_config,  only: output_1st_tstep_rst, first_kdt, nsout
+use module_fv3_config,  only: output_1st_tstep_rst, first_kdt, nsout,    &
+                              frestart, restart_endfcst
 
 !-----------------------------------------------------------------------
 
@@ -221,7 +222,8 @@ module atmos_model_mod
   logical,parameter :: flip_vc = .true.
 #endif
 
-  real(kind=IPD_kind_phys), parameter :: zero=0.0, one=1.0
+  real(kind=IPD_kind_phys), parameter :: zero = 0.0_IPD_kind_phys, &
+                                         one  = 1.0_IPD_kind_phys
 
 contains
 
@@ -944,7 +946,7 @@ end subroutine update_atmos_model_state
 subroutine atmos_model_end (Atmos)
   type (atmos_data_type), intent(inout) :: Atmos
 !---local variables
-  integer :: idx
+  integer :: idx, seconds
 #ifdef CCPP
   integer :: ierr
 #endif
@@ -952,9 +954,11 @@ subroutine atmos_model_end (Atmos)
 !-----------------------------------------------------------------------
 !---- termination routine for atmospheric model ----
                                               
-    call atmosphere_end (Atmos % Time, Atmos%grid)
-    call FV3GFS_restart_write (IPD_Data, IPD_Restart, Atm_block, &
-                               IPD_Control, Atmos%domain)
+    call atmosphere_end (Atmos % Time, Atmos%grid, restart_endfcst)
+    if(restart_endfcst) then
+      call FV3GFS_restart_write (IPD_Data, IPD_Restart, Atm_block, &
+                                 IPD_Control, Atmos%domain)
+    endif
 
 #ifdef CCPP
 !   Fast physics (from dynamics) are finalized in atmosphere_end above;
@@ -1457,6 +1461,24 @@ subroutine update_atmos_chemistry(state, rc)
         enddo
       enddo
 
+      ! -- zero out accumulated fields
+!$OMP parallel do default (none) &
+!$OMP             shared  (nj, ni, Atm_block, IPD_Control, IPD_Data) &
+!$OMP             private (j, jb, i, ib, nb, ix)
+      do j = 1, nj
+        jb = j + Atm_block%jsc - 1
+        do i = 1, ni
+          ib = i + Atm_block%isc - 1
+          nb = Atm_block%blkno(ib,jb)
+          ix = Atm_block%ixp(ib,jb)
+          IPD_Data(nb)%coupling%rainc_cpl(ix)  = zero
+          if (.not.IPD_Control%cplflx) then
+            IPD_Data(nb)%coupling%rain_cpl(ix) = zero
+            IPD_Data(nb)%coupling%snow_cpl(ix) = zero
+          end if
+        enddo
+      enddo
+
       if (IPD_Control%debug) then
         ! -- diagnostics
         write(6,'("update_atmos: prsi   - min/max/avg",3g16.6)') minval(prsi),   maxval(prsi),   sum(prsi)/size(prsi)
@@ -1698,7 +1720,7 @@ subroutine assign_importdata(rc)
                   IPD_Data(nb)%Coupling%ficein_cpl(ix)   = zero
                   if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) > zero) then
                     if (datar8(i,j) >= IPD_control%min_seaice*IPD_Data(nb)%Sfcprop%oceanfrac(ix)) then
-                      IPD_Data(nb)%Coupling%ficein_cpl(ix) = datar8(i,j)
+                      IPD_Data(nb)%Coupling%ficein_cpl(ix) = max(zero, min(datar8(i,j),one))
 !                     if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) == one) IPD_Data(nb)%Sfcprop%slmsk(ix) = 2. !slmsk=2 crashes in gcycle on partial land points
                       IPD_Data(nb)%Sfcprop%slmsk(ix)         = 2.                                        !slmsk=2 crashes in gcycle on partial land points
                       IPD_Data(nb)%Coupling%slimskin_cpl(ix) = 4.
diff --git a/ccpp/framework b/ccpp/framework
index 42596226f..1c41a614d 160000
--- a/ccpp/framework
+++ b/ccpp/framework
@@ -1 +1 @@
-Subproject commit 42596226fcbbf4c151829ec99ff0649adc09a5df
+Subproject commit 1c41a614d5cc68a33e418c9c8aaac9d7a2b971e7
diff --git a/ccpp/physics b/ccpp/physics
index 02812f6b2..01ed01fb0 160000
--- a/ccpp/physics
+++ b/ccpp/physics
@@ -1 +1 @@
-Subproject commit 02812f6b250fe637e37089259c051d2018ad43c3
+Subproject commit 01ed01fb0b3112e96eb619e0339d88fb0201982f
diff --git a/fv3_cap.F90 b/fv3_cap.F90
index d04abdce2..a099ae32e 100644
--- a/fv3_cap.F90
+++ b/fv3_cap.F90
@@ -30,7 +30,7 @@ module fv3gfs_cap_mod
                                     calendar, calendar_type, cpl,            &
                                     force_date_from_configure,               &
                                     cplprint_flag,output_1st_tstep_rst,      &
-                                    first_kdt                            
+                                    first_kdt,num_restart_interval       
 
   use module_fv3_io_def,      only: num_pes_fcst,write_groups,app_domain,    &
                                     num_files, filename_base,                &
@@ -278,9 +278,16 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
     CALL ESMF_ConfigLoadFile(config=CF ,filename='model_configure' ,rc=RC)
     if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
 !
-    CALL ESMF_ConfigGetAttribute(config=CF,value=restart_interval, &
-                                 label ='restart_interval:',rc=rc)
+    num_restart_interval = ESMF_ConfigGetLen(config=CF, label ='restart_interval:',rc=rc)
     if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+    if(mype == 0) print *,'af nems config,num_restart_interval=',num_restart_interval
+    if (num_restart_interval<=0) num_restart_interval = 1
+    allocate(restart_interval(num_restart_interval))
+    restart_interval = 0
+    CALL  ESMF_ConfigGetAttribute(CF,valueList=restart_interval,label='restart_interval:', &
+      count=num_restart_interval, rc=RC)
+    if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+    if(mype == 0) print *,'af nems config,restart_interval=',restart_interval
 !
     CALL ESMF_ConfigGetAttribute(config=CF,value=calendar, &
                                  label ='calendar:',rc=rc)
@@ -326,9 +333,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
                                    label ='app_domain:',rc=rc)
       if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
 
-      if(mype == 0) print *,'af nems config,restart_interval=',restart_interval, &
-      'quilting=',quilting,'write_groups=',write_groups,wrttasks_per_group,      &
-      'calendar=',trim(calendar),'calendar_type=',calendar_type
+      if(mype == 0) print *,'af nems config,quilting=',quilting,'write_groups=', &
+        write_groups,wrttasks_per_group,'calendar=',trim(calendar),'calendar_type=',calendar_type
 !
       CALL ESMF_ConfigGetAttribute(config=CF,value=num_files, &
                                    label ='num_files:',rc=rc)
diff --git a/gfsphysics/GFS_layer/GFS_physics_driver.F90 b/gfsphysics/GFS_layer/GFS_physics_driver.F90
index 4e34d5b39..36266aab8 100644
--- a/gfsphysics/GFS_layer/GFS_physics_driver.F90
+++ b/gfsphysics/GFS_layer/GFS_physics_driver.F90
@@ -1103,7 +1103,7 @@ subroutine GFS_physics_driver                                 &
 !*## CCPP ##
       enddo
 !
-!## CCPP ##* note: this block is not yet in CCPP
+!## CCPP ##* GFS_surface_generic.F90/GFS_surface_generic_pre_run
       if (Model%cplflx) then
         do i=1,im
           islmsk_cice(i) = nint(Coupling%slimskin_cpl(i))
@@ -1273,7 +1273,7 @@ subroutine GFS_physics_driver                                 &
           dtdt(i,k)  = zero
           dtdtc(i,k) = zero
 
-!## CCPP ##* note: this block is not yet in CCPP
+!## CCPP ##* GFS_typedefs.F90/interstitial_phys_reset
 !vay-2018
 ! Pure tendency arrays w/o accumulation of Phys-tendencies from each
 !      chain of GFS-physics (later add container for species)
@@ -1911,14 +1911,16 @@ subroutine GFS_physics_driver                                 &
 !    &,' stsoil=',stsoil(ipr,:)
 
 !  --- ...  surface energy balance over seaice
-!## CCPP ##* This section is not in the CCPP yet
+!## CCPP ##* sfc_sice.f/sfc_sice_run (local adjustment to avoid resetting islmsk after call to sfc_sice_run)
         if (Model%cplflx) then
           do i=1,im
             if (flag_cice(i)) then
                islmsk (i) = islmsk_cice(i)
             endif
           enddo
+!*## CCPP ##
 
+!## CCPP ##* sfc_cice.f/sfc_cice_run
 ! call sfc_cice for sea ice points in the coupled model (i.e. islmsk=4)
 !
           call sfc_cice                                                  &
@@ -1954,7 +1956,7 @@ subroutine GFS_physics_driver                                 &
             snowd3(:,2), qss3(:,2), snowmt, gflx3(:,2), cmm3(:,2), chh3(:,2),    &
             evap3(:,2),  hflx3(:,2))
 !*## CCPP ##
-!## CCPP ##* This section is not in the CCPP yet.
+!## CCPP ##* This section is not needed for CCPP.
         if (Model%cplflx) then
           do i = 1, im
             if (flag_cice(i)) then
@@ -2805,7 +2807,7 @@ subroutine GFS_physics_driver                                 &
       endif
 !*## CCPP ##
 
-!## CCPP ##* This block is not yet in CCPP
+!## CCPP ##* GFS_PBL_generic.F90/GFS_PBL_generic_post_run
       if (Model%cplchm) then
         do i = 1, im
           tem1 = max(Diag%q1(i), 1.e-8)
@@ -2814,7 +2816,6 @@ subroutine GFS_physics_driver                                 &
         enddo
         Coupling%dkt     (:,:) = dkt (:,:)
       endif
-!*## CCPP ##
 
 !     if (lprnt) then
 !       write(0,*) ' dusfc1=',dusfc1(ipr),' kdt=',kdt,' lat=',lat
@@ -2827,8 +2828,6 @@ subroutine GFS_physics_driver                                 &
 
 !  --- ...  coupling insertion
 
-!## CCPP ## This block is not in the CCPP yet. It should probably be put in 
-! GFS_PBL_generic.F90/GFS_PBL_generic_post_run.
       if (Model%cplflx) then
         do i=1,im
           if (Sfcprop%oceanfrac(i) > zero) then               ! Ocean only, NO LAKES
@@ -3182,7 +3181,7 @@ subroutine GFS_physics_driver                                 &
       Stateout%gq0(1:im,:,:) = Statein%qgrs(1:im,:,:) + dqdt(1:im,:,:) * dtp
 !*## CCPP ##
 
-! DH* TODO - WHERE IS THIS IN CCPP?
+!## CCPP ##* This is not in the CCPP yet.
 !================================================================================
 !     above: updates of the state by UGWP oro-GWS and RF-damp
 !  Diag%tav_ugwp & Diag%uav_ugwp(i,k)-Updated U-T state before moist/micro !  physics
@@ -3197,7 +3196,7 @@ subroutine GFS_physics_driver                                 &
           enddo
         enddo
       endif
-! *DH
+!*## CCPP ##
 
 !================================================================================
 ! It is not clear Do we need it, "ideaca_up", having stability check inside UGWP-module
@@ -3308,9 +3307,13 @@ subroutine GFS_physics_driver                                 &
         dtdt(1:im,:) = Stateout%gt0(1:im,:)
       endif   ! end if_ldiag3d/cnvgwd
 
-      if (Model%ldiag3d) then
+      if (Model%ldiag3d .or. Model%cplchm) then
         dqdt(1:im,:,1) = Stateout%gq0(1:im,:,1)
-      endif   ! end if_ldiag3d
+      endif   ! end if_ldiag3d/cplchm
+
+      if (Model%cplchm) then
+        Coupling%dqdti(1:im,:) = zero
+      endif   ! end if_cplchm
 !*## CCPP ##
 
 !## CCPP ## Only get_prs_fv3.F90/get_phi_fv3_run is a scheme (GFS_HYDRO is assumed to be undefined)
diff --git a/gfsphysics/GFS_layer/GFS_radiation_driver.F90 b/gfsphysics/GFS_layer/GFS_radiation_driver.F90
index bacf23d6c..8df014231 100644
--- a/gfsphysics/GFS_layer/GFS_radiation_driver.F90
+++ b/gfsphysics/GFS_layer/GFS_radiation_driver.F90
@@ -2135,18 +2135,40 @@ subroutine GFS_radiation_driver                             &
               Diag%fluxr(i,11-j) = Diag%fluxr(i,11-j) + tem0d * Statein%prsi(i,itop+kt)
               Diag%fluxr(i,14-j) = Diag%fluxr(i,14-j) + tem0d * Statein%prsi(i,ibtc+kb)
               Diag%fluxr(i,17-j) = Diag%fluxr(i,17-j) + tem0d * Statein%tgrs(i,itop)
+            enddo
+          enddo
 
 !       Anning adds optical depth and emissivity output
-              tem1 = 0.
-              tem2 = 0.
-              do k=ibtc,itop
-                tem1 = tem1 + cldtausw(i,k)      ! approx .55 mu channel
-                tem2 = tem2 + cldtaulw(i,k)      ! approx 10. mu channel
+          if (Model%lsswr .and. (nday > 0)) then
+            do j = 1, 3
+              do i = 1, IM
+                tem0d = raddt * cldsa(i,j)
+                itop  = mtopa(i,j) - kd
+                ibtc  = mbota(i,j) - kd
+                tem1 = 0.
+                do k=ibtc,itop
+                  tem1 = tem1 + cldtausw(i,k)      ! approx .55 um channel
+                enddo
+                Diag%fluxr(i,43-j) = Diag%fluxr(i,43-j) + tem0d * tem1
               enddo
-              Diag%fluxr(i,43-j) = Diag%fluxr(i,43-j) + tem0d * tem1
-              Diag%fluxr(i,46-j) = Diag%fluxr(i,46-j) + tem0d * (1.0-exp(-tem2))
             enddo
-          enddo
+          endif
+
+          if (Model%lslwr) then
+            do j = 1, 3
+              do i = 1, IM
+                tem0d = raddt * cldsa(i,j)
+                itop  = mtopa(i,j) - kd
+                ibtc  = mbota(i,j) - kd
+                tem2 = 0.
+                do k=ibtc,itop
+                  tem2 = tem2 + cldtaulw(i,k)      ! approx 10. um channel
+                enddo
+                Diag%fluxr(i,46-j) = Diag%fluxr(i,46-j) + tem0d * (1.0-exp(-tem2))
+              enddo
+            enddo
+          endif
+
         endif
 
       endif                                ! end_if_lssav
diff --git a/gfsphysics/physics/gfdl_cloud_microphys.F90 b/gfsphysics/physics/gfdl_cloud_microphys.F90
index ba4c814d6..f01486db0 100644
--- a/gfsphysics/physics/gfdl_cloud_microphys.F90
+++ b/gfsphysics/physics/gfdl_cloud_microphys.F90
@@ -3266,7 +3266,7 @@ subroutine fall_speed (ktop, kbot, den, qs, qi, qg, ql, tk, vts, vti, vtg)
             else
                 tc (k) = tk (k) - tice
                 vti (k) = (3. + log10 (qi (k) * den (k))) * (tc (k) * (aa * tc (k) + bb) + cc) + dd * tc (k) + ee
-                vti (k) = vi0 * exp (log_10 * vti (k)) * 0.8 
+                vti (k) = vi0 * exp (log_10 * vti (k)) * 0.9
                 vti (k) = min (vi_max, max (vf_min, vti (k)))
             endif
         enddo
diff --git a/gfsphysics/physics/module_sf_noahmp_glacier.f90 b/gfsphysics/physics/module_sf_noahmp_glacier.f90
old mode 100755
new mode 100644
diff --git a/gfsphysics/physics/module_sf_noahmplsm.f90 b/gfsphysics/physics/module_sf_noahmplsm.f90
old mode 100755
new mode 100644
diff --git a/gfsphysics/physics/module_wrf_utl.f90 b/gfsphysics/physics/module_wrf_utl.f90
old mode 100755
new mode 100644
diff --git a/gfsphysics/physics/moninedmf_hafs.f b/gfsphysics/physics/moninedmf_hafs.f
index 42de545f2..a035ad3d6 100644
--- a/gfsphysics/physics/moninedmf_hafs.f
+++ b/gfsphysics/physics/moninedmf_hafs.f
@@ -1360,7 +1360,11 @@ subroutine moninedmf_hafs(ix,im,km,ntrac,ntcw,dv,du,tau,rtg,      &
           tem = 0.5 * (diss(i,k-1)+diss(i,k))
           tem  = max(tem, 0.)
           ttend = tem / cp
-          tau(i,k) = tau(i,k) + 0.5*ttend
+          if (alpha .gt. 0.0) then
+            tau(i,k) = tau(i,k) + 0.5*ttend
+          else
+            tau(i,k) = tau(i,k) + 0.7*ttend  ! in HWRF/HMON, use 0.7
+          endif
         enddo
       enddo
 !
diff --git a/gfsphysics/physics/noahmp_tables.f90 b/gfsphysics/physics/noahmp_tables.f90
old mode 100755
new mode 100644
diff --git a/gfsphysics/physics/samfdeepcnv.f b/gfsphysics/physics/samfdeepcnv.f
index 76204ebb4..d29410d2f 100644
--- a/gfsphysics/physics/samfdeepcnv.f
+++ b/gfsphysics/physics/samfdeepcnv.f
@@ -1547,22 +1547,22 @@ subroutine samfdeepcnv(im,ix,km,delt,itc,ntc,ntk,ntr,delp,
       enddo
       enddo
       do i = 1, im
-        betamn = betas
-        if(islimsk(i) == 1) betamn = betal
-        if(ntk > 0) then
-          betamx = betamn + dbeta
-          if(tkemean(i) > tkemx) then
-            beta = betamn
-          else if(tkemean(i) < tkemn) then
-            beta = betamx
+        if(cnvflg(i)) then
+          betamn = betas
+          if(islimsk(i) == 1) betamn = betal
+          if(ntk > 0) then
+            betamx = betamn + dbeta
+            if(tkemean(i) > tkemx) then
+              beta = betamn
+            else if(tkemean(i) < tkemn) then
+              beta = betamx
+            else
+              tem = (betamx - betamn) * (tkemean(i) - tkemn)
+              beta = betamx - tem  / dtke
+            endif
           else
-            tem = (betamx - betamn) * (tkemean(i) - tkemn)
-            beta = betamx - tem  / dtke
+            beta = betamn
           endif
-        else
-          beta = betamn
-        endif
-        if(cnvflg(i)) then
           dz  = (sumx(i)+zi(i,1))/float(kbcon(i))
           tem = 1./float(kbcon(i))
           xlamd(i) = (1.-beta**tem)/dz
diff --git a/gfsphysics/physics/satmedmfvdifq.f b/gfsphysics/physics/satmedmfvdifq.f
index 11c047fd0..1cc0bbe89 100644
--- a/gfsphysics/physics/satmedmfvdifq.f
+++ b/gfsphysics/physics/satmedmfvdifq.f
@@ -148,7 +148,8 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke,
      &                     epsi,    beta,   chx,    cqx,
      &                     rdt,     rdz,    qmin,   qlmin,
      &                     rimin,   rbcr,   rbint,  tdzmin,
-     &                     rlmn,    rlmn1,  rlmx,   elmx,
+     &                     rlmn,    rlmn1,  rlmn2,  
+     &                     rlmx,    elmx,
      &                     ttend,   utend,  vtend,  qtend,
      &                     zfac,    zfmin,  vk,     spdk2,
      &                     tkmin,   tkminx, xkzinv, xkgdx,
@@ -172,7 +173,8 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke,
       parameter(gamcrt=3.,gamcrq=0.,sfcfrac=0.1)
       parameter(vk=0.4,rimin=-100.)
       parameter(rbcr=0.25,zolcru=-0.02,tdzmin=1.e-3)
-      parameter(rlmn=30.,rlmn1=5.,rlmx=300.,elmx=300.)
+      parameter(rlmn=30.,rlmn1=5.,rlmn2=10.)
+      parameter(rlmx=300.,elmx=300.)
       parameter(prmin=0.25,prmax=4.0)
       parameter(pr0=1.0,prtke=1.0,prscu=0.67)
       parameter(f0=1.e-4,crbmin=0.15,crbmax=0.35)
@@ -698,8 +700,9 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke,
 !         if(tem1 > 1.e-5) then
           tem1 = tvx(i,k+1)-tvx(i,k)
           if(tem1 > 0.) then
-             xkzo(i,k)  = min(xkzo(i,k),xkzinv)
-             xkzmo(i,k) = min(xkzmo(i,k),xkzinv)
+             xkzo(i,k)  = min(xkzo(i,k), xkzinv)
+             xkzmo(i,k) = min(xkzmo(i,k), xkzinv)
+             rlmnz(i,k) = min(rlmnz(i,k), rlmn2)
           endif
         enddo
       enddo
diff --git a/gfsphysics/physics/sfc_noahmp_drv.f b/gfsphysics/physics/sfc_noahmp_drv.f
old mode 100755
new mode 100644
diff --git a/gfsphysics/physics/sflx.f b/gfsphysics/physics/sflx.f
index 29467fe37..f84006be9 100644
--- a/gfsphysics/physics/sflx.f
+++ b/gfsphysics/physics/sflx.f
@@ -251,6 +251,7 @@ subroutine sflx                                                   &
       runoff2 = 0.0
       runoff3 = 0.0
       snomlt  = 0.0
+      rc      = 0.0
 
 !  --- ...  define local variable ice to achieve:
 !             sea-ice case,          ice =  1
diff --git a/io/FV3GFS_io.F90 b/io/FV3GFS_io.F90
index 7ba17b379..972d43fde 100644
--- a/io/FV3GFS_io.F90
+++ b/io/FV3GFS_io.F90
@@ -157,10 +157,10 @@ subroutine FV3GFS_restart_write (IPD_Data, IPD_Restart, Atm_block, Model, fv_dom
     type(domain2d),              intent(in)    :: fv_domain
     character(len=32), optional, intent(in)    :: timestamp
  
-    !--- read in surface data from chgres 
+    !--- write surface data from chgres 
     call sfc_prop_restart_write (IPD_Data%Sfcprop, Atm_block, Model, fv_domain, timestamp)
  
-    !--- read in physics restart data
+    !--- write physics restart data
     call phys_restart_write (IPD_Restart, Atm_block, Model, fv_domain, timestamp)
 
   end subroutine FV3GFS_restart_write
@@ -285,7 +285,7 @@ subroutine FV3GFS_IPD_checksum (Model, IPD_Data, Atm_block)
        temp2d(i,j,55) = IPD_Data(nb)%Coupling%visbmui(ix)
        temp2d(i,j,56) = IPD_Data(nb)%Coupling%visdfui(ix)
        temp2d(i,j,57) = IPD_Data(nb)%Coupling%sfcdsw(ix)
-       temp2d(i,j,59) = IPD_Data(nb)%Coupling%sfcnsw(ix)
+       temp2d(i,j,58) = IPD_Data(nb)%Coupling%sfcnsw(ix)
        temp2d(i,j,59) = IPD_Data(nb)%Coupling%sfcdlw(ix)
        temp2d(i,j,60) = IPD_Data(nb)%Grid%xlon(ix)
        temp2d(i,j,61) = IPD_Data(nb)%Grid%xlat(ix)
diff --git a/io/module_write_nemsio.F90 b/io/module_write_nemsio.F90
index 3afd66789..e51f64a52 100644
--- a/io/module_write_nemsio.F90
+++ b/io/module_write_nemsio.F90
@@ -51,7 +51,7 @@ subroutine nemsio_first_call(fieldbundle, imo, jmo, &
     integer, intent(in)                    :: wrt_mype, wrt_ntasks, wrt_mpi_comm
     integer, intent(in)                    :: wrt_nbdl, mybdl
     integer, intent(in)                    :: inidate(7)
-    real, intent(in)                       :: lat(:), lon(:)
+    real(8), intent(in)                    :: lat(:), lon(:)
     integer, optional,intent(out)          :: rc
 
 !** local vars
diff --git a/io/module_write_netcdf.F90 b/io/module_write_netcdf.F90
index 1fce3d8b9..fc59c75c9 100644
--- a/io/module_write_netcdf.F90
+++ b/io/module_write_netcdf.F90
@@ -42,7 +42,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
     real(4), dimension(:,:,:), allocatable :: arrayr4_3d,arrayr4_3d_save
     real(8), dimension(:,:,:), allocatable :: arrayr8_3d
 
-    real(8) rad2dg,x(im),y(jm)
+    real(8) x(im),y(jm)
     integer :: fieldCount, fieldDimCount, gridDimCount
     integer, dimension(:), allocatable   :: ungriddedLBound, ungriddedUBound
 
@@ -56,7 +56,8 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
     character(len=ESMF_MAXSTR) :: attName, fldName
 
     integer :: varival
-    real(4) :: varr4val, scale_fact, compress_err, offset, dataMin, dataMax
+    real(4) :: varr4val, scale_fact, offset, dataMin, dataMax
+    real(4), allocatable, dimension(:) :: compress_err
     real(8) :: varr8val
     character(len=ESMF_MAXSTR) :: varcval
 
@@ -71,10 +72,10 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
 !
 !!
 !
-    rad2dg = 45./atan(1.0)
 
     call ESMF_FieldBundleGet(fieldbundle, fieldCount=fieldCount, rc=rc); ESMF_ERR_RETURN(rc)
 
+    allocate(compress_err(fieldCount)); compress_err=-999.
     allocate(fldlev(fieldCount)) ; fldlev = 0
     allocate(fcstField(fieldCount))
     allocate(varids(fieldCount))
@@ -117,13 +118,13 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
     if (mype==0) then
 
     if (ideflate == 0) then
-        ncerr = nf90_create(trim(filename), cmode=IOR(NF90_CLOBBER,NF90_64BIT_OFFSET), &
+        ncerr = nf90_create(trim(filename), cmode=IOR(IOR(NF90_CLOBBER,NF90_64BIT_OFFSET),NF90_SHARE), &
         ncid=ncid); NC_ERR_STOP(ncerr)
         ncerr = nf90_set_fill(ncid, NF90_NOFILL, oldMode); NC_ERR_STOP(ncerr)
     else
         ncerr = nf90_create(trim(filename), cmode=IOR(IOR(NF90_CLOBBER,NF90_NETCDF4),NF90_CLASSIC_MODEL), &
         ncid=ncid); NC_ERR_STOP(ncerr)
-        ! if compression on use HDF5 format with default _FillValue
+        ncerr = nf90_set_fill(ncid, NF90_NOFILL, oldMode); NC_ERR_STOP(ncerr)
     endif
 
     ! define dimensions
@@ -156,28 +157,32 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
       ! define variables
       if (fldlev(i) == 1) then
         if (typekind == ESMF_TYPEKIND_R4) then
-          ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, &
-                               (/im_dimid,jm_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr)
           if (ideflate > 0) then
-             ! shuffle filter on for lossless compression
-             ncerr = nf90_def_var_deflate(ncid, varids(i), 1, 1, ideflate)
-             NC_ERR_STOP(ncerr)
+            ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, &
+                                (/im_dimid,jm_dimid,time_dimid/), varids(i), &
+                                 shuffle=.true.,deflate_level=ideflate, &
+                                 chunksizes=(/im,jm,1/)); NC_ERR_STOP(ncerr)
+          else
+            ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, &
+                                (/im_dimid,jm_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr)
           endif
         else if (typekind == ESMF_TYPEKIND_R8) then
           ncerr = nf90_def_var(ncid, trim(fldName), NF90_DOUBLE, &
-                               (/im_dimid,jm_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr)
+                              (/im_dimid,jm_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr)
         else
           write(0,*)'Unsupported typekind ', typekind
           stop
         end if
       else if (fldlev(i) > 1) then
          if (typekind == ESMF_TYPEKIND_R4) then
-           ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, &
-                                (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr)
            if (ideflate > 0) then
-              ! shuffle filter off since lossy compression used
-              ncerr = nf90_def_var_deflate(ncid, varids(i), 0, 1, ideflate)
-              NC_ERR_STOP(ncerr)
+             ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, &
+                                 (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i), &
+                                 shuffle=.false.,deflate_level=ideflate, &
+                                 chunksizes=(/im,jm,1,1/)); NC_ERR_STOP(ncerr)
+           else
+             ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, &
+                                 (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr)
            endif
          else if (typekind == ESMF_TYPEKIND_R8) then
            ncerr = nf90_def_var(ncid, trim(fldName), NF90_DOUBLE, &
@@ -219,8 +224,8 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
            call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", &
                                   name=trim(attName), value=varr8val, &
                                   rc=rc); ESMF_ERR_RETURN(rc)
-           if (trim(attName) /= '_FillValue' .or. ideflate == 0) then
-              ! FIXME:  _FillValue must be cast to var type when using NF90_NETCDF4
+           if (trim(attName) /= '_FillValue' ) then
+              ! FIXME:  _FillValue must be cast to var type for recent versions of netcdf
               ncerr = nf90_put_att(ncid, varids(i), trim(attName), varr8val); NC_ERR_STOP(ncerr)
            endif
 
@@ -236,6 +241,25 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
 
     end do   ! i=1,fieldCount
 
+! write grid_xt, grid_yt attributes
+    if (trim(output_grid) == 'gaussian_grid' .or. &
+        trim(output_grid) == 'regional_latlon') then
+      ncerr = nf90_put_att(ncid, im_varid, "long_name", "T-cell longitude"); NC_ERR_STOP(ncerr)
+      ncerr = nf90_put_att(ncid, im_varid, "units", "degrees_E"); NC_ERR_STOP(ncerr)
+      ncerr = nf90_put_att(ncid, jm_varid, "long_name", "T-cell latiitude"); NC_ERR_STOP(ncerr)
+      ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees_N"); NC_ERR_STOP(ncerr)
+    else if (trim(output_grid) == 'rotated_latlon') then
+      ncerr = nf90_put_att(ncid, im_varid, "long_name", "rotated T-cell longiitude"); NC_ERR_STOP(ncerr)
+      ncerr = nf90_put_att(ncid, im_varid, "units", "degrees"); NC_ERR_STOP(ncerr)
+      ncerr = nf90_put_att(ncid, jm_varid, "long_name", "rotated T-cell latiitude"); NC_ERR_STOP(ncerr)
+      ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees"); NC_ERR_STOP(ncerr)
+    else if (trim(output_grid) == 'lambert_conformal') then
+      ncerr = nf90_put_att(ncid, im_varid, "long_name", "x-coordinate of projection"); NC_ERR_STOP(ncerr)
+      ncerr = nf90_put_att(ncid, im_varid, "units", "meters"); NC_ERR_STOP(ncerr)
+      ncerr = nf90_put_att(ncid, jm_varid, "long_name", "y-coordinate of projection"); NC_ERR_STOP(ncerr)
+      ncerr = nf90_put_att(ncid, jm_varid, "units", "meters"); NC_ERR_STOP(ncerr)
+    endif
+
     ncerr = nf90_enddef(ncid); NC_ERR_STOP(ncerr)
     end if
 
@@ -247,31 +271,19 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
     if (mype==0) then
        if (trim(output_grid) == 'gaussian_grid' .or. &
            trim(output_grid) == 'regional_latlon') then
-          ncerr = nf90_put_var(ncid, im_varid, values=rad2dg*arrayr8(:,1)  ); NC_ERR_STOP(ncerr)
-          ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr)
-          ncerr = nf90_put_att(ncid, im_varid, "long_name", "T-cell longitude"); NC_ERR_STOP(ncerr)
-          ncerr = nf90_put_att(ncid, im_varid, "units", "degrees_E"); NC_ERR_STOP(ncerr)
-          ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr)
+          ncerr = nf90_put_var(ncid, im_varid, values=arrayr8(:,1)  ); NC_ERR_STOP(ncerr)
        else if (trim(output_grid) == 'rotated_latlon') then
           do i=1,im
              x(i) = lon1 + (lon2-lon1)/(im-1) * (i-1)
           enddo
           ncerr = nf90_put_var(ncid, im_varid, values=x  ); NC_ERR_STOP(ncerr)
-          ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr)
-          ncerr = nf90_put_att(ncid, im_varid, "long_name", "rotated T-cell longiitude"); NC_ERR_STOP(ncerr)
-          ncerr = nf90_put_att(ncid, im_varid, "units", "degrees"); NC_ERR_STOP(ncerr)
-          ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr)
        else if (trim(output_grid) == 'lambert_conformal') then
           do i=1,im
              x(i) = dx * (i-1)
           enddo
           ncerr = nf90_put_var(ncid, im_varid, values=x  ); NC_ERR_STOP(ncerr)
-          ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr)
-          ncerr = nf90_put_att(ncid, im_varid, "long_name", "x-coordinate of projection"); NC_ERR_STOP(ncerr)
-          ncerr = nf90_put_att(ncid, im_varid, "units", "meters"); NC_ERR_STOP(ncerr)
-          ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr)
        endif
-       ncerr = nf90_put_var(ncid, lon_varid, values=rad2dg*arrayr8 ); NC_ERR_STOP(ncerr)
+       ncerr = nf90_put_var(ncid, lon_varid, values=arrayr8 ); NC_ERR_STOP(ncerr)
     endif
 
     call ESMF_GridGetCoord(wrtGrid, coordDim=2, array=array, rc=rc); ESMF_ERR_RETURN(rc)
@@ -279,31 +291,19 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
     if (mype==0) then
        if (trim(output_grid) == 'gaussian_grid' .or. &
            trim(output_grid) == 'regional_latlon') then
-          ncerr = nf90_put_var(ncid, jm_varid, values=rad2dg*arrayr8(1,:)  ); NC_ERR_STOP(ncerr)
-          ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr)
-          ncerr = nf90_put_att(ncid, jm_varid, "long_name", "T-cell latiitude"); NC_ERR_STOP(ncerr)
-          ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees_N"); NC_ERR_STOP(ncerr)
-          ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr)
+          ncerr = nf90_put_var(ncid, jm_varid, values=arrayr8(1,:)  ); NC_ERR_STOP(ncerr)
        else if (trim(output_grid) == 'rotated_latlon') then
           do j=1,jm
              y(j) = lat1 + (lat2-lat1)/(jm-1) * (j-1)
           enddo
           ncerr = nf90_put_var(ncid, jm_varid, values=y  ); NC_ERR_STOP(ncerr)
-          ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr)
-          ncerr = nf90_put_att(ncid, jm_varid, "long_name", "rotated T-cell latiitude"); NC_ERR_STOP(ncerr)
-          ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees"); NC_ERR_STOP(ncerr)
-          ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr)
        else if (trim(output_grid) == 'lambert_conformal') then
           do j=1,jm
              y(j) = dy * (j-1)
           enddo
           ncerr = nf90_put_var(ncid, jm_varid, values=y  ); NC_ERR_STOP(ncerr)
-          ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr)
-          ncerr = nf90_put_att(ncid, jm_varid, "long_name", "y-coordinate of projection"); NC_ERR_STOP(ncerr)
-          ncerr = nf90_put_att(ncid, jm_varid, "units", "meters"); NC_ERR_STOP(ncerr)
-          ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr)
        endif
-       ncerr = nf90_put_var(ncid, lat_varid, values=rad2dg*arrayr8 ); NC_ERR_STOP(ncerr)
+       ncerr = nf90_put_var(ncid, lat_varid, values=arrayr8 ); NC_ERR_STOP(ncerr)
     endif
 
     do i=1, fieldCount
@@ -344,11 +344,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
                 arrayr4_3d = quantized(arrayr4_3d_save, nbits, dataMin, dataMax)
                 ! compute max abs compression error, save as a variable
                 ! attribute.
-                compress_err = maxval(abs(arrayr4_3d_save-arrayr4_3d))
-                ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr)
-                ncerr = nf90_put_att(ncid, varids(i), 'max_abs_compression_error', compress_err); NC_ERR_STOP(ncerr)
-                ncerr = nf90_put_att(ncid, varids(i), 'nbits', nbits); NC_ERR_STOP(ncerr)
-                ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr)
+                compress_err(i) = maxval(abs(arrayr4_3d_save-arrayr4_3d))
              endif
              ncerr = nf90_put_var(ncid, varids(i), values=arrayr4_3d, start=(/1,1,1/),count=(/im,jm,lm,1/) ); NC_ERR_STOP(ncerr)
            end if
@@ -363,6 +359,17 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
 
     end do
 
+    if (ideflate > 0 .and. nbits > 0 .and. mype == 0) then
+       ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr)
+       do i=1, fieldCount
+          if (compress_err(i) > 0) then
+             ncerr = nf90_put_att(ncid, varids(i), 'max_abs_compression_error', compress_err(i)); NC_ERR_STOP(ncerr)
+             ncerr = nf90_put_att(ncid, varids(i), 'nbits', nbits); NC_ERR_STOP(ncerr)
+          endif
+       enddo
+       ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr)
+    endif
+
     deallocate(arrayr4)
     deallocate(arrayr8)
     deallocate(arrayr4_3d,arrayr4_3d_save)
@@ -484,9 +491,9 @@ subroutine get_grid_attr(grid, prefix, ncid, varid, rc)
          else if (typekind==ESMF_TYPEKIND_R8) then
             call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", &
                                    name=trim(attName), value=varr8val, rc=rc); ESMF_ERR_RETURN(rc)
-            if (trim(attName) /= '_FillValue' .or. ideflate == 0) then
-              ! FIXME:  _FillValue must be cast to var type when using
-              ! NF90_NETCDF4. Until this is fixed, using netCDF default _FillValue.
+            if (trim(attName) /= '_FillValue') then
+              ! FIXME:  _FillValue must be cast to var type for recent versions
+              ! of netcdf
               ncerr = nf90_put_att(ncid, varid, trim(attName(ind+1:len(attName))), varr8val); NC_ERR_STOP(ncerr)
             endif
 
diff --git a/io/module_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90
index aad3991d0..84769eaea 100644
--- a/io/module_wrt_grid_comp.F90
+++ b/io/module_wrt_grid_comp.F90
@@ -183,7 +183,7 @@ subroutine wrt_initialize(wrt_comp, imp_state_write, exp_state_write, clock, rc)
       character(128) :: FBlist_outfilename(100), outfile_name
       character(128),dimension(:,:), allocatable    :: outfilename
       real(8), dimension(:),         allocatable    :: slat
-      real, dimension(:),            allocatable    :: lat, lon, axesdata
+      real(8), dimension(:),         allocatable    :: lat, lon
       real(ESMF_KIND_R8), dimension(:,:), pointer   :: lonPtr, latPtr
       real(ESMF_KIND_R8)                            :: rot_lon, rot_lat
       real(ESMF_KIND_R8)                            :: geo_lon, geo_lat
@@ -358,19 +358,20 @@ subroutine wrt_initialize(wrt_comp, imp_state_write, exp_state_write, clock, rc)
         wrt_int_state%latstart = lat(1)
         wrt_int_state%latlast  = lat(jmo)
         do j=1,imo
-          lon(j) = 360./real(imo) *real(j-1)
+          lon(j) = 360.d0/real(imo,8) *real(j-1,8)
         enddo
         wrt_int_state%lonstart = lon(1)
         wrt_int_state%lonlast  = lon(imo)
         do j=lbound(latPtr,2),ubound(latPtr,2)
           do i=lbound(lonPtr,1),ubound(lonPtr,1)
-            lonPtr(i,j) = 360./real(imo) * (i-1)
+            lonPtr(i,j) = 360.d0/real(imo,8) * real(i-1,8)
             latPtr(i,j) = lat(j)
           enddo
         enddo 
 !        print *,'aft wrtgrd, Gaussian, dimi,i=',lbound(lonPtr,1),ubound(lonPtr,1), &
 !         ' j=',lbound(lonPtr,2),ubound(lonPtr,2),'imo=',imo,'jmo=',jmo
-!        print *,'aft wrtgrd, lon=',lonPtr(lbound(lonPtr,1),lbound(lonPtr,2)), &
+!       if(wrt_int_state%mype==0) print *,'aft wrtgrd, lon=',lonPtr(1:5,1), &
+!        'lat=',latPtr(1,1:5),'imo,jmo=',imo,jmo
 !        lonPtr(lbound(lonPtr,1),ubound(lonPtr,2)),'lat=',latPtr(lbound(lonPtr,1),lbound(lonPtr,2)), &
 !        latPtr(lbound(lonPtr,1),ubound(lonPtr,2))
         wrt_int_state%lat_start = lbound(latPtr,2)
@@ -1622,13 +1623,14 @@ subroutine recover_fields(file_bundle,rc)
      character(100) fieldName,uwindname,vwindname
      type(ESMF_Field),   allocatable  :: fcstField(:)
      real(ESMF_KIND_R8), dimension(:,:),     pointer  :: lon, lat
+     real(ESMF_KIND_R8), dimension(:,:),     pointer  :: lonloc, latloc
      real(ESMF_KIND_R4), dimension(:,:),     pointer  :: pressfc
      real(ESMF_KIND_R4), dimension(:,:),     pointer  :: uwind2dr4,vwind2dr4
      real(ESMF_KIND_R4), dimension(:,:,:),   pointer  :: uwind3dr4,vwind3dr4
      real(ESMF_KIND_R4), dimension(:,:,:),   pointer  :: cart3dPtr2dr4
      real(ESMF_KIND_R4), dimension(:,:,:,:), pointer  :: cart3dPtr3dr4
      real(ESMF_KIND_R8), dimension(:,:,:,:), pointer  :: cart3dPtr3dr8
-     save lon, lat
+     save lonloc, latloc
      real(ESMF_KIND_R8) :: coslon, sinlon, sinlat
 !
 ! get filed count
@@ -1648,9 +1650,18 @@ subroutine recover_fields(file_bundle,rc)
 
        if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
 
-       lon = lon * pi/180.
-!     print *,'in 3DCartesian2wind, lon dim=',lbound(lon,1),ubound(lon,1),lbound(lon,2),ubound(lon,2), &
-!       'lon=',lon(lbound(lon,1),lbound(lon,2)), lon(ubound(lon,1),ubound(lon,2))
+       allocate(lonloc(lbound(lon,1):ubound(lon,1),lbound(lon,2):ubound(lon,2)))
+       istart = lbound(lon,1)
+       iend   = ubound(lon,1)
+       jstart = lbound(lon,2)
+       jend   = ubound(lon,2)
+!$omp parallel do default(none) shared(lon,lonloc,jstart,jend,istart,iend) &
+!$omp             private(i,j)
+       do j=jstart,jend
+        do i=istart,iend
+          lonloc(i,j) = lon(i,j) * pi/180.
+        enddo
+       enddo
 
        CALL ESMF_LogWrite("call recover field get coord 2",ESMF_LOGMSG_INFO,rc=RC)
 
@@ -1658,9 +1669,18 @@ subroutine recover_fields(file_bundle,rc)
 
        if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
 
-       lat = lat * pi/180.
-!     print *,'in 3DCartesian2wind, lat dim=',lbound(lat,1),ubound(lat,1),lbound(lat,2),ubound(lat,2), &
-!       'lat=',lat(lbound(lon,1),lbound(lon,2)), lat(ubound(lon,1),ubound(lon,2))
+       allocate(latloc(lbound(lat,1):ubound(lat,1),lbound(lat,2):ubound(lat,2)))
+       istart = lbound(lat,1)
+       iend   = ubound(lat,1)
+       jstart = lbound(lat,2)
+       jend   = ubound(lat,2)
+!$omp parallel do default(none) shared(lat,latloc,jstart,jend,istart,iend) &
+!$omp             private(i,j)
+       do j=jstart,jend
+        do i=istart,iend
+          latloc(i,j) = lat(i,j) * pi/180.d0
+        enddo
+       enddo
        first_getlatlon = .false.
      endif
 !
@@ -1718,18 +1738,18 @@ subroutine recover_fields(file_bundle,rc)
 ! update u , v wind
 !$omp parallel do default(shared) private(i,j,k,coslon,sinlon,sinlat)
              do k=kstart,kend
-!!$omp parallel do default(none) shared(uwind3dr4,vwind3dr4,lon,lat,cart3dPtr3dr4,jstart,jend,istart,iend,k) &
-!!$omp             private(i,j,coslon,sinlon,sinlat)
+!$omp parallel do default(none) shared(uwind3dr4,vwind3dr4,lonloc,latloc,cart3dPtr3dr4,jstart,jend,istart,iend,k) &
+!$omp             private(i,j,coslon,sinlon,sinlat)
                do j=jstart, jend
                  do i=istart, iend
-                  coslon = cos(lon(i,j))
-                  sinlon = sin(lon(i,j))
-                  sinlat = sin(lat(i,j))
+                  coslon = cos(lonloc(i,j))
+                  sinlon = sin(lonloc(i,j))
+                  sinlat = sin(latloc(i,j))
                   uwind3dr4(i,j,k) = cart3dPtr3dr4(1,i,j,k) * coslon           &
                                    + cart3dPtr3dr4(2,i,j,k) * sinlon
                   vwind3dr4(i,j,k) =-cart3dPtr3dr4(1,i,j,k) * sinlat*sinlon    &
                                    + cart3dPtr3dr4(2,i,j,k) * sinlat*coslon    &
-                                   + cart3dPtr3dr4(3,i,j,k) * cos(lat(i,j))
+                                   + cart3dPtr3dr4(3,i,j,k) * cos(latloc(i,j))
                  enddo
                enddo
              enddo
@@ -1749,18 +1769,18 @@ subroutine recover_fields(file_bundle,rc)
              call ESMF_FieldGet(ufield, localDe=0, farrayPtr=uwind2dr4,rc=rc)
              call ESMF_FieldGet(vfield, localDe=0, farrayPtr=vwind2dr4,rc=rc)
               ! update u , v wind
-!$omp parallel do default(none) shared(uwind2dr4,vwind2dr4,lon,lat,cart3dPtr2dr4,jstart,jend,istart,iend) &
+!$omp parallel do default(none) shared(uwind2dr4,vwind2dr4,lonloc,latloc,cart3dPtr2dr4,jstart,jend,istart,iend) &
 !$omp             private(i,j,k,coslon,sinlon,sinlat)
              do j=jstart, jend
                do i=istart, iend
-                  coslon = cos(lon(i,j))
-                  sinlon = sin(lon(i,j))
-                  sinlat = sin(lat(i,j))
+                  coslon = cos(lonloc(i,j))
+                  sinlon = sin(lonloc(i,j))
+                  sinlat = sin(latloc(i,j))
                   uwind2dr4(i,j) = cart3dPtr2dr4(1,i,j) * coslon         &
                                  + cart3dPtr2dr4(2,i,j) * sinlon
                   vwind2dr4(i,j) =-cart3dPtr2dr4(1,i,j) * sinlat*sinlon  &
                                  + cart3dPtr2dr4(2,i,j) * sinlat*coslon  &
-                                 + cart3dPtr2dr4(3,i,j) * cos(lat(i,j))
+                                 + cart3dPtr2dr4(3,i,j) * cos(latloc(i,j))
                enddo
              enddo
            endif
diff --git a/io/post_gfs.F90 b/io/post_gfs.F90
index 9f3a0a80f..daa6551bd 100644
--- a/io/post_gfs.F90
+++ b/io/post_gfs.F90
@@ -12,6 +12,7 @@ module post_gfs
   include 'mpif.h'
 
   integer mype, nbdl
+  logical setvar_atmfile, setvar_sfcfile, read_postcntrl
   public  post_run_gfs, post_getattr_gfs
 
   contains
@@ -28,9 +29,10 @@ subroutine post_run_gfs(wrt_int_state,mypei,mpicomp,lead_write,      &
 !
       use ctlblk_mod, only : komax,ifhr,ifmin,modelname,datapd,fld_info, &
                              npset,grib,gocart_on,icount_calmict, jsta,  &
-                             jend,im, nsoil
+                             jend,im, nsoil, filenameflat
       use gridspec_mod, only : maptype, gridtype
       use grib2_module, only : gribit2,num_pset,nrecout,first_grbtbl
+      use xml_perl_data,only : paramset
 !
 !-----------------------------------------------------------------------
 !
@@ -53,9 +55,8 @@ subroutine post_run_gfs(wrt_int_state,mypei,mpicomp,lead_write,      &
       integer n,nwtpg,ieof,lcntrl,ierr,i,j,k,jts,jte,mynsoil
       integer,allocatable  :: jstagrp(:),jendgrp(:)
       integer,save         :: kpo,kth,kpv
+      logical,save         :: log_postalct=.false.
       real,dimension(komax),save :: po, th, pv
-      logical,save   :: log_postalct=.false.
-      logical,save   :: setvar_atmfile=.false.,setvar_sfcfile=.false.
       logical        :: Log_runpost
       character(255) :: post_fname*255
 
@@ -124,6 +125,7 @@ subroutine post_run_gfs(wrt_int_state,mypei,mpicomp,lead_write,      &
 !
         log_postalct = .true.
         first_grbtbl = .true.
+        read_postcntrl = .true.
 !
       ENDIF
 !
@@ -135,6 +137,8 @@ subroutine post_run_gfs(wrt_int_state,mypei,mpicomp,lead_write,      &
       ifmin = mynfmin
       if (ifhr == 0 ) ifmin = 0
       if(mype==0) print *,'bf set_postvars,ifmin=',ifmin,'ifhr=',ifhr
+      setvar_atmfile=.false.
+      setvar_sfcfile=.false.
       call set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile,   &
            setvar_sfcfile)
 
@@ -145,8 +149,28 @@ subroutine post_run_gfs(wrt_int_state,mypei,mpicomp,lead_write,      &
 ! 20190807 no need to call microinit for GFDLMP
 !        call MICROINIT
 !
-        if(grib=="grib2" .and. first_grbtbl) then
-          call read_xml()
+        if(grib=="grib2" .and. read_postcntrl) then
+          if (ifhr == 0) then
+            filenameflat = 'postxconfig-NT_FH00.txt'
+            call read_xml()
+            if(mype==0) print *,'af read_xml at fh00,name=',trim(filenameflat)
+          else if(ifhr > 0) then
+            filenameflat = 'postxconfig-NT.txt'
+            if(size(paramset)>0) then
+              do i=1,size(paramset)
+                if (size(paramset(i)%param)>0) then
+                  deallocate(paramset(i)%param)
+                  nullify(paramset(i)%param)
+                endif
+              enddo
+              deallocate(paramset)
+              nullify(paramset)
+            endif
+            num_pset = 0
+            call read_xml()
+            if(mype==0) print *,'af read_xml,name=',trim(filenameflat),'ifhr=',ifhr
+            read_postcntrl = .false.
+          endif
         endif
 !
         IEOF  = 0
@@ -181,9 +205,6 @@ subroutine post_run_gfs(wrt_int_state,mypei,mpicomp,lead_write,      &
           endif
 !
         enddo
-!
-        setvar_atmfile = .false.
-        setvar_sfcfile = .false.
 !
       endif
 
@@ -335,7 +356,7 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile,   &
                              avgetrans, avgesnow, avgprec_cont, avgcprate_cont,&
                              avisbeamswin, avisdiffswin, airbeamswin, airdiffswin, &
                              alwoutc, alwtoac, aswoutc, aswtoac, alwinc, aswinc,& 
-                             avgpotevp, snoavg, si, cuppt
+                             avgpotevp, snoavg, ti, si, cuppt
       use soil,        only: sldpth, sh2o, smc, stc
       use masks,       only: lmv, lmh, htm, vtm, gdlat, gdlon, dx, dy, hbm2, sm, sice
       use ctlblk_mod,  only: im, jm, lm, lp1, jsta, jend, jsta_2l, jend_2u, jsta_m,jend_m, &
@@ -477,6 +498,7 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile,   &
           qs(i,j) = SPVAL
           twbs(i,j) = SPVAL
           qwbs(i,j) = SPVAL
+          ths(i,j) = SPVAL
         enddo
       enddo
 
@@ -1122,11 +1144,32 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile,   &
               !$omp parallel do private(i,j)
               do j=jsta,jend
                 do i=ista, iend
-                  sr(i,j) = arrayr42d(i,j)
+                  if (arrayr42d(i,j) /= spval) then
+                  !set range within (0,1)
+                    sr(i,j) = min(1.,max(0.,arrayr42d(i,j)))
+                  else
+                    sr(i,j) = spval
+                  endif
                 enddo
               enddo
             endif
 
+            ! sea ice skin temperature
+            if(trim(fieldname)=='tisfc') then
+              !$omp parallel do private(i,j)
+              do j=jsta,jend
+                do i=1,im
+                  if (arrayr42d(i,j) /= spval) then
+                    ti(i,j) = arrayr42d(i,j)
+                    if (sice(i,j) == spval .or. sice(i,j) == 0.) ti(i,j)=spval
+                  else
+                    ti(i,j) = spval
+                  endif
+                enddo
+              enddo
+!              print *,'in gfs_post, get tisfc=',maxval(ti), minval(ti)
+            endif
+
             ! vegetation fraction
             if(trim(fieldname)=='veg') then
               !$omp parallel do private(i,j)
@@ -1237,7 +1280,8 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile,   &
               do j=jsta,jend
                 do i=ista, iend
                   stc(i,j,1) = arrayr42d(i,j)
-                  if (sm(i,j) /= 0.0) stc(i,j,1) = spval
+                  !mask open water areas, combine with sea ice tmp
+                  if (sm(i,j) /= 0.0 .and. sice(i,j) ==0.) stc(i,j,1) = spval
                 enddo
               enddo
             endif
@@ -1248,7 +1292,8 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile,   &
               do j=jsta,jend
                 do i=ista, iend
                   stc(i,j,2) = arrayr42d(i,j)
-                  if (sm(i,j) /= 0.0) stc(i,j,2) = spval
+                  !mask open water areas, combine with sea ice tmp
+                  if (sm(i,j) /= 0.0 .and. sice(i,j) ==0.) stc(i,j,2) = spval
                 enddo
               enddo
             endif
@@ -1259,7 +1304,8 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile,   &
               do j=jsta,jend
                 do i=ista, iend
                   stc(i,j,3) = arrayr42d(i,j)
-                  if (sm(i,j) /= 0.0) stc(i,j,3) = spval
+                  !mask open water areas, combine with sea ice tmp
+                  if (sm(i,j) /= 0.0 .and. sice(i,j) ==0.) stc(i,j,3) = spval
                 enddo
               enddo
             endif
@@ -1270,7 +1316,8 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile,   &
               do j=jsta,jend
                 do i=ista, iend
                   stc(i,j,4) = arrayr42d(i,j)
-                  if (sm(i,j) /= 0.0) stc(i,j,4) = spval
+                  !mask open water areas, combine with sea ice tmp
+                  if (sm(i,j) /= 0.0 .and. sice(i,j) ==0.) stc(i,j,4) = spval
                 enddo
               enddo
             endif
@@ -2248,7 +2295,6 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile,   &
         do j=jsta,jend
           do i=1,im
             omga(i,j,l) = (-1.) * wh(i,j,l) * dpres(i,j,l)/zint(i,j,l)
-            pmid(i,j,l) = rgas*dpres(i,j,l) * t(i,j,l)*(q(i,j,l)*fv+1.0)/grav/zint(i,j,l)
             zint(i,j,l) = zint(i,j,l) + zint(i,j,l+1)
           enddo
         enddo
@@ -2271,6 +2317,15 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile,   &
         enddo
       end do
 
+!compute pmid from averaged two layer pint
+      do l=lm,1,-1
+        do j=jsta,jend
+          do i=1,im
+            pmid(i,j,l) = 0.5*(pint(i,j,l)+pint(i,j,l+1))
+          enddo
+        enddo
+      enddo
+
 !$omp parallel do private(i,j)
       do j=jsta,jend
         do i=1,im
@@ -2290,7 +2345,7 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile,   &
         end do
       end do
 
-! compute zmid  ??? where is definition of alpint(1)
+! compute zmid  
       do l=lm,1,-1
 !$omp parallel do private(i,j)
         do j=jsta,jend
@@ -2313,6 +2368,12 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile,   &
 !$omp parallel do private(i,j)
       do j=jsta,jend
         do i=ista, iend
+          !assign sst
+          if (sm(i,j) /= 0.0 .and. ths(i,j) /= spval) then
+             sst(i,j) = ths(i,j)
+          else
+             sst(i,j) = spval
+          endif
           if (ths(i,j) /= spval) then
             ths(i,j)  = ths(i,j)* (p1000/pint(i,j,lp1))**capa
             thz0(i,j) = ths(i,j)
diff --git a/module_fcst_grid_comp.F90 b/module_fcst_grid_comp.F90
index 6ff1f39c7..fef9698ab 100644
--- a/module_fcst_grid_comp.F90
+++ b/module_fcst_grid_comp.F90
@@ -22,11 +22,12 @@ module module_fcst_grid_comp
 !
   use time_manager_mod,   only: time_type, set_calendar_type, set_time,    &
                                 set_date, days_in_month, month_name,       &
-                                operator(+), operator (<), operator (>),   &
-                                operator (/=), operator (/), operator (==),&
-                                operator (*), THIRTY_DAY_MONTHS, JULIAN,   &
-                                NOLEAP, NO_CALENDAR, date_to_string,       &
-                                get_date
+                                operator(+), operator(-), operator (<),    &
+                                operator (>), operator (/=), operator (/), &
+                                operator (==), operator (*),               &
+                                THIRTY_DAY_MONTHS, JULIAN, NOLEAP,         &
+                                NO_CALENDAR, date_to_string, get_date,     &
+                                get_time
 
   use  atmos_model_mod,   only: atmos_model_init, atmos_model_end,         &
                                 get_atmos_model_ungridded_dim,             &
@@ -70,7 +71,8 @@ module module_fcst_grid_comp
                                 iau_offset
   use module_fv3_config, only:  dt_atmos, calendar, restart_interval,             &
                                 quilting, calendar_type, cpl,                     &
-                                cplprint_flag, force_date_from_configure
+                                cplprint_flag, force_date_from_configure,         &
+                                num_restart_interval, frestart, restart_endfcst
 !
 !-----------------------------------------------------------------------
 !
@@ -88,7 +90,8 @@ module module_fcst_grid_comp
     type(atmos_data_type)  :: Atm
     type(time_type)        :: Time_atmos, Time_init, Time_end,  &
                               Time_step_atmos, Time_step_ocean, &
-                              Time_restart, Time_step_restart
+                              Time_restart, Time_step_restart,  &
+                              Time_atstart
     integer :: num_atmos_calls, ret, intrm_rst
   end type
 
@@ -179,12 +182,11 @@ subroutine fcst_initialize(fcst_comp, importState, exportState, clock, rc)
 
     integer                                :: Run_length
     integer,dimension(6)                   :: date, date_end
-    integer                                :: res_intvl
     integer                                :: mpi_comm_comp
 !
     logical,save                           :: first=.true.
     character(len=9) :: month
-    integer :: initClock, unit, nfhour
+    integer :: initClock, unit, nfhour, total_inttime
     integer :: mype, ntasks
     character(3) cfhour
     character(4) dateSY
@@ -203,7 +205,8 @@ subroutine fcst_initialize(fcst_comp, importState, exportState, clock, rc)
     real(ESMF_KIND_R8),parameter :: dtor = 180.0_ESMF_KIND_R8 / 3.1415926535897931_ESMF_KIND_R8
     integer :: jsc, jec, isc, iec, nlev
     type(domain2D) :: domain
-    integer :: n, fcstNpes
+    integer :: n, fcstNpes, tmpvar
+    logical :: single_restart
     integer, allocatable, dimension(:) :: isl, iel, jsl, jel
     integer, allocatable, dimension(:,:,:) :: deBlockList
 
@@ -271,8 +274,9 @@ subroutine fcst_initialize(fcst_comp, importState, exportState, clock, rc)
 !----------------------------------------------------------------------- 
 !
     call ESMF_ClockGet(clock, CurrTime=CurrTime, StartTime=StartTime, &
-                       StopTime=StopTime, RunDuration=RunDuration, rc=rc)
+                       StopTime=StopTime, rc=rc)
     if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+    RunDuration = StopTime - CurrTime
 
     date_init = 0
     call ESMF_TimeGet (StartTime,                      &
@@ -317,16 +321,46 @@ subroutine fcst_initialize(fcst_comp, importState, exportState, clock, rc)
 !
     atm_int_state%Time_step_atmos = set_time (dt_atmos,0)
     atm_int_state%num_atmos_calls = Run_length / dt_atmos
+    atm_int_state%Time_atstart = atm_int_state%Time_atmos
     if (mype == 0) write(0,*)'num_atmos_calls=',atm_int_state%num_atmos_calls,'time_init=', &
                     date_init,'time_atmos=',date,'time_end=',date_end,'dt_atmos=',dt_atmos, &
                     'Run_length=',Run_length
-     res_intvl = restart_interval*3600
-     atm_int_state%Time_step_restart = set_time (res_intvl, 0)
-     atm_int_state%Time_restart      = atm_int_state%Time_atmos + atm_int_state%Time_step_restart
-     atm_int_state%intrm_rst         = 0
-     if (res_intvl>0) atm_int_state%intrm_rst = 1
-     atm_int_state%Atm%iau_offset    = iau_offset
-!
+    frestart = 0
+    single_restart = .false.
+    call get_time(atm_int_state%Time_end - atm_int_state%Time_atstart,total_inttime)
+    if(num_restart_interval == 2) then
+      if(restart_interval(2)== -1) single_restart = .true.
+    endif
+    if(single_restart) then
+      frestart(1) =  restart_interval(1) * 3600
+    elseif ( num_restart_interval == 1) then
+      if(restart_interval(1) == 0) then
+        frestart(1) = total_inttime
+      else if(restart_interval(1) > 0) then
+        tmpvar = restart_interval(1) * 3600
+        frestart(1) = tmpvar
+        atm_int_state%Time_step_restart = set_time (tmpvar, 0)
+        atm_int_state%Time_restart      = atm_int_state%Time_atstart + atm_int_state%Time_step_restart
+        i = 2
+        do while ( atm_int_state%Time_restart < atm_int_state%Time_end ) 
+          frestart(i) = frestart(i-1) + tmpvar
+          atm_int_state%Time_restart = atm_int_state%Time_restart + atm_int_state%Time_step_restart
+           i = i + 1
+        enddo
+      endif
+    else if(num_restart_interval > 1) then
+      do i=1,num_restart_interval
+        frestart(i) = restart_interval(i) * 3600
+      enddo
+    endif
+    restart_endfcst = .false.
+    if ( ANY(frestart(:) == total_inttime) ) restart_endfcst = .true. 
+    if (mype == 0) print *,'frestart=',frestart(1:10)/3600, 'restart_endfcst=',restart_endfcst, &
+      'total_inttime=',total_inttime
+       
+    atm_int_state%intrm_rst         = 0
+    if (frestart(1)>0) atm_int_state%intrm_rst = 1
+    atm_int_state%Atm%iau_offset    = iau_offset
 !
 !----- write time stamps (for start time and end time) ------
 
@@ -737,9 +771,10 @@ subroutine fcst_run_phase_2(fcst_comp, importState, exportState,clock,rc)
 !-----------------------------------------------------------------------
 !***  local variables
 !
-      integer                    :: i,j, mype, na, date(6)
+      integer                    :: i,j, mype, na, date(6), seconds
       character(20)              :: compname
-
+ 
+      type(time_type)            :: restart_inctime
       type(ESMF_Time)            :: currtime
       integer(kind=ESMF_KIND_I8) :: ntimestep_esmf
       character(len=64)          :: timestamp
@@ -776,13 +811,16 @@ subroutine fcst_run_phase_2(fcst_comp, importState, exportState,clock,rc)
 
 !--- intermediate restart
       if (atm_int_state%intrm_rst>0) then
-        if ((na /= atm_int_state%num_atmos_calls) .and.   &
-           (atm_int_state%Time_atmos == atm_int_state%Time_restart)) then
-          timestamp = date_to_string (atm_int_state%Time_restart)
-          call atmos_model_restart(atm_int_state%Atm, timestamp)
-
-          call wrt_atmres_timestamp(atm_int_state,timestamp)
-          atm_int_state%Time_restart = atm_int_state%Time_restart + atm_int_state%Time_step_restart
+        if (na /= atm_int_state%num_atmos_calls-1) then 
+          call get_time(atm_int_state%Time_atmos - atm_int_state%Time_atstart, seconds)
+          if (ANY(frestart(:) == seconds)) then
+            restart_inctime = set_time(seconds, 0)
+            atm_int_state%Time_restart = atm_int_state%Time_atstart + restart_inctime
+            timestamp = date_to_string (atm_int_state%Time_restart)
+            call atmos_model_restart(atm_int_state%Atm, timestamp)
+
+            call wrt_atmres_timestamp(atm_int_state,timestamp)
+          endif
         endif
       endif
 !
@@ -847,20 +885,21 @@ subroutine fcst_finalize(fcst_comp, importState, exportState,clock,rc)
                          'final time does not match expected ending time', WARNING)
 
 !*** write restart file
-
-      call get_date (atm_int_state%Time_atmos, date(1), date(2), date(3),  &
+      if( restart_endfcst ) then
+        call get_date (atm_int_state%Time_atmos, date(1), date(2), date(3),  &
                                date(4), date(5), date(6))
-      call mpp_open( unit, 'RESTART/coupler.res', nohdrs=.TRUE. )
-      if (mpp_pe() == mpp_root_pe())then
-         write( unit, '(i6,8x,a)' )calendar_type, &
+        call mpp_open( unit, 'RESTART/coupler.res', nohdrs=.TRUE. )
+        if (mpp_pe() == mpp_root_pe())then
+          write( unit, '(i6,8x,a)' )calendar_type, &
               '(Calendar: no_calendar=0, thirty_day_months=1, julian=2, gregorian=3, noleap=4)'
 
-         write( unit, '(6i6,8x,a)' )date_init, &
+          write( unit, '(6i6,8x,a)' )date_init, &
               'Model start time:   year, month, day, hour, minute, second'
-         write( unit, '(6i6,8x,a)' )date, &
+          write( unit, '(6i6,8x,a)' )date, &
               'Current model time: year, month, day, hour, minute, second'
+        endif
+        call mpp_close(unit)
       endif
-      call mpp_close(unit)
 !
       call diag_manager_end(atm_int_state%Time_atmos )
 
diff --git a/module_fv3_config.F90 b/module_fv3_config.F90
index eb1bb2036..2e0fd8883 100644
--- a/module_fv3_config.F90
+++ b/module_fv3_config.F90
@@ -13,11 +13,10 @@ module module_fv3_config
 
   implicit none
 !
-  
-  integer                  :: restart_interval
-!
   integer                  :: nfhout, nfhout_hf, nsout, dt_atmos
   integer                  :: nfhmax, nfhmax_hf, first_kdt
+  integer                  :: num_restart_interval
+  integer                  :: frestart(1000)
   type(ESMF_Alarm)         :: alarm_output_hf, alarm_output
   type(ESMF_TimeInterval)  :: output_hfmax
   type(ESMF_TimeInterval)  :: output_interval,output_interval_hf
@@ -25,7 +24,9 @@ module module_fv3_config
   logical                  :: cpl, cplprint_flag
   logical                  :: quilting, output_1st_tstep_rst
   logical                  :: force_date_from_configure
+  logical                  :: restart_endfcst
 !
+  integer,dimension(:),allocatable                :: restart_interval
   character(esmf_maxstr),dimension(:),allocatable :: filename_base
   character(17)            :: calendar='                 '
   integer                  :: calendar_type = -99