From c92e0eee7d297f6050465672dcd6f9f88bd1ca58 Mon Sep 17 00:00:00 2001
From: jswhit2 <Jeffrey.S.Whitaker@noaa.gov>
Date: Thu, 6 Feb 2020 20:58:06 +0000
Subject: [PATCH 01/10] add dry mass fixer for IAU (requires change to
 GFS_typedefs.F90)

---
 driver/fvGFS/atmosphere.F90 | 75 +++++++++++++++++++++++++++++++++++--
 1 file changed, 72 insertions(+), 3 deletions(-)

diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90
index 0d81b6cb1..49b2ac072 100644
--- a/driver/fvGFS/atmosphere.F90
+++ b/driver/fvGFS/atmosphere.F90
@@ -179,11 +179,11 @@ module atmosphere_mod
 use fv_fill_mod,        only: fill_gfs
 use fv_dynamics_mod,    only: fv_dynamics
 use fv_nesting_mod,     only: twoway_nesting
-use fv_diagnostics_mod, only: fv_diag_init, fv_diag, fv_time, prt_maxmin, prt_height
+use fv_diagnostics_mod, only: fv_diag_init, fv_diag, fv_time, prt_maxmin, prt_height, prt_mass
 use fv_nggps_diags_mod, only: fv_nggps_diag_init, fv_nggps_diag, fv_nggps_tavg
 use fv_restart_mod,     only: fv_restart, fv_write_restart
 use fv_timing_mod,      only: timing_on, timing_off
-use fv_mp_mod,          only: switch_current_Atm
+use fv_mp_mod,          only: switch_current_Atm, is_master
 use fv_sg_mod,          only: fv_subgrid_z
 use fv_update_phys_mod, only: fv_update_phys
 use fv_nwp_nudge_mod,   only: fv_nwp_nudge_init, fv_nwp_nudge_end, do_adiabatic_init
@@ -194,6 +194,7 @@ module atmosphere_mod
                               a_step, p_step, current_time_in_seconds
 
 use mpp_domains_mod,    only:  mpp_get_data_domain, mpp_get_compute_domain
+use fv_grid_utils_mod,  only: g_sum
 
 implicit none
 private
@@ -1395,6 +1396,8 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
    integer :: i, j, ix, k, k1, n, w_diff, nt_dyn, iq
    integer :: nb, blen, nwat, dnats, nq_adv
    real(kind=kind_phys):: rcp, q0, qwat(nq), qt, rdt
+   real psum, qsum, psumb, qsumb, psuma, qsuma, betad
+   real, allocatable, dimension(:,:,:,:) :: qtmp
    Time_prev = Time
    Time_next = Time + Time_step_atmos
    rdt = 1.d0 / dt_atmos
@@ -1407,6 +1410,32 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
    if( nq<3 ) call mpp_error(FATAL, 'GFS phys must have 3 interactive tracers')
 
    if (IAU_Data%in_interval) then
+
+      if (IAU_Data%drymassfixer) then
+         allocate ( qtmp(isd:ied  ,jsd:jed  ,npz, nwat) )
+         ! global mean total pressure and water before IAU
+         psumb = g_sum(Atm(n)%domain,sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz),dim=3),&
+                 isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) 
+         do nb = 1,Atm_block%nblks
+            blen = Atm_block%blksz(nb)
+            do k=1,npz
+               if(flip_vc) then
+                 k1 = npz+1-k !reverse the k direction 
+               else
+                 k1 = k
+               endif
+               do ix = 1, blen
+                   i = Atm_block%index(nb)%ii(ix)
+                   j = Atm_block%index(nb)%jj(ix)
+                   qtmp(i,j,k1,1:nwat) = IPD_Data(nb)%Stateout%gq0(ix,k,1:nwat)
+               enddo
+            enddo
+         enddo
+         qsumb = g_sum(Atm(n)%domain,&
+                 sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(qtmp(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
+                 isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) 
+      endif
+
 !     IAU increments are in units of 1/sec
 
 !     add analysis increment to u,v,t tendencies
@@ -1447,6 +1476,35 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
             enddo
          enddo
       enddo
+
+      if (IAU_data%drymassfixer) then
+         ! global mean total pressure and water after IAU
+         psuma = g_sum(Atm(n)%domain,sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz),dim=3),&
+                 isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) 
+         do nb = 1,Atm_block%nblks
+            blen = Atm_block%blksz(nb)
+            do k=1,npz
+               if(flip_vc) then
+                 k1 = npz+1-k !reverse the k direction 
+               else
+                 k1 = k
+               endif
+               do ix = 1, blen
+                   i = Atm_block%index(nb)%ii(ix)
+                   j = Atm_block%index(nb)%jj(ix)
+                   qtmp(i,j,k1,1:nwat) = IPD_Data(nb)%Stateout%gq0(ix,k,1:nwat)
+               enddo
+            enddo
+         enddo
+         qsuma = g_sum(Atm(n)%domain,&
+                 sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(qtmp(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
+                 isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) 
+         ! rescale water consitituents to ensure IAU doesn't change dry mass
+         betad = (psuma - (psumb - qsumb))/qsuma
+         IPD_Data(nb)%Stateout%gq0 = betad*IPD_Data(nb)%Stateout%gq0
+         deallocate(qtmp)
+      endif
+
    endif
 
    call set_domain ( Atm(mytile)%domain )
@@ -1501,7 +1559,7 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
 #ifdef MULTI_GASES
          q0 = Atm(n)%delp(i,j,k1)*(1.-sum(Atm(n)%q(i,j,k1,1:max(nwat,num_gas)))) + sum(qwat(1:max(nwat,num_gas)))
 #else
-         qt = sum(qwat(1:nwat))
+         qt = sum(qwat(1:nwat)) 
          q0 = Atm(n)%delp(i,j,k1)*(1.-sum(Atm(n)%q(i,j,k1,1:nwat))) + qt 
 #endif
          Atm(n)%delp(i,j,k1) = q0
@@ -1575,6 +1633,17 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
        call timing_off('FV_UPDATE_PHYS')
    call mpp_clock_end (id_dynam)
 
+   ! global mean total pressure
+   !psum = g_sum(Atm(n)%domain,sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz),dim=3),&
+   !       isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) + Atm(n)%ptop
+   !! global mean total water
+   !qsum = g_sum(Atm(n)%domain,&
+   !       sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
+   !       isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) 
+   !if( is_master() ) then
+   !     print *,'dry surface pressure after IAU/physics = ',psum-qsum
+   !endif
+
 !--- nesting update after updating atmospheric variables with
 !--- physics tendencies
     if (ngrids > 1 .and. p_split > 0) then

From 6f5f192d9822a9aed5bb745cc23e95966a296780 Mon Sep 17 00:00:00 2001
From: jswhit2 <Jeffrey.S.Whitaker@noaa.gov>
Date: Thu, 6 Feb 2020 23:58:14 +0000
Subject: [PATCH 02/10] update - still not quite working

---
 driver/fvGFS/atmosphere.F90 | 68 ++++++++++++++++++++++++-------------
 1 file changed, 45 insertions(+), 23 deletions(-)

diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90
index 49b2ac072..115c4e02b 100644
--- a/driver/fvGFS/atmosphere.F90
+++ b/driver/fvGFS/atmosphere.F90
@@ -1398,6 +1398,7 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
    real(kind=kind_phys):: rcp, q0, qwat(nq), qt, rdt
    real psum, qsum, psumb, qsumb, psuma, qsuma, betad
    real, allocatable, dimension(:,:,:,:) :: qtmp
+   logical in_iau_interval,use_iau_massfixer
    Time_prev = Time
    Time_next = Time + Time_step_atmos
    rdt = 1.d0 / dt_atmos
@@ -1406,9 +1407,18 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
    nwat = Atm(n)%flagstruct%nwat
    dnats = Atm(mytile)%flagstruct%dnats
    nq_adv = nq - dnats
+   in_iau_interval=IAU_Data%in_interval
+   use_iau_massfixer=IAU_Data%drymassfixer
 
    if( nq<3 ) call mpp_error(FATAL, 'GFS phys must have 3 interactive tracers')
 
+!SJL: perform vertical filling to fix the negative humidity if the SAS convection scheme is used
+!     This call may be commented out if RAS or other positivity-preserving CPS is used.
+!  do nb = 1,Atm_block%nblks
+!    blen = Atm_block%blksz(nb)
+!    call fill_gfs(blen, npz, IPD_Data(nb)%Statein%prsi, IPD_Data(nb)%Stateout%gq0, 1.e-9_kind_phys)
+!  enddo
+
    if (IAU_Data%in_interval) then
 
       if (IAU_Data%drymassfixer) then
@@ -1500,8 +1510,13 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
                  sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(qtmp(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
                  isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) 
          ! rescale water consitituents to ensure IAU doesn't change dry mass
+         ! following
+         ! https://onlinelibrary.wiley.com/doi/full/10.1111/j.1600-0870.2007.00299.x
          betad = (psuma - (psumb - qsumb))/qsuma
-         IPD_Data(nb)%Stateout%gq0 = betad*IPD_Data(nb)%Stateout%gq0
+         if (is_master()) print *,'betad = ',betad
+         do nb = 1,Atm_block%nblks
+            IPD_Data(nb)%Stateout%gq0(:,:,1:nwat) = betad*IPD_Data(nb)%Stateout%gq0(:,:,1:nwat)
+         enddo
          deallocate(qtmp)
       endif
 
@@ -1517,7 +1532,7 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
 #ifdef MULTI_GASES
 !$OMP                      num_gas,                                                      &
 #endif
-!$OMP                      snowwat, graupel, nq_adv, flip_vc)   &
+!$OMP                      in_iau_interval, use_iau_massfixer,snowwat, graupel, nq_adv, flip_vc)   &
 !$OMP             private (nb, blen, i, j, k, k1, ix, q0, qwat, qt)
    do nb = 1,Atm_block%nblks
 
@@ -1527,11 +1542,11 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
      call fill_gfs(blen, npz, IPD_Data(nb)%Statein%prsi, IPD_Data(nb)%Stateout%gq0, 1.e-9_kind_phys)
 
      do k = 1, npz
-           if(flip_vc) then
-             k1 = npz+1-k !reverse the k direction 
-           else
-             k1 = k
-           endif
+       if(flip_vc) then
+         k1 = npz+1-k !reverse the k direction 
+       else
+         k1 = k
+       endif
        do ix = 1, blen
          i = Atm_block%index(nb)%ii(ix)
          j = Atm_block%index(nb)%jj(ix)
@@ -1551,6 +1566,9 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
            q0 = IPD_Data(nb)%Statein%prsi(ix,k+1) - IPD_Data(nb)%Statein%prsi(ix,k)
          endif
          qwat(1:nq_adv) = q0*IPD_Data(nb)%Stateout%gq0(ix,k,1:nq_adv)
+         if (in_iau_interval .and. use_iau_massfixer) then
+         Atm(n)%q(i,j,k1,1:nq_adv) = IPD_Data(nb)%Stateout%gq0(ix,k,1:nq_adv)
+         else
 ! **********************************************************************************************************
 ! Dry mass: the following way of updating delp is key to mass conservation with hybrid 32-64 bit computation
 ! **********************************************************************************************************
@@ -1565,6 +1583,7 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
          Atm(n)%delp(i,j,k1) = q0
          Atm(n)%q(i,j,k1,1:nq_adv) = qwat(1:nq_adv) / q0
 !        if (dnats .gt. 0) Atm(n)%q(i,j,k1,nq_adv+1:nq) = IPD_Data(nb)%Stateout%gq0(ix,k,nq_adv+1:nq)
+         endif
        enddo
      enddo
 
@@ -1573,11 +1592,11 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
      !--- See Note in statein...
      do iq = nq+1, ncnst
        do k = 1, npz
-           if(flip_vc) then
-             k1 = npz+1-k !reverse the k direction 
-           else
-             k1 = k
-           endif
+         if(flip_vc) then
+           k1 = npz+1-k !reverse the k direction 
+         else
+           k1 = k
+         endif
          do ix = 1, blen
            i = Atm_block%index(nb)%ii(ix)
            j = Atm_block%index(nb)%jj(ix)
@@ -1588,6 +1607,20 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
 
    enddo  ! nb-loop
 
+! print out dry mass inside IAU interval.
+   if (IAU_Data%in_interval) then
+      ! global mean total pressure
+      psum = g_sum(Atm(n)%domain,sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz),dim=3),&
+             isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) + Atm(n)%ptop
+      ! global mean total water
+      qsum = g_sum(Atm(n)%domain,&
+             sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
+             isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) 
+      if( is_master() ) then
+           print *,'dry surface pressure in IAU interval = ',psum-qsum
+      endif
+   endif
+
    call timing_off('GFS_TENDENCIES')
 
    w_diff = get_tracer_index (MODEL_ATMOS, 'w_diff' )
@@ -1633,17 +1666,6 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
        call timing_off('FV_UPDATE_PHYS')
    call mpp_clock_end (id_dynam)
 
-   ! global mean total pressure
-   !psum = g_sum(Atm(n)%domain,sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz),dim=3),&
-   !       isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) + Atm(n)%ptop
-   !! global mean total water
-   !qsum = g_sum(Atm(n)%domain,&
-   !       sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
-   !       isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) 
-   !if( is_master() ) then
-   !     print *,'dry surface pressure after IAU/physics = ',psum-qsum
-   !endif
-
 !--- nesting update after updating atmospheric variables with
 !--- physics tendencies
     if (ngrids > 1 .and. p_split > 0) then

From fd00e4abfc7658649525bbc957bc9ef1f0f6debe Mon Sep 17 00:00:00 2001
From: jswhit2 <Jeffrey.S.Whitaker@noaa.gov>
Date: Fri, 7 Feb 2020 01:47:38 +0000
Subject: [PATCH 03/10] update 2 - now working

---
 driver/fvGFS/atmosphere.F90 | 95 ++++++++++---------------------------
 1 file changed, 24 insertions(+), 71 deletions(-)

diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90
index 115c4e02b..48a023125 100644
--- a/driver/fvGFS/atmosphere.F90
+++ b/driver/fvGFS/atmosphere.F90
@@ -1396,9 +1396,7 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
    integer :: i, j, ix, k, k1, n, w_diff, nt_dyn, iq
    integer :: nb, blen, nwat, dnats, nq_adv
    real(kind=kind_phys):: rcp, q0, qwat(nq), qt, rdt
-   real psum, qsum, psumb, qsumb, psuma, qsuma, betad
-   real, allocatable, dimension(:,:,:,:) :: qtmp
-   logical in_iau_interval,use_iau_massfixer
+   real psum, qsum, psumb, qsumb, betad
    Time_prev = Time
    Time_next = Time + Time_step_atmos
    rdt = 1.d0 / dt_atmos
@@ -1407,43 +1405,22 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
    nwat = Atm(n)%flagstruct%nwat
    dnats = Atm(mytile)%flagstruct%dnats
    nq_adv = nq - dnats
-   in_iau_interval=IAU_Data%in_interval
-   use_iau_massfixer=IAU_Data%drymassfixer
 
    if( nq<3 ) call mpp_error(FATAL, 'GFS phys must have 3 interactive tracers')
 
-!SJL: perform vertical filling to fix the negative humidity if the SAS convection scheme is used
-!     This call may be commented out if RAS or other positivity-preserving CPS is used.
-!  do nb = 1,Atm_block%nblks
-!    blen = Atm_block%blksz(nb)
-!    call fill_gfs(blen, npz, IPD_Data(nb)%Statein%prsi, IPD_Data(nb)%Stateout%gq0, 1.e-9_kind_phys)
-!  enddo
 
    if (IAU_Data%in_interval) then
 
       if (IAU_Data%drymassfixer) then
-         allocate ( qtmp(isd:ied  ,jsd:jed  ,npz, nwat) )
          ! global mean total pressure and water before IAU
          psumb = g_sum(Atm(n)%domain,sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz),dim=3),&
                  isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) 
-         do nb = 1,Atm_block%nblks
-            blen = Atm_block%blksz(nb)
-            do k=1,npz
-               if(flip_vc) then
-                 k1 = npz+1-k !reverse the k direction 
-               else
-                 k1 = k
-               endif
-               do ix = 1, blen
-                   i = Atm_block%index(nb)%ii(ix)
-                   j = Atm_block%index(nb)%jj(ix)
-                   qtmp(i,j,k1,1:nwat) = IPD_Data(nb)%Stateout%gq0(ix,k,1:nwat)
-               enddo
-            enddo
-         enddo
          qsumb = g_sum(Atm(n)%domain,&
-                 sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(qtmp(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
+                 sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
                  isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) 
+         if( is_master() ) then
+            print *,'dry surface pressure in IAU interval (before) = ',psumb++Atm%ptop-qsumb
+         endif
       endif
 
 !     IAU increments are in units of 1/sec
@@ -1487,39 +1464,6 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
          enddo
       enddo
 
-      if (IAU_data%drymassfixer) then
-         ! global mean total pressure and water after IAU
-         psuma = g_sum(Atm(n)%domain,sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz),dim=3),&
-                 isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) 
-         do nb = 1,Atm_block%nblks
-            blen = Atm_block%blksz(nb)
-            do k=1,npz
-               if(flip_vc) then
-                 k1 = npz+1-k !reverse the k direction 
-               else
-                 k1 = k
-               endif
-               do ix = 1, blen
-                   i = Atm_block%index(nb)%ii(ix)
-                   j = Atm_block%index(nb)%jj(ix)
-                   qtmp(i,j,k1,1:nwat) = IPD_Data(nb)%Stateout%gq0(ix,k,1:nwat)
-               enddo
-            enddo
-         enddo
-         qsuma = g_sum(Atm(n)%domain,&
-                 sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(qtmp(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
-                 isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) 
-         ! rescale water consitituents to ensure IAU doesn't change dry mass
-         ! following
-         ! https://onlinelibrary.wiley.com/doi/full/10.1111/j.1600-0870.2007.00299.x
-         betad = (psuma - (psumb - qsumb))/qsuma
-         if (is_master()) print *,'betad = ',betad
-         do nb = 1,Atm_block%nblks
-            IPD_Data(nb)%Stateout%gq0(:,:,1:nwat) = betad*IPD_Data(nb)%Stateout%gq0(:,:,1:nwat)
-         enddo
-         deallocate(qtmp)
-      endif
-
    endif
 
    call set_domain ( Atm(mytile)%domain )
@@ -1532,7 +1476,7 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
 #ifdef MULTI_GASES
 !$OMP                      num_gas,                                                      &
 #endif
-!$OMP                      in_iau_interval, use_iau_massfixer,snowwat, graupel, nq_adv, flip_vc)   &
+!$OMP                      snowwat, graupel, nq_adv, flip_vc)   &
 !$OMP             private (nb, blen, i, j, k, k1, ix, q0, qwat, qt)
    do nb = 1,Atm_block%nblks
 
@@ -1566,9 +1510,6 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
            q0 = IPD_Data(nb)%Statein%prsi(ix,k+1) - IPD_Data(nb)%Statein%prsi(ix,k)
          endif
          qwat(1:nq_adv) = q0*IPD_Data(nb)%Stateout%gq0(ix,k,1:nq_adv)
-         if (in_iau_interval .and. use_iau_massfixer) then
-         Atm(n)%q(i,j,k1,1:nq_adv) = IPD_Data(nb)%Stateout%gq0(ix,k,1:nq_adv)
-         else
 ! **********************************************************************************************************
 ! Dry mass: the following way of updating delp is key to mass conservation with hybrid 32-64 bit computation
 ! **********************************************************************************************************
@@ -1583,7 +1524,6 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
          Atm(n)%delp(i,j,k1) = q0
          Atm(n)%q(i,j,k1,1:nq_adv) = qwat(1:nq_adv) / q0
 !        if (dnats .gt. 0) Atm(n)%q(i,j,k1,nq_adv+1:nq) = IPD_Data(nb)%Stateout%gq0(ix,k,nq_adv+1:nq)
-         endif
        enddo
      enddo
 
@@ -1607,18 +1547,31 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
 
    enddo  ! nb-loop
 
-! print out dry mass inside IAU interval.
-   if (IAU_Data%in_interval) then
+! dry mass fixer in IAU interval following
+! https://onlinelibrary.wiley.com/doi/full/10.1111/j.1600-0870.2007.00299.x
+   if (IAU_Data%in_interval .and. IAU_data%drymassfixer) then
       ! global mean total pressure
       psum = g_sum(Atm(n)%domain,sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz),dim=3),&
-             isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) + Atm(n)%ptop
-      ! global mean total water
+             isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) 
+      ! global mean total water (before adjustment)
       qsum = g_sum(Atm(n)%domain,&
              sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
              isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) 
       if( is_master() ) then
-           print *,'dry surface pressure in IAU interval = ',psum-qsum
+         print *,'dry surface pressure in IAU interval (after) =',&
+         psum+Atm%ptop-qsum
       endif
+      betad = (psum - (psumb - qsumb))/qsum
+      !if (is_master()) print *,'betad = ',betad
+      Atm(n)%q(:,:,:,1:nwat) = betad*Atm(n)%q(:,:,:,1:nwat)
+      ! global mean total water after adjustment
+      !qsum = g_sum(Atm(n)%domain,&
+      !       sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
+      !       isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) 
+      !if( is_master() ) then
+      !   print *,'dry surface pressure in IAU interval (after 2) =',&
+      !   psum+Atm%ptop-qsum
+      !endif
    endif
 
    call timing_off('GFS_TENDENCIES')

From a66f3a28c9cbaf8804247889bd9cd10c384ca023 Mon Sep 17 00:00:00 2001
From: jswhit2 <Jeffrey.S.Whitaker@noaa.gov>
Date: Fri, 7 Feb 2020 05:02:32 +0000
Subject: [PATCH 04/10] remove whitespace changes

---
 driver/fvGFS/atmosphere.F90 | 7 ++-----
 1 file changed, 2 insertions(+), 5 deletions(-)

diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90
index 48a023125..74c342d47 100644
--- a/driver/fvGFS/atmosphere.F90
+++ b/driver/fvGFS/atmosphere.F90
@@ -179,7 +179,7 @@ module atmosphere_mod
 use fv_fill_mod,        only: fill_gfs
 use fv_dynamics_mod,    only: fv_dynamics
 use fv_nesting_mod,     only: twoway_nesting
-use fv_diagnostics_mod, only: fv_diag_init, fv_diag, fv_time, prt_maxmin, prt_height, prt_mass
+use fv_diagnostics_mod, only: fv_diag_init, fv_diag, fv_time, prt_maxmin, prt_height
 use fv_nggps_diags_mod, only: fv_nggps_diag_init, fv_nggps_diag, fv_nggps_tavg
 use fv_restart_mod,     only: fv_restart, fv_write_restart
 use fv_timing_mod,      only: timing_on, timing_off
@@ -1408,9 +1408,7 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
 
    if( nq<3 ) call mpp_error(FATAL, 'GFS phys must have 3 interactive tracers')
 
-
    if (IAU_Data%in_interval) then
-
       if (IAU_Data%drymassfixer) then
          ! global mean total pressure and water before IAU
          psumb = g_sum(Atm(n)%domain,sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz),dim=3),&
@@ -1463,7 +1461,6 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
             enddo
          enddo
       enddo
-
    endif
 
    call set_domain ( Atm(mytile)%domain )
@@ -1518,7 +1515,7 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
 #ifdef MULTI_GASES
          q0 = Atm(n)%delp(i,j,k1)*(1.-sum(Atm(n)%q(i,j,k1,1:max(nwat,num_gas)))) + sum(qwat(1:max(nwat,num_gas)))
 #else
-         qt = sum(qwat(1:nwat)) 
+         qt = sum(qwat(1:nwat))
          q0 = Atm(n)%delp(i,j,k1)*(1.-sum(Atm(n)%q(i,j,k1,1:nwat))) + qt 
 #endif
          Atm(n)%delp(i,j,k1) = q0

From d4c23401834035bb3aaf64ce8b3c646e0aa819a7 Mon Sep 17 00:00:00 2001
From: jswhit2 <Jeffrey.S.Whitaker@noaa.gov>
Date: Fri, 7 Feb 2020 13:52:11 +0000
Subject: [PATCH 05/10] remove debug prints

---
 driver/fvGFS/atmosphere.F90 | 18 +-----------------
 1 file changed, 1 insertion(+), 17 deletions(-)

diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90
index 74c342d47..bfef8f290 100644
--- a/driver/fvGFS/atmosphere.F90
+++ b/driver/fvGFS/atmosphere.F90
@@ -183,7 +183,7 @@ module atmosphere_mod
 use fv_nggps_diags_mod, only: fv_nggps_diag_init, fv_nggps_diag, fv_nggps_tavg
 use fv_restart_mod,     only: fv_restart, fv_write_restart
 use fv_timing_mod,      only: timing_on, timing_off
-use fv_mp_mod,          only: switch_current_Atm, is_master
+use fv_mp_mod,          only: switch_current_Atm
 use fv_sg_mod,          only: fv_subgrid_z
 use fv_update_phys_mod, only: fv_update_phys
 use fv_nwp_nudge_mod,   only: fv_nwp_nudge_init, fv_nwp_nudge_end, do_adiabatic_init
@@ -1416,9 +1416,6 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
          qsumb = g_sum(Atm(n)%domain,&
                  sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
                  isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) 
-         if( is_master() ) then
-            print *,'dry surface pressure in IAU interval (before) = ',psumb++Atm%ptop-qsumb
-         endif
       endif
 
 !     IAU increments are in units of 1/sec
@@ -1554,21 +1551,8 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
       qsum = g_sum(Atm(n)%domain,&
              sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
              isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) 
-      if( is_master() ) then
-         print *,'dry surface pressure in IAU interval (after) =',&
-         psum+Atm%ptop-qsum
-      endif
       betad = (psum - (psumb - qsumb))/qsum
-      !if (is_master()) print *,'betad = ',betad
       Atm(n)%q(:,:,:,1:nwat) = betad*Atm(n)%q(:,:,:,1:nwat)
-      ! global mean total water after adjustment
-      !qsum = g_sum(Atm(n)%domain,&
-      !       sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
-      !       isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) 
-      !if( is_master() ) then
-      !   print *,'dry surface pressure in IAU interval (after 2) =',&
-      !   psum+Atm%ptop-qsum
-      !endif
    endif
 
    call timing_off('GFS_TENDENCIES')

From 7a03030c627db5654d4bfa2a7ac5957e41643629 Mon Sep 17 00:00:00 2001
From: jswhit2 <Jeffrey.S.Whitaker@noaa.gov>
Date: Fri, 7 Feb 2020 13:56:47 +0000
Subject: [PATCH 06/10] add drymassfixer to IAU_data structure

---
 tools/fv_iau_mod.F90 | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/tools/fv_iau_mod.F90 b/tools/fv_iau_mod.F90
index e9b961cbb..bf2a73ea1 100644
--- a/tools/fv_iau_mod.F90
+++ b/tools/fv_iau_mod.F90
@@ -93,6 +93,7 @@ module fv_iau_mod
     real,allocatable :: delz_inc(:,:,:)
     real,allocatable :: tracer_inc(:,:,:,:)
     logical          :: in_interval = .false.
+    logical          :: drymassfixer = .false.
   end type iau_external_data_type
   type iau_state_type
       type(iau_internal_data_type):: inc1
@@ -280,6 +281,7 @@ subroutine IAU_initialize (IPD_Control, IAU_Data,Init_parm)
        call read_iau_forcing(IPD_Control,iau_state%inc2,'INPUT/'//trim(IPD_Control%iau_inc_files(2)))
     endif
 !   print*,'in IAU init',dt,rdt
+    IAU_data%drymassfixer = IPD_control%iau_drymassfixer
 
 end subroutine IAU_initialize
 

From 15abb3d9030a66078f2f07c14f30eac1898c8205 Mon Sep 17 00:00:00 2001
From: jswhit2 <Jeffrey.S.Whitaker@noaa.gov>
Date: Sat, 8 Feb 2020 14:58:57 +0000
Subject: [PATCH 07/10] use area_64 in g_sum call

---
 driver/fvGFS/atmosphere.F90 | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90
index bfef8f290..685146e07 100644
--- a/driver/fvGFS/atmosphere.F90
+++ b/driver/fvGFS/atmosphere.F90
@@ -1412,10 +1412,10 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
       if (IAU_Data%drymassfixer) then
          ! global mean total pressure and water before IAU
          psumb = g_sum(Atm(n)%domain,sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz),dim=3),&
-                 isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) 
+                 isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1) 
          qsumb = g_sum(Atm(n)%domain,&
                  sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
-                 isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) 
+                 isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1) 
       endif
 
 !     IAU increments are in units of 1/sec
@@ -1546,11 +1546,11 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
    if (IAU_Data%in_interval .and. IAU_data%drymassfixer) then
       ! global mean total pressure
       psum = g_sum(Atm(n)%domain,sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz),dim=3),&
-             isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) 
+             isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1) 
       ! global mean total water (before adjustment)
       qsum = g_sum(Atm(n)%domain,&
              sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
-             isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) 
+             isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1) 
       betad = (psum - (psumb - qsumb))/qsum
       Atm(n)%q(:,:,:,1:nwat) = betad*Atm(n)%q(:,:,:,1:nwat)
    endif

From 0c3016190c2851d3db54a7878279ec2869579f1b Mon Sep 17 00:00:00 2001
From: jswhit2 <Jeffrey.S.Whitaker@noaa.gov>
Date: Sat, 8 Feb 2020 19:37:31 +0000
Subject: [PATCH 08/10] add debug prints

---
 driver/fvGFS/atmosphere.F90 | 14 +++++++++++++-
 1 file changed, 13 insertions(+), 1 deletion(-)

diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90
index 685146e07..a044895dc 100644
--- a/driver/fvGFS/atmosphere.F90
+++ b/driver/fvGFS/atmosphere.F90
@@ -183,7 +183,7 @@ module atmosphere_mod
 use fv_nggps_diags_mod, only: fv_nggps_diag_init, fv_nggps_diag, fv_nggps_tavg
 use fv_restart_mod,     only: fv_restart, fv_write_restart
 use fv_timing_mod,      only: timing_on, timing_off
-use fv_mp_mod,          only: switch_current_Atm
+use fv_mp_mod,          only: switch_current_Atm, is_master
 use fv_sg_mod,          only: fv_subgrid_z
 use fv_update_phys_mod, only: fv_update_phys
 use fv_nwp_nudge_mod,   only: fv_nwp_nudge_init, fv_nwp_nudge_end, do_adiabatic_init
@@ -1416,6 +1416,9 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
          qsumb = g_sum(Atm(n)%domain,&
                  sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
                  isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1) 
+         if (is_master()) then
+           print *,'dry ps before IAU',psumb+Atm(n)%ptop-qsumb
+         endif
       endif
 
 !     IAU increments are in units of 1/sec
@@ -1552,7 +1555,16 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
              sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
              isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1) 
       betad = (psum - (psumb - qsumb))/qsum
+      if (is_master()) then
+        print *,'dry ps after IAU+physics',psum+Atm(n)%ptop-qsum
+      endif
       Atm(n)%q(:,:,:,1:nwat) = betad*Atm(n)%q(:,:,:,1:nwat)
+      qsum = g_sum(Atm(n)%domain,&
+             sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
+             isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1) 
+      if (is_master()) then
+        print *,'dry ps after iau_drymassfixer',psum+Atm(n)%ptop-qsum
+      endif
    endif
 
    call timing_off('GFS_TENDENCIES')

From 0e30846b5ddacdc8585a288bf8d38dfd30df6a12 Mon Sep 17 00:00:00 2001
From: jswhit2 <Jeffrey.S.Whitaker@noaa.gov>
Date: Sun, 9 Feb 2020 03:09:33 +0000
Subject: [PATCH 09/10] remove one debug print

---
 driver/fvGFS/atmosphere.F90 | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90
index a044895dc..39ba8c4c8 100644
--- a/driver/fvGFS/atmosphere.F90
+++ b/driver/fvGFS/atmosphere.F90
@@ -1417,7 +1417,7 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
                  sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
                  isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1) 
          if (is_master()) then
-           print *,'dry ps before IAU',psumb+Atm(n)%ptop-qsumb
+           print *,'dry ps before IAU/physics',psumb+Atm(n)%ptop-qsumb
          endif
       endif
 
@@ -1556,15 +1556,15 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
              isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1) 
       betad = (psum - (psumb - qsumb))/qsum
       if (is_master()) then
-        print *,'dry ps after IAU+physics',psum+Atm(n)%ptop-qsum
+        print *,'dry ps after IAU/physics',psum+Atm(n)%ptop-qsum
       endif
       Atm(n)%q(:,:,:,1:nwat) = betad*Atm(n)%q(:,:,:,1:nwat)
-      qsum = g_sum(Atm(n)%domain,&
-             sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
-             isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1) 
-      if (is_master()) then
-        print *,'dry ps after iau_drymassfixer',psum+Atm(n)%ptop-qsum
-      endif
+      !qsum = g_sum(Atm(n)%domain,&
+      !       sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
+      !       isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1) 
+      !if (is_master()) then
+      !  print *,'dry ps after iau_drymassfixer',psum+Atm(n)%ptop-qsum
+      !endif
    endif
 
    call timing_off('GFS_TENDENCIES')

From da9077d09a6e294d5b383710c6c5fdddecff46e6 Mon Sep 17 00:00:00 2001
From: "Jun.Wang" <Jun.Wang@noaa.gov>
Date: Mon, 17 Feb 2020 21:36:06 +0000
Subject: [PATCH 10/10] add reproduce in g_sum call

---
 driver/fvGFS/atmosphere.F90 | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90
index 39ba8c4c8..0fc5ba08e 100644
--- a/driver/fvGFS/atmosphere.F90
+++ b/driver/fvGFS/atmosphere.F90
@@ -1412,10 +1412,10 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
       if (IAU_Data%drymassfixer) then
          ! global mean total pressure and water before IAU
          psumb = g_sum(Atm(n)%domain,sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz),dim=3),&
-                 isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1) 
+                 isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1,reproduce=.true.) 
          qsumb = g_sum(Atm(n)%domain,&
                  sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
-                 isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1) 
+                 isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1,reproduce=.true.) 
          if (is_master()) then
            print *,'dry ps before IAU/physics',psumb+Atm(n)%ptop-qsumb
          endif
@@ -1549,11 +1549,11 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
    if (IAU_Data%in_interval .and. IAU_data%drymassfixer) then
       ! global mean total pressure
       psum = g_sum(Atm(n)%domain,sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz),dim=3),&
-             isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1) 
+             isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1,reproduce=.true.) 
       ! global mean total water (before adjustment)
       qsum = g_sum(Atm(n)%domain,&
              sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
-             isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1) 
+             isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1,reproduce=.true.) 
       betad = (psum - (psumb - qsumb))/qsum
       if (is_master()) then
         print *,'dry ps after IAU/physics',psum+Atm(n)%ptop-qsum