From 25d8b7ba407f28f591960e261aa00df23fce7e36 Mon Sep 17 00:00:00 2001
From: Wen Meng <wen.meng@noaa.gov>
Date: Mon, 12 Jun 2023 22:42:50 +0000
Subject: [PATCH 1/9] Update upp revision and clean up inline post read
 interface.

---
 io/post_fv3.F90 | 1134 +++++++++++++++++++++--------------------------
 upp             |    2 +-
 2 files changed, 509 insertions(+), 627 deletions(-)

diff --git a/io/post_fv3.F90 b/io/post_fv3.F90
index 7cc6ab45e..c8c6fa5c3 100644
--- a/io/post_fv3.F90
+++ b/io/post_fv3.F90
@@ -187,7 +187,7 @@ subroutine post_run_fv3(wrt_int_state,grid_id,mype,mpicomp,lead_write, &
         call set_outflds(kth,th,kpv,pv)
         if(allocated(datapd))deallocate(datapd)
         allocate(datapd(ite-its+1,jte-jts+1,nrecout+100))
-!$omp parallel do default(none),private(i,j,k),shared(nrecout,jend,jsta,im,datapd,ista,iend)
+!$omp parallel do default(none),private(i,j,k),shared(nrecout,jend,jsta,datapd,ista,iend)
         do k=1,nrecout+100
           do j=1,jend+1-jsta
             do i=1,iend+1-ista
@@ -494,6 +494,8 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
 !     Apr 2022    W. Meng      Unify set_postvars_gfs and
 !                               set_postvars_regional to set_postvars_fv3
 !     Apr 2023    W. Meng      Sync RRFS and GFS changes from off-line post
+!     Jun 2023    W. Meng      Remove duplicate initialization;
+!                              relocate computation of aerosol fields
 !
 !-----------------------------------------------------------------------
 !*** set up post fields from nmint_state
@@ -670,7 +672,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
         bk5(i) = wrt_int_state%bk(i)
       enddo
 
-!$omp parallel do default(none) private(i,j) shared(jsta,jend,im,f,gdlat,ista,iend)
+!$omp parallel do default(none) private(i,j) shared(jsta,jend,f,gdlat,ista,iend)
       do j=jsta,jend
         do i=ista,iend
           f(I,J) = 1.454441e-4*sin(gdlat(i,j)*dtr)   ! 2*omeg*sin(phi)
@@ -679,18 +681,6 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
 !
       pt    = ak5(1)
 
-! GFS does not have surface specific humidity
-!                   inst sensible heat flux
-!                   inst latent heat flux
-!$omp parallel do default(none),private(i,j),shared(jsta,jend,im,spval,qs,twbs,qwbs,ths,ista,iend)
-      do j=jsta,jend
-        do i=ista,iend
-          qs(i,j) = SPVAL
-          twbs(i,j) = SPVAL
-          qwbs(i,j) = SPVAL
-        enddo
-      enddo
-
 ! GFS set up DT to compute accumulated fields, set it to one
       dtq2 = wrt_int_state%dtp
       nphs = 2.
@@ -706,49 +696,15 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
         allocate(snacc_land(ista:iend,jsta:jend))
       endif
 
-!
-! GFS does not have convective cloud efficiency
-!                   similated precip
-!                   10 m theta
-!                   10 m humidity
-!                   snow free albedo
-!$omp parallel do default(none), private(i,j), shared(jsta,jend,im,spval,ista,iend), &
-!$omp& shared(cldefi,lspa,th10,q10,albase)
-      do j=jsta,jend
-        do i=ista,iend
-          cldefi(i,j) = SPVAL
-          lspa(i,j) = SPVAL
-          th10(i,j) = SPVAL
-          q10(i,j) = SPVAL
-          albase(i,j) = SPVAL
-        enddo
-      enddo
-
-! GFS does not have convective precip
-!$omp parallel do default(none) private(i,j) shared(jsta,jend,im,cprate,ista,iend)
-      do j=jsta,jend
-        do i=ista,iend
-          cprate(i,j) = 0.
-        enddo
-      enddo
+      !wm
+!! GFS does not have convective precip
+!!$omp parallel do default(none) private(i,j) shared(jsta,jend,cprate,ista,iend)
+!      do j=jsta,jend
+!        do i=ista,iend
+!          cprate(i,j) = 0.
+!        enddo
+!      enddo
 
-! GFS probably does not use zenith angle, czen, czmean
-!                       inst surface outgoing longwave, radot
-!                       inst cloud fraction for high, middle, and low cloud,
-!                            cfrach
-!                       inst ground heat flux, grnflx
-!$omp parallel do default(none) private(i,j) shared(jsta,jend,im,spval,ista,iend), &
-!$omp& shared(czen,czmean,radot,cfrach,cfracl,cfracm,grnflx)
-      do j=jsta,jend
-        do i=ista,iend
-          czen(i,j)   = SPVAL
-          czmean(i,j) = SPVAL
-          cfrach(i,j) = SPVAL
-          cfracl(i,j) = SPVAL
-          cfracm(i,j) = SPVAL
-          grnflx(i,j) = SPVAL
-        enddo
-      enddo
 !
 ! GFS doesn not yet output soil layer thickness, assign SLDPTH to be the same as nam
       sldpth(1) = 0.10
@@ -766,17 +722,12 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
 !                     inst outgoing sfc shortwave, rswout
 !                     snow phase change heat flux, snopcx
 ! GFS does not use total momentum flux,sfcuvx
-!$omp parallel do default(none),private(i,j),shared(jsta,jend,im,spval,ista,iend), &
-!$omp& shared(acfrcv,ncfrcv,acfrst,ncfrst,bgroff,rswin,rswinc,rswout,snopcx,sfcuvx,&
-!$omp& ltg1_max,ltg2_max,ltg3_max)
+!$omp parallel do default(none),private(i,j),shared(jsta,jend,spval,ista,iend), &
+!$omp& shared(acfrcv,ncfrcv,acfrst,ncfrst,bgroff,rswinc)
       do j=jsta,jend
         do i=ista,iend
-          acfrcv(i,j) = spval
           ncfrcv(i,j) = 1.0
-          acfrst(i,j) = spval
           ncfrst(i,j) = 1.0
-          bgroff(i,j) = spval
-          rswinc(i,j) = spval
         enddo
       enddo
 
@@ -787,77 +738,11 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
 ! GFS surface flux has been averaged, set  ASRFC to 1
       asrfc = 1.0
 
-! GFS does not have temperature tendency due to long wave radiation
-!                   temperature tendency due to short wave radiation
-!                   temperature tendency due to latent heating from convection
-!                   temperature tendency due to latent heating from grid scale
-      do l=1,lm
-!$omp parallel do default(none),private(i,j),shared(jsta_2l,jend_2u,im,spval,l,ista_2l,iend_2u), &
-!$omp& shared(rlwtt,rswtt,tcucn,tcucns,train)
-        do j=jsta_2l,jend_2u
-          do i=ista_2l,iend_2u
-            rlwtt(i,j,l) = spval
-            rswtt(i,j,l)  = spval
-            tcucn(i,j,l)  = spval
-            tcucns(i,j,l) = spval
-            train(i,j,l)  = spval
-          enddo
-        enddo
-      enddo
-
 ! set avrain to 1
       avrain = 1.0
       avcnvc = 1.0
       theat  = 6.0 ! just in case GFS decides to output T tendency
 
-! GFS does not have temperature tendency due to latent heating from grid scale
-      train  = spval
-
-! GFS does not have soil moisture availability, smstav
-!                   accumulated surface evaporatio, sfcevp
-!                   averaged accumulated snow, acsnow
-!                   snow melt,acsnom
-!                   humidity at roughness length, qz0
-!                   u at roughness length, uz0
-!                   v at roughness length, vz0
-!                   shelter rh max, maxrhshltr
-!                   shelter rh min, minrhshltr
-!$omp parallel do default(none),private(i,j),shared(jsta_2l,jend_2u,im,spval,ista_2l,iend_2u), &
-!$omp& shared(sfcevp,acsnom,qz0,uz0,vz0,maxrhshltr,minrhshltr)
-      do j=jsta_2l,jend_2u
-        do i=ista_2l,iend_2u
-          sfcevp(i,j) = spval
-          acsnom(i,j) = spval
-          qz0(i,j)    = spval
-          uz0(i,j)    = spval
-          vz0(i,j)    = spval
-        enddo
-      enddo
-
-! GFS does not have mixing length,el_pbl
-!                   exchange coefficient, exch_h
-      do l=1,lm
-!$omp parallel do default(none),private(i,j),shared(jsta_2l,jend_2u,im,l,spval,el_pbl,exch_h,ista_2l,iend_2u)
-        do j=jsta_2l,jend_2u
-          do i=ista_2l,iend_2u
-            el_pbl(i,j,l) = spval
-            exch_h(i,j,l) = spval
-          enddo
-        enddo
-      enddo
-
-! GFS does not have deep convective cloud top and bottom fields
-!$omp parallel do default(none),private(i,j),shared(jsta_2l,jend_2u,im,spval,ista_2l,iend_2u), &
-!$omp& shared(htopd,hbotd,htops,hbots,cuppt)
-      do j=jsta_2l,jend_2u
-        do i=ista_2l,iend_2u
-          htopd(i,j) = SPVAL
-          hbotd(i,j) = SPVAL
-          htops(i,j) = SPVAL
-          hbots(i,j) = SPVAL
-          cuppt(i,j) = SPVAL
-        enddo
-      enddo
 !
 ! get inital date
       sdat(1)  = wrt_int_state%idate(2)   !month
@@ -2457,9 +2342,9 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
             endif !end rdaod
 
             if ((gocart_on .or. gccpp_on) .and. d2d_chem) then
+
               do K = 1, nbin_du
                 write(VarName, '(A,I3.3)') 'duem', k
-
                 if(trim(fieldname)==VarName) then
                   !$omp parallel do default(none) private(i,j,K) shared(jsta,jend,ista,iend,spval,duem,arrayr42d,fillvalue)
                   do j=jsta,jend
@@ -2595,7 +2480,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
                     enddo
                   enddo
                  endif
-               enddo
+              enddo
 
               do K = 1, nbin_ss
                 if ( K == 1) VarName='seas1wt'
@@ -2613,7 +2498,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
                     enddo
                   enddo
                  endif
-               enddo
+              enddo
 
               do K = 1, nbin_ss
                 if ( K == 1) VarName='seas1wtc'
@@ -2631,7 +2516,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
                     enddo
                   enddo
                  endif
-               enddo
+              enddo
 
               do K = 1, nbin_bc
                 if ( K == 1) VarName='bceman'
@@ -2646,7 +2531,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
                     enddo
                   enddo
                  endif
-               enddo
+              enddo
 
               do K = 1, nbin_bc
                 if ( K == 1) VarName='bc1sd'
@@ -2784,15 +2669,15 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
               enddo
 
 
-                if(trim(fieldname)=='maod') then
-                  !$omp parallel do default(none) private(i,j,K) shared(jsta,jend,ista,iend,spval,maod,arrayr42d,fillvalue)
-                  do j=jsta,jend
-                    do i=ista, iend
-                      maod(i,j) = arrayr42d(i,j)
-                      if( abs(arrayr42d(i,j)-fillValue) < small) maod(i,j) = spval
-                    enddo
+              if(trim(fieldname)=='maod') then
+                !$omp parallel do default(none) private(i,j,K) shared(jsta,jend,ista,iend,spval,maod,arrayr42d,fillvalue)
+                do j=jsta,jend
+                  do i=ista, iend
+                    maod(i,j) = arrayr42d(i,j)
+                    if( abs(arrayr42d(i,j)-fillValue) < small) maod(i,j) = spval
                   enddo
-                endif
+                enddo
+              endif
 
             endif !end gocart_on
 
@@ -3775,7 +3660,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
                 enddo
               enddo
             endif
-              if(mype==0) print *,'in gfs_post, get ref_10cm=',maxval(ref_10cm), minval(ref_10cm)
+!              if(mype==0) print *,'in gfs_post, get ref_10cm=',maxval(ref_10cm), minval(ref_10cm)
             else
             if(trim(fieldname)=='refl_10cm') then
               !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,ref_10cm,arrayr43d,fillvalue,spval)
@@ -3952,467 +3837,302 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
               enddo
             endif
 
-! read chemical fields
-      if(gocart_on .or. gccpp_on .or. nasa_on) then
-
-        if(trim(fieldname)=='dust1') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,dust,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                dust(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) dust(i,j,l,1) = spval
-              enddo
-            enddo
-          enddo
-        endif
+            ! read chemical fields
+            if(gocart_on .or. gccpp_on .or. nasa_on) then
 
-        if(trim(fieldname)=='dust2') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,dust,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                dust(i,j,l,2) = max(arrayr43d(i,j,l),0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) dust(i,j,l,2) = spval
-              enddo
-            enddo
-          enddo
-        endif
+              if(trim(fieldname)=='dust1') then
+              !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,dust,arrayr43d,fillvalue,spval)
+                do l=1,lm
+                  do j=jsta,jend
+                    do i=ista, iend
+                      dust(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
+                      if(abs(arrayr43d(i,j,l)-fillvalue) < small) dust(i,j,l,1) = spval
+                    enddo
+                  enddo
+                enddo
+              endif
 
-        if(trim(fieldname)=='dust3') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,dust,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                dust(i,j,l,3) = max(arrayr43d(i,j,l), 0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) dust(i,j,l,3) = spval
-              enddo
-            enddo
-          enddo
-        endif
+              if(trim(fieldname)=='dust2') then
+              !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,dust,arrayr43d,fillvalue,spval)
+                do l=1,lm
+                  do j=jsta,jend
+                    do i=ista, iend
+                      dust(i,j,l,2) = max(arrayr43d(i,j,l),0.0)
+                      if(abs(arrayr43d(i,j,l)-fillvalue) < small) dust(i,j,l,2) = spval
+                    enddo
+                  enddo
+                enddo
+              endif
 
-        if(trim(fieldname)=='dust4') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,dust,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                dust(i,j,l,4) = max(arrayr43d(i,j,l), 0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) dust(i,j,l,4) = spval
-              enddo
-            enddo
-          enddo
-        endif
+              if(trim(fieldname)=='dust3') then
+              !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,dust,arrayr43d,fillvalue,spval)
+                do l=1,lm
+                  do j=jsta,jend
+                    do i=ista, iend
+                      dust(i,j,l,3) = max(arrayr43d(i,j,l), 0.0)
+                      if(abs(arrayr43d(i,j,l)-fillvalue) < small) dust(i,j,l,3) = spval
+                    enddo
+                  enddo
+                enddo
+              endif
 
-        if(trim(fieldname)=='dust5') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,dust,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                dust(i,j,l,5) = max(arrayr43d(i,j,l), 0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) dust(i,j,l,5) = spval
-              enddo
-            enddo
-          enddo
-        endif
+              if(trim(fieldname)=='dust4') then
+              !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,dust,arrayr43d,fillvalue,spval)
+                do l=1,lm
+                  do j=jsta,jend
+                    do i=ista, iend
+                      dust(i,j,l,4) = max(arrayr43d(i,j,l), 0.0)
+                      if(abs(arrayr43d(i,j,l)-fillvalue) < small) dust(i,j,l,4) = spval
+                    enddo
+                  enddo
+                enddo
+              endif
 
-        dustcb=0.0
-        dustallcb=0.0
-        do l=1,lm
-          do j=jsta,jend
-            do i=ista, iend
-              if(dust(i,j,l,1)<spval.and.dust(i,j,l,2)<spval.and. &
-                 dust(i,j,l,3)<spval.and.dust(i,j,l,4)<spval.and. &
-                 dust(i,j,l,5)<spval) then
-                dustcb(i,j)=dustcb(i,j)+ &
-                  (dust(i,j,l,1)+0.38*dust(i,j,l,2))* &
-                  dpres(i,j,l)/grav
-                dustallcb(i,j)=dustallcb(i,j)+ &
-                  (dust(i,j,l,1)+dust(i,j,l,2)+ &
-                  dust(i,j,l,3)+0.74*dust(i,j,l,4))* &
-                  dpres(i,j,l)/grav
+              if(trim(fieldname)=='dust5') then
+              !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,dust,arrayr43d,fillvalue,spval)
+                do l=1,lm
+                  do j=jsta,jend
+                    do i=ista, iend
+                      dust(i,j,l,5) = max(arrayr43d(i,j,l), 0.0)
+                      if(abs(arrayr43d(i,j,l)-fillvalue) < small) dust(i,j,l,5) = spval
+                    enddo
+                  enddo
+                enddo
               endif
-            enddo
-          enddo
-        enddo
 
-        if(trim(fieldname)=='seas1') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,salt,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                salt(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) salt(i,j,l,1) = spval
-              enddo
-            enddo
-          enddo
-        endif
+              if(trim(fieldname)=='seas1') then
+              !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,salt,arrayr43d,fillvalue,spval)
+                do l=1,lm
+                  do j=jsta,jend
+                    do i=ista, iend
+                      salt(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
+                      if(abs(arrayr43d(i,j,l)-fillvalue) < small) salt(i,j,l,1) = spval
+                    enddo
+                  enddo
+                enddo
+              endif
 
-        if(trim(fieldname)=='seas2') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,salt,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                salt(i,j,l,2) = max(arrayr43d(i,j,l),0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) salt(i,j,l,2) = spval
-              enddo
-            enddo
-          enddo
-        endif
+              if(trim(fieldname)=='seas2') then
+              !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,salt,arrayr43d,fillvalue,spval)
+                do l=1,lm
+                  do j=jsta,jend
+                    do i=ista, iend
+                      salt(i,j,l,2) = max(arrayr43d(i,j,l),0.0)
+                      if(abs(arrayr43d(i,j,l)-fillvalue) < small) salt(i,j,l,2) = spval
+                    enddo
+                  enddo
+                enddo
+              endif
 
+              if(trim(fieldname)=='seas3') then
+              !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,salt,arrayr43d,fillvalue,spval)
+                do l=1,lm
+                  do j=jsta,jend
+                    do i=ista, iend
+                      salt(i,j,l,3) = max(arrayr43d(i,j,l),0.0)
+                      if(abs(arrayr43d(i,j,l)-fillvalue) < small) salt(i,j,l,3) = spval
+                    enddo
+                  enddo
+                enddo
+              endif
 
-        if(trim(fieldname)=='seas3') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,salt,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                salt(i,j,l,3) = max(arrayr43d(i,j,l),0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) salt(i,j,l,3) = spval
-              enddo
-            enddo
-          enddo
-        endif
+              if(trim(fieldname)=='seas4') then
+              !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,salt,arrayr43d,fillvalue,spval)
+                do l=1,lm
+                  do j=jsta,jend
+                    do i=ista, iend
+                      salt(i,j,l,4) = max(arrayr43d(i,j,l),0.0)
+                      if(abs(arrayr43d(i,j,l)-fillvalue) < small) salt(i,j,l,4) = spval
+                    enddo
+                  enddo
+                enddo
+              endif
 
-        if(trim(fieldname)=='seas4') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,salt,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                salt(i,j,l,4) = max(arrayr43d(i,j,l),0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) salt(i,j,l,4) = spval
-              enddo
-            enddo
-          enddo
-        endif
+              if(trim(fieldname)=='seas5') then
+              !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,salt,arrayr43d,fillvalue,spval)
+                do l=1,lm
+                  do j=jsta,jend
+                    do i=ista, iend
+                      salt(i,j,l,5) = max(arrayr43d(i,j,l),0.0)
+                      if(abs(arrayr43d(i,j,l)-fillvalue) < small) salt(i,j,l,5) = spval
+                    enddo
+                  enddo
+                enddo
+              endif
 
-        if(trim(fieldname)=='seas5') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,salt,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                salt(i,j,l,5) = max(arrayr43d(i,j,l),0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) salt(i,j,l,5) = spval
-              enddo
-            enddo
-          enddo
-        endif
+              if(trim(fieldname)=='bc1') then
+              !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,soot,arrayr43d,fillvalue,spval)
+                do l=1,lm
+                  do j=jsta,jend
+                    do i=ista, iend
+                      soot(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
+                      if(abs(arrayr43d(i,j,l)-fillvalue) < small) soot(i,j,l,1) = spval
+                    enddo
+                  enddo
+                enddo
+              endif
 
-        sscb=0.0
-        ssallcb=0.0
-        do l=1,lm
-          do j=jsta,jend
-            do i=ista,iend
-              if(salt(i,j,l,1)<spval.and.salt(i,j,l,2)<spval.and. &
-                 salt(i,j,l,3)<spval.and.salt(i,j,l,4)<spval.and. & 
-                 salt(i,j,l,5)<spval) then
-            sscb(i,j)=sscb(i,j)+ &
-         (salt(i,j,l,1)+salt(i,j,l,2)+0.83*salt(i,j,l,3))*  &
-           dpres(i,j,l)/grav
+              if(trim(fieldname)=='bc2') then
+              !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,soot,arrayr43d,fillvalue,spval)
+                do l=1,lm
+                  do j=jsta,jend
+                    do i=ista, iend
+                      soot(i,j,l,2) = max(arrayr43d(i,j,l),0.0)
+                      if(abs(arrayr43d(i,j,l)-fillvalue) < small) soot(i,j,l,2) = spval
+                    enddo
+                  enddo
+                enddo
+              endif
 
+              if(trim(fieldname)=='oc1') then
+              !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,waso,arrayr43d,fillvalue,spval)
+                do l=1,lm
+                  do j=jsta,jend
+                    do i=ista, iend
+                      waso(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
+                      if(abs(arrayr43d(i,j,l)-fillvalue) < small) waso(i,j,l,1) = spval
+                    enddo
+                  enddo
+                enddo
+              endif
 
-          ssallcb(i,j)=ssallcb(i,j)+ &
-         (salt(i,j,l,1)+salt(i,j,l,2)+salt(i,j,l,3)+salt(i,j,l,4))* &
-           dpres(i,j,l)/grav
+              if(trim(fieldname)=='oc2') then
+              !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,waso,arrayr43d,fillvalue,spval)
+                do l=1,lm
+                  do j=jsta,jend
+                    do i=ista, iend
+                      waso(i,j,l,2) = max(arrayr43d(i,j,l),0.0)
+                      if(abs(arrayr43d(i,j,l)-fillvalue) < small) waso(i,j,l,2) = spval
+                    enddo
+                  enddo
+                enddo
               endif
-            enddo
-          enddo
-        end do ! do loop for l       
 
-        if(trim(fieldname)=='bc1') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,soot,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                soot(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) soot(i,j,l,1) = spval
-              enddo
-            enddo
-          enddo
-        endif
+              if (gocart_on .or. gccpp_on) then
+                if(trim(fieldname)=='sulf') then
+                 !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,suso,arrayr43d,fillvalue,spval)
+                  do l=1,lm
+                    do j=jsta,jend
+                      do i=ista, iend
+                        suso(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
+                        if(abs(arrayr43d(i,j,l)-fillvalue) < small) suso(i,j,l,1) = spval
+                      enddo
+                    enddo
+                  enddo
+                endif
 
-        if(trim(fieldname)=='bc2') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,soot,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                soot(i,j,l,2) = max(arrayr43d(i,j,l),0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) soot(i,j,l,2) = spval
-              enddo
-            enddo
-          enddo
-        endif
+                if(trim(fieldname)=='pp25') then
+                 !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,pp25,arrayr43d,fillvalue,spval)
+                  do l=1,lm
+                    do j=jsta,jend
+                      do i=ista, iend
+                        pp25(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
+                        if(abs(arrayr43d(i,j,l)-fillvalue) < small) pp25(i,j,l,1) = spval
+                      enddo
+                    enddo
+                  enddo
+                endif
 
-        if(trim(fieldname)=='oc1') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,waso,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                waso(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) waso(i,j,l,1) = spval
-              enddo
-            enddo
-          enddo
-        endif
-
-        if(trim(fieldname)=='oc2') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,waso,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                waso(i,j,l,2) = max(arrayr43d(i,j,l),0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) waso(i,j,l,2) = spval
-              enddo
-            enddo
-          enddo
-        endif
-
-        bccb=0.0
-        occb=0.0
-        do l=1,lm
-          do j=jsta,jend
-            do i=ista,iend
-              if(soot(i,j,l,1)<spval.and.soot(i,j,l,2)<spval)then
-               bccb(i,j)=bccb(i,j)+(soot(i,j,l,1)+soot(i,j,l,2))* &
-               dpres(i,j,l)/grav
-              endif
-
-              if(waso(i,j,l,1)<spval.and.waso(i,j,l,2)<spval)then
-               occb(i,j)=occb(i,j)+ (waso(i,j,l,1)+waso(i,j,l,2))* &
-               dpres(i,j,l)/grav
-              endif
-            enddo
-          enddo
-        end do ! do loop for l
-
-
-
-       if (gocart_on .or. gccpp_on) then
-        if(trim(fieldname)=='sulf') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,suso,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                suso(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) suso(i,j,l,1) = spval
-              enddo
-            enddo
-          enddo
-        endif
-
-        if(trim(fieldname)=='pp25') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,pp25,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                pp25(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) pp25(i,j,l,1) = spval
-              enddo
-            enddo
-          enddo
-        endif
-
-        if(trim(fieldname)=='pp10') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,pp10,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                pp10(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) pp10(i,j,l,1) = spval
-              enddo
-            enddo
-          enddo
-        endif
-
-       else if (nasa_on) then 
-        if(trim(fieldname)=='so4') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,suso,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                suso(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) suso(i,j,l,1) = spval
-              enddo
-            enddo
-          enddo
-        endif
-
-        if(trim(fieldname)=='no3an1') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,no3,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                no3(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) no3(i,j,l,1) = spval
-              enddo
-            enddo
-          enddo
-        endif
-
-        if(trim(fieldname)=='no3an2') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,no3,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                no3(i,j,l,2) = max(arrayr43d(i,j,l),0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) no3(i,j,l,2) = spval
-              enddo
-            enddo
-          enddo
-        endif
-
-        if(trim(fieldname)=='no3an3') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,no3,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                no3(i,j,l,3) = max(arrayr43d(i,j,l),0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) no3(i,j,l,3) = spval
-              enddo
-            enddo
-          enddo
-        endif
-
-        if(trim(fieldname)=='nh4a') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,nh4,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                nh4(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) nh4(i,j,l,1) = spval
-              enddo
-            enddo
-          enddo
-        endif
-
-       no3cb=0.0
-       nh4cb=0.0
-        do l=1,lm
-          do j=jsta,jend
-            do i=ista,iend
-              if(no3(i,j,l,1)<spval .and. no3(i,j,l,2)<spval .and. &
-                 no3(i,j,l,3)<spval) then
-              no3cb(i,j)=no3cb(i,j)+ (no3(i,j,l,1)+no3(i,j,l,2)+ &
-              no3(i,j,l,3) ) * dpres(i,j,l)/grav
-              endif
-              if(nh4(i,j,l,1)<spval)then
-              nh4cb(i,j)=nh4cb(i,j)+ nh4(i,j,l,1)* &
-              dpres(i,j,l)/grav
-              endif
-            enddo
-          enddo
-        end do ! do loop for l
-
-     if(trim(fieldname)=='pm25') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,pp25,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                pp25(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) pp25(i,j,l,1) = spval
-              enddo
-            enddo
-          enddo
-        endif
-
-        if(trim(fieldname)=='pm10') then
-          !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,pp10,arrayr43d,fillvalue,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                pp10(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
-                if(abs(arrayr43d(i,j,l)-fillvalue) < small) pp10(i,j,l,1) = spval
-              enddo
-            enddo
-          enddo
-        endif
-
-       endif !nasa_on
-
-       sulfcb=0.0
-       pp25cb=0.0
-       pp10cb=0.0
-        do l=1,lm
-          do j=jsta,jend 
-            do i=ista,iend
-              if(suso(i,j,l,1)<spval)then
-              sulfcb(i,j)=sulfcb(i,j)+ suso(i,j,l,1)* &
-              dpres(i,j,l)/grav
-              endif
-              if(pp25(i,j,l,1)<spval)then
-              pp25cb(i,j)=pp25cb(i,j)+ pp25(i,j,l,1)* &
-              dpres(i,j,l)/grav
-              endif
-              if(pp10(i,j,l,1)<spval)then
-              pp10cb(i,j)=pp10cb(i,j)+ pp10(i,j,l,1)* &
-              dpres(i,j,l)/grav
-              endif
-            enddo
-          enddo
-        end do ! do loop for l
-
-        l=lm
-        do j=jsta,jend
-          do i=ista,iend
+                if(trim(fieldname)=='pp10') then
+                 !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,pp10,arrayr43d,fillvalue,spval)
+                  do l=1,lm
+                    do j=jsta,jend
+                      do i=ista, iend
+                        pp10(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
+                        if(abs(arrayr43d(i,j,l)-fillvalue) < small) pp10(i,j,l,1) = spval
+                      enddo
+                    enddo
+                  enddo
+                endif
 
-            TV = T(I,J,L) * (H1+D608*MAX(Q(I,J,L),QMIN))
-            RHOMID(I,J,L) = PMID(I,J,L) / (RD*TV)
+              else if (nasa_on) then 
+                if(trim(fieldname)=='so4') then
+                 !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,suso,arrayr43d,fillvalue,spval)
+                  do l=1,lm
+                    do j=jsta,jend
+                      do i=ista, iend
+                        suso(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
+                        if(abs(arrayr43d(i,j,l)-fillvalue) < small) suso(i,j,l,1) = spval
+                      enddo
+                    enddo
+                  enddo
+                endif
 
-            dustcb(i,j) = MAX(dustcb(i,j), 0.0)
-            dustallcb(i,j) = MAX(dustallcb(i,j), 0.0)
-            sscb(i,j) = MAX(sscb(i,j), 0.0)
-            ssallcb(i,j) = MAX(ssallcb(i,j), 0.0)
-            bccb(i,j) = MAX(bccb(i,j), 0.0)
-            occb(i,j) = MAX(occb(i,j), 0.0)
-            sulfcb(i,j) = MAX(sulfcb(i,j), 0.0)
-            if (nasa_on) then
-            no3cb(i,j) = MAX(no3cb(i,j), 0.0)
-            nh4cb(i,j) = MAX(nh4cb(i,j), 0.0)
-            endif
-            pp25cb(i,j) = MAX(pp25cb(i,j), 0.0)
-            pp10cb(i,j) = MAX(pp10cb(i,j), 0.0)
+                if(trim(fieldname)=='no3an1') then
+                 !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,no3,arrayr43d,fillvalue,spval)
+                  do l=1,lm
+                    do j=jsta,jend
+                      do i=ista, iend
+                        no3(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
+                        if(abs(arrayr43d(i,j,l)-fillvalue) < small) no3(i,j,l,1) = spval
+                      enddo
+                    enddo
+                  enddo
+                endif
 
-!      Surface PM25 dust and seasalt      
-           dustpm(i,j)=(dust(i,j,l,1)+0.38*dust(i,j,l,2))*RHOMID(i,j,l) !ug/m3
-           dustpm10(i,j)=(dust(i,j,l,1)+dust(i,j,l,2)+dust(i,j,l,3)+ &
-            0.74*dust(i,j,l,4))*RHOMID(i,j,l) !ug/m3
-           sspm(i,j)=(salt(i,j,l,1)+salt(i,j,l,2)+ &
-            0.83*salt(i,j,l,3))*RHOMID(i,j,l)  !ug/m3 
+                if(trim(fieldname)=='no3an2') then
+                 !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,no3,arrayr43d,fillvalue,spval)
+                  do l=1,lm
+                    do j=jsta,jend
+                      do i=ista, iend
+                        no3(i,j,l,2) = max(arrayr43d(i,j,l),0.0)
+                        if(abs(arrayr43d(i,j,l)-fillvalue) < small) no3(i,j,l,2) = spval
+                      enddo
+                    enddo
+                  enddo
+                endif
 
-           if (gocart_on .or. gccpp_on) then
-             !Surface PM10 concentration
-             dusmass(i,j)=(dust(i,j,l,1)+dust(i,j,l,2)+dust(i,j,l,3)+ &
-              0.74*dust(i,j,l,4)+salt(i,j,l,1)+salt(i,j,l,2)+salt(i,j,l,3)+ &
-              salt(i,j,l,4) + soot(i,j,l,1)+soot(i,j,l,2)+waso(i,j,l,1)+ &
-              waso(i,j,l,2) +suso(i,j,l,1)+pp25(i,j,l,1)+pp10(i,j,l,1)) &
-              *RHOMID(i,j,l)  !ug/m3
-             !Surface PM25 concentration       
-             dusmass25(i,j)=(dust(i,j,l,1)+0.38*dust(i,j,l,2)+ &
-             salt(i,j,l,1)+salt(i,j,l,2)+0.83*salt(i,j,l,3) + &
-             soot(i,j,l,1)+soot(i,j,l,2)+waso(i,j,l,1)+ &
-             waso(i,j,l,2) +suso(i,j,l,1)+pp25(i,j,l,1))*RHOMID(i,j,l)  !ug/m3
+                if(trim(fieldname)=='no3an3') then
+                 !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,no3,arrayr43d,fillvalue,spval)
+                  do l=1,lm
+                    do j=jsta,jend
+                      do i=ista, iend
+                        no3(i,j,l,3) = max(arrayr43d(i,j,l),0.0)
+                        if(abs(arrayr43d(i,j,l)-fillvalue) < small) no3(i,j,l,3) = spval
+                      enddo
+                    enddo
+                  enddo
+                endif
 
-             !PM10 column
-             ducmass(i,j)=dustallcb(i,j)+ssallcb(i,j)+bccb(i,j)+ &
-             occb(i,j)+sulfcb(i,j)+pp25cb(i,j)+pp10cb(i,j)
-             !PM25 column
-             ducmass25(i,j)=dustcb(i,j)+sscb(i,j)+bccb(i,j)+occb(i,j) &
-             +sulfcb(i,j)+pp25cb(i,j)
+                if(trim(fieldname)=='nh4a') then
+                 !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,nh4,arrayr43d,fillvalue,spval)
+                  do l=1,lm
+                    do j=jsta,jend
+                      do i=ista, iend
+                        nh4(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
+                        if(abs(arrayr43d(i,j,l)-fillvalue) < small) nh4(i,j,l,1) = spval
+                      enddo
+                    enddo
+                  enddo
+                endif
 
-           elseif (nasa_on) then
-             !Surface PM10 concentration
-             dusmass(i,j)=pp10(i,j,l,1)*RHOMID(i,j,l)  !ug/m3
-             !Surface PM25 concentration       
-             dusmass25(i,j)=pp25(i,j,l,1)*RHOMID(i,j,l)  !ug/m3
+                if(trim(fieldname)=='pm25') then
+                 !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,pp25,arrayr43d,fillvalue,spval)
+                  do l=1,lm
+                    do j=jsta,jend
+                      do i=ista, iend
+                        pp25(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
+                        if(abs(arrayr43d(i,j,l)-fillvalue) < small) pp25(i,j,l,1) = spval
+                      enddo
+                    enddo
+                  enddo
+                endif
 
-             !PM10 column
-             ducmass(i,j)=pp10cb(i,j)
-             !PM25 column
-             ducmass25(i,j)=pp25cb(i,j)
-           endif !nasa_on
+                if(trim(fieldname)=='pm10') then
+                 !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,pp10,arrayr43d,fillvalue,spval)
+                  do l=1,lm
+                    do j=jsta,jend
+                      do i=ista, iend
+                        pp10(i,j,l,1) = max(arrayr43d(i,j,l),0.0)
+                        if(abs(arrayr43d(i,j,l)-fillvalue) < small) pp10(i,j,l,1) = spval
+                      enddo
+                    enddo
+                  enddo
+                endif
 
-          end do
-        end do
+              endif !nasa_on
 
-        endif !end gocart_on, nasa_on
+            endif !end gocart_on, gccpp_on, nasa_on
 
 
 !3d fields
@@ -4427,7 +4147,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
       enddo file_loop_all
 
 ! recompute full layer of zint
-!$omp parallel do default(none) private(i,j) shared(jsta,jend,im,lp1,spval,zint,fis,ista,iend)
+!$omp parallel do default(none) private(i,j) shared(jsta,jend,lp1,spval,zint,fis,ista,iend)
       do j=jsta,jend
         do i=ista,iend
           if (fis(i,j) /= spval) then
@@ -4441,7 +4161,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
       enddo
 
       do l=lm,1,-1
-!$omp parallel do default(none) private(i,j) shared(l,jsta,jend,im,omga,wh,dpres,zint,spval,ista,iend)
+!$omp parallel do default(none) private(i,j) shared(l,jsta,jend,omga,wh,dpres,zint,spval,ista,iend)
         do j=jsta,jend
           do i=ista,iend
             if(wh(i,j,l) /= spval) then
@@ -4458,7 +4178,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
 !           'lm=',maxval(omga(ista:iend,jsta:jend,lm)),minval(omga(ista:iend,jsta:jend,lm))
 
 ! compute pint from top down
-!$omp parallel do default(none) private(i,j) shared(jsta,jend,im,ak5,pint,ista,iend)
+!$omp parallel do default(none) private(i,j) shared(jsta,jend,ak5,pint,ista,iend)
       do j=jsta,jend
         do i=ista,iend
           pint(i,j,1) = ak5(1)
@@ -4466,7 +4186,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
       end do
 
       do l=2,lp1
-!$omp parallel do default(none) private(i,j) shared(l,jsta,jend,im,pint,dpres,spval,ista,iend)
+!$omp parallel do default(none) private(i,j) shared(l,jsta,jend,pint,dpres,spval,ista,iend)
         do j=jsta,jend
           do i=ista,iend
             if(dpres(i,j,l-1) /= spval) then
@@ -4480,7 +4200,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
 
 !compute pmid from averaged two layer pint
       do l=lm,1,-1
-!$omp parallel do default(none) private(i,j) shared(l,jsta,jend,im,pmid,pint,spval,ista,iend)
+!$omp parallel do default(none) private(i,j) shared(l,jsta,jend,pmid,pint,spval,ista,iend)
         do j=jsta,jend
           do i=ista,iend
             if(pint(i,j,l+1) /= spval) then
@@ -4492,7 +4212,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
         enddo
       enddo
 
-!$omp parallel do default(none) private(i,j) shared(jsta,jend,im,spval,pt,pd,pint,ista,iend)
+!$omp parallel do default(none) private(i,j) shared(jsta,jend,spval,pt,pd,pint,ista,iend)
       do j=jsta,jend
         do i=ista,iend
           pd(i,j)     = spval
@@ -4503,7 +4223,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
 
 ! compute alpint
       do l=lp1,1,-1
-!$omp parallel do default(none) private(i,j) shared(l,jsta,jend,im,alpint,pint,spval,ista,iend)
+!$omp parallel do default(none) private(i,j) shared(l,jsta,jend,alpint,pint,spval,ista,iend)
         do j=jsta,jend
           do i=ista,iend
             if(pint(i,j,l) /= spval) then
@@ -4517,7 +4237,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
 
 ! compute zmid
       do l=lm,1,-1
-!$omp parallel do default(none) private(i,j) shared(l,jsta,jend,im,zmid,zint,pmid,alpint,spval,ista,iend)
+!$omp parallel do default(none) private(i,j) shared(l,jsta,jend,zmid,zint,pmid,alpint,spval,ista,iend)
         do j=jsta,jend
           do i=ista,iend
             if( zint(i,j,l+1)/=spval .and. zint(i,j,l)/=spval .and. pmid(i,j,l) /= spval) then
@@ -4530,13 +4250,6 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
           end do
         end do
       end do
-!        print *,'in post_gfs,zmid=',maxval(zmid(1:im,jsta:jend,1)), &
-!          minval(zmid(1:im,jsta:jend,1)),maxloc(zmid(1:im,jsta:jend,1)), &
-!          'zint=',maxval(zint(1:im,jsta:jend,2)),minval(zint(1:im,jsta:jend,1)),  &
-!          'pmid=',maxval(pmid(1:im,jsta:jend,1)),minval(pmid(1:im,jsta:jend,1)),  &
-!          'alpint=',maxval(alpint(1:im,jsta:jend,2)),minval(alpint(1:im,jsta:jend,2))
-!        print *,'in post_gfs,alpint=',maxval(alpint(1:im,jsta:jend,1)), &
-!          minval(alpint(1:im,jsta:jend,1))
 
 ! surface potential T, and potential T at roughness length
 !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,lp1,sm,ths,sst,thz0,sice,pint)
@@ -4628,51 +4341,220 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
       if(modelname=='FV3R') then
         ! smoke and dust extinction
         !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,zint,taod5503d,aextc55,extsmoke,extdust,spval)
-          do l=1,lm
-            do j=jsta,jend
-              do i=ista, iend
-                if(taod5503d(i,j,l)<spval.and.aextc55(i,j,l)<spval) then
-                  taod5503d(i,j,l)=extsmoke(i,j,l)+extdust(i,j,l)
-                  aextc55(i,j,l)=taod5503d(i,j,l)/(zint(i,j,l)-zint(i,j,l+1))
-                endif
-              enddo
+        do l=1,lm
+          do j=jsta,jend
+            do i=ista, iend
+              if(taod5503d(i,j,l)<spval.and.aextc55(i,j,l)<spval) then
+                taod5503d(i,j,l)=extsmoke(i,j,l)+extdust(i,j,l)
+                aextc55(i,j,l)=taod5503d(i,j,l)/(zint(i,j,l)-zint(i,j,l+1))
+              endif
             enddo
           enddo
+        enddo
 
-             !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,accswe_ice,accswe_land,acsnow)
-             do j=jsta,jend
-               do i=ista, iend
-                 if(accswe_land(i,j)<spval) then
-                   acsnow(i,j) = accswe_land(i,j) 
-                 elseif(accswe_ice(i,j)<spval) then
-                   acsnow(i,j) = accswe_ice(i,j)
-                 else
-                   acsnow(i,j) = spval
-                 endif
-               enddo
-             enddo
+        !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,accswe_ice,accswe_land,acsnow)
+        do j=jsta,jend
+          do i=ista, iend
+            if(accswe_land(i,j)<spval) then
+              acsnow(i,j) = accswe_land(i,j) 
+            elseif(accswe_ice(i,j)<spval) then
+              acsnow(i,j) = accswe_ice(i,j)
+            else
+              acsnow(i,j) = spval
+            endif
+          enddo
+        enddo
+
+        !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,snacc_ice,snacc_land,sndepac)
+        do j=jsta,jend
+          do i=ista, iend
+            if(snacc_land(i,j)<spval) then
+              sndepac(i,j) = snacc_land(i,j) 
+            elseif(snacc_ice(i,j)<spval) then
+              sndepac(i,j) = snacc_ice(i,j)
+            else
+              sndepac(i,j) = spval
+            endif
+          enddo
+        enddo
+
+        deallocate(extsmoke)
+        deallocate(extdust)
+        deallocate(accswe_ice)
+        deallocate(accswe_land)
+        deallocate(snacc_ice)
+        deallocate(snacc_land)
+
+      endif !end FV3R
+
+      ! chmical field computation
+      if(gocart_on .or. gccpp_on .or. nasa_on) then
+        dustcb=0.0
+        dustallcb=0.0
+        do l=1,lm
+          do j=jsta,jend
+            do i=ista, iend
+              if(dust(i,j,l,1)<spval.and.dust(i,j,l,2)<spval.and. &
+                 dust(i,j,l,3)<spval.and.dust(i,j,l,4)<spval.and. &
+                 dust(i,j,l,5)<spval) then
+                dustcb(i,j)=dustcb(i,j)+ &
+                  (dust(i,j,l,1)+0.38*dust(i,j,l,2))* &
+                  dpres(i,j,l)/grav
+                dustallcb(i,j)=dustallcb(i,j)+ &
+                  (dust(i,j,l,1)+dust(i,j,l,2)+ &
+                  dust(i,j,l,3)+0.74*dust(i,j,l,4))* &
+                  dpres(i,j,l)/grav
+              endif
+            enddo
+          enddo
+        enddo
+ 
+        sscb=0.0
+        ssallcb=0.0
+        do l=1,lm
+          do j=jsta,jend
+            do i=ista,iend
+              if(salt(i,j,l,1)<spval.and.salt(i,j,l,2)<spval.and. &
+                 salt(i,j,l,3)<spval.and.salt(i,j,l,4)<spval.and. &
+                 salt(i,j,l,5)<spval) then
+            sscb(i,j)=sscb(i,j)+ &
+         (salt(i,j,l,1)+salt(i,j,l,2)+0.83*salt(i,j,l,3))*  &
+           dpres(i,j,l)/grav
+
+
+          ssallcb(i,j)=ssallcb(i,j)+ &
+         (salt(i,j,l,1)+salt(i,j,l,2)+salt(i,j,l,3)+salt(i,j,l,4))* &
+           dpres(i,j,l)/grav
+              endif
+            enddo
+          enddo
+        end do 
+
+        bccb=0.0
+        occb=0.0
+        do l=1,lm
+          do j=jsta,jend
+            do i=ista,iend
+              if(soot(i,j,l,1)<spval.and.soot(i,j,l,2)<spval)then
+               bccb(i,j)=bccb(i,j)+(soot(i,j,l,1)+soot(i,j,l,2))* &
+               dpres(i,j,l)/grav
+              endif
+
+              if(waso(i,j,l,1)<spval.and.waso(i,j,l,2)<spval)then
+               occb(i,j)=occb(i,j)+ (waso(i,j,l,1)+waso(i,j,l,2))* &
+               dpres(i,j,l)/grav
+              endif
+            enddo
+          enddo
+        end do
 
-            !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,snacc_ice,snacc_land,sndepac)
+        if(nasa_on) then
+          no3cb=0.0
+          nh4cb=0.0
+          do l=1,lm
             do j=jsta,jend
-              do i=ista, iend
-                if(snacc_land(i,j)<spval) then
-                  sndepac(i,j) = snacc_land(i,j) 
-                elseif(snacc_ice(i,j)<spval) then
-                  sndepac(i,j) = snacc_ice(i,j)
-                else
-                  sndepac(i,j) = spval
+              do i=ista,iend
+                if(no3(i,j,l,1)<spval .and. no3(i,j,l,2)<spval .and. &
+                   no3(i,j,l,3)<spval) then
+                   no3cb(i,j)=no3cb(i,j)+ (no3(i,j,l,1)+no3(i,j,l,2)+ &
+                   no3(i,j,l,3) ) * dpres(i,j,l)/grav
+                endif
+                if(nh4(i,j,l,1)<spval)then
+                   nh4cb(i,j)=nh4cb(i,j)+ nh4(i,j,l,1)* &
+                   dpres(i,j,l)/grav
                 endif
               enddo
             enddo
+          end do 
+        endif !end nasa_on
 
-      deallocate(extsmoke)
-      deallocate(extdust)
-      deallocate(accswe_ice)
-      deallocate(accswe_land)
-      deallocate(snacc_ice)
-      deallocate(snacc_land)
+        sulfcb=0.0
+        pp25cb=0.0
+        pp10cb=0.0
+        do l=1,lm
+          do j=jsta,jend
+            do i=ista,iend
+              if(suso(i,j,l,1)<spval)then
+                sulfcb(i,j)=sulfcb(i,j)+ suso(i,j,l,1)* &
+                  dpres(i,j,l)/grav
+              endif
+              if(pp25(i,j,l,1)<spval)then
+                pp25cb(i,j)=pp25cb(i,j)+ pp25(i,j,l,1)* &
+                  dpres(i,j,l)/grav
+              endif
+              if(pp10(i,j,l,1)<spval)then
+                pp10cb(i,j)=pp10cb(i,j)+ pp10(i,j,l,1)* &
+                  dpres(i,j,l)/grav
+              endif
+            enddo
+          enddo
+        enddo ! do loop for l
+ 
+        l=lm
+        do j=jsta,jend
+          do i=ista,iend
 
-      endif !end FV3R
+            tv = t(i,j,l) * (h1+d608*MAX(q(I,J,L),qmin))
+            rhomid(i,j,l) = pmid(i,j,l) / (rd*tv)
+
+            dustcb(i,j) = MAX(dustcb(i,j), 0.0)
+            dustallcb(i,j) = MAX(dustallcb(i,j), 0.0)
+            sscb(i,j) = MAX(sscb(i,j), 0.0)
+            ssallcb(i,j) = MAX(ssallcb(i,j), 0.0)
+            bccb(i,j) = MAX(bccb(i,j), 0.0)
+            occb(i,j) = MAX(occb(i,j), 0.0)
+            sulfcb(i,j) = MAX(sulfcb(i,j), 0.0)
+            if(nasa_on) then
+              no3cb(i,j) = MAX(no3cb(i,j), 0.0)
+              nh4cb(i,j) = MAX(nh4cb(i,j), 0.0)
+            endif
+            pp25cb(i,j) = MAX(pp25cb(i,j), 0.0)
+            pp10cb(i,j) = MAX(pp10cb(i,j), 0.0)
+
+           ! Surface PM25 dust and seasalt      
+           dustpm(i,j)=(dust(i,j,l,1)+0.38*dust(i,j,l,2))*rhomid(i,j,l) !ug/m3
+           dustpm10(i,j)=(dust(i,j,l,1)+dust(i,j,l,2)+dust(i,j,l,3)+ &
+             0.74*dust(i,j,l,4))*rhomid(i,j,l) !ug/m3
+           sspm(i,j)=(salt(i,j,l,1)+salt(i,j,l,2)+ &
+             0.83*salt(i,j,l,3))*rhomid(i,j,l)  !ug/m3 
+            
+           if(gocart_on .or. gccpp_on) then
+
+             !Surface PM10 concentration 
+             dusmass(i,j)=(dust(i,j,l,1)+dust(i,j,l,2)+dust(i,j,l,3)+ &
+               0.74*dust(i,j,l,4)+salt(i,j,l,1)+salt(i,j,l,2)+salt(i,j,l,3)+ &
+               salt(i,j,l,4) + soot(i,j,l,1)+soot(i,j,l,2)+waso(i,j,l,1)+ &
+               waso(i,j,l,2) +suso(i,j,l,1)+pp25(i,j,l,1)+pp10(i,j,l,1)) &
+               *rhomid(i,j,l)  !ug/m3
+             !Surface PM25 concentration       
+             dusmass25(i,j)=(dust(i,j,l,1)+0.38*dust(i,j,l,2)+ &
+               salt(i,j,l,1)+salt(i,j,l,2)+0.83*salt(i,j,l,3) + &
+               soot(i,j,l,1)+soot(i,j,l,2)+waso(i,j,l,1)+ &
+               waso(i,j,l,2) +suso(i,j,l,1)+pp25(i,j,l,1))*rhomid(i,j,l)  !ug/m3
+
+             !PM10 column
+             ducmass(i,j)=dustallcb(i,j)+ssallcb(i,j)+bccb(i,j)+ &
+               occb(i,j)+sulfcb(i,j)+pp25cb(i,j)+pp10cb(i,j)
+             !PM25 column
+             ducmass25(i,j)=dustcb(i,j)+sscb(i,j)+bccb(i,j)+occb(i,j) &
+               +sulfcb(i,j)+pp25cb(i,j)
+
+           elseif(nasa_on) then
+             !Surface PM10 concentration
+             dusmass(i,j)=pp10(i,j,l,1)*rhomid(i,j,l)  !ug/m3
+             !Surface PM25 concentration       
+             dusmass25(i,j)=pp25(i,j,l,1)*rhomid(i,j,l)  !ug/m3
+
+             !PM10 column
+             ducmass(i,j)=pp10cb(i,j)
+             !PM25 column
+             ducmass25(i,j)=pp25cb(i,j)
+
+           endif !nasa_on
+         end do
+       end do
+
+      endif
 
 ! generate look up table for lifted parcel calculations
       thl    = 210.
diff --git a/upp b/upp
index dccb32176..7529528c9 160000
--- a/upp
+++ b/upp
@@ -1 +1 @@
-Subproject commit dccb32176930676ef2a258eb65571ab4e3f7e7a4
+Subproject commit 7529528c9a5bc99188349e7d2d418ec68eaedbce

From 9430e2402e78e44d567f66b39c2b71b5615ea223 Mon Sep 17 00:00:00 2001
From: Wen Meng <wen.meng@noaa.gov>
Date: Tue, 13 Jun 2023 01:34:13 +0000
Subject: [PATCH 2/9] sync changes from offline post interface

---
 io/post_fv3.F90 | 136 ++++++++++++++++++------------------------------
 1 file changed, 50 insertions(+), 86 deletions(-)

diff --git a/io/post_fv3.F90 b/io/post_fv3.F90
index c8c6fa5c3..eff4f3cd5 100644
--- a/io/post_fv3.F90
+++ b/io/post_fv3.F90
@@ -544,12 +544,13 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
                              sfcvxi, t10m, t10avg, psfcavg, akhsavg, akmsavg,  &
                              albedo, tg, prate_max, pwat, snow_acm, snow_bkt,  &
                              acgraup, graup_bucket, acfrain, frzrn_bucket,     &
-                             ltg1_max, ltg2_max, ltg3_max, aodtot, ebb, hwp,   &
+                             ltg1_max, ltg2_max, ltg3_max, ebb, hwp,           &
                              aod550,du_aod550,ss_aod550,su_aod550,oc_aod550,   &
                              bc_aod550,maod,                                   &
                              dustpm10, dustcb, bccb, occb, sulfcb, sscb,       &
                              dustallcb, ssallcb, dustpm, sspm, pp25cb, pp10cb, &
-                             no3cb, nh4cb, dusmass, ducmass, dusmass25,ducmass25
+                             no3cb, nh4cb, dusmass, ducmass, dusmass25,ducmass25, &
+                             snownc, graupelnc, qrmax
       use soil,        only: sldpth, sh2o, smc, stc, sllevel
       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, &
@@ -605,9 +606,8 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
       real,dimension(:),    allocatable :: slat,qstl
       real,external::FPVSNEW
       real,dimension(:,:),allocatable :: dummy, p2d, t2d, q2d,  qs2d,  &
-                             cw2d, cfr2d, accswe_ice, accswe_land, &
-                             snacc_land, snacc_ice
-      real,dimension(:,:,:),allocatable :: extsmoke, extdust
+                             cw2d, cfr2d, snacc_land, snacc_ice
+      real,dimension(:,:,:),allocatable :: ext550
       character(len=80)              :: fieldname, wrtFBName, flatlon, &
                                         VarName
       type(ESMF_Grid)                :: wrtGrid
@@ -688,10 +688,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
 
       !Allocate for regional models only
       if(modelname=='FV3R') then
-        allocate(extsmoke(ista:iend,jsta:jend,lm))
-        allocate(extdust(ista:iend,jsta:jend,lm))
-        allocate(accswe_ice(ista:iend,jsta:jend))
-        allocate(accswe_land(ista:iend,jsta:jend))
+        allocate(ext550(ista:iend,jsta:jend,lm))
         allocate(snacc_ice(ista:iend,jsta:jend))
         allocate(snacc_land(ista:iend,jsta:jend))
       endif
@@ -999,17 +996,6 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
               enddo
             endif
 
-            ! total aod
-            if(trim(fieldname)=='aodtot') then
-              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,aodtot,arrayr42d,fillValue,spval)
-              do j=jsta,jend
-                do i=ista, iend
-                  aodtot(i,j)=arrayr42d(i,j)
-                  if(abs(arrayr42d(i,j)-fillValue) < small) aodtot(i,j)=spval
-                enddo
-              enddo
-            endif
-
             ! biomass burning emissions
             if(trim(fieldname)=='ebb_smoke_hr') then
               !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,ebb,arrayr42d,fillValue,spval)
@@ -1252,6 +1238,28 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
               enddo
             endif
 
+            !time step snow (in m)
+            if(trim(fieldname)=='snow') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,snownc,arrayr42d,sm,fillValue)
+              do j=jsta,jend
+                do i=ista, iend
+                  snownc(i,j) = arrayr42d(i,j)
+                  if (abs(arrayr42d(i,j)-fillValue) < small) snownc(i,j) = spval
+                enddo
+              enddo
+            endif
+
+            !time step graupel (in m)
+            if(trim(fieldname)=='graupel') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,graupelnc,arrayr42d,sm,fillValue)
+              do j=jsta,jend
+                do i=ista, iend
+                  graupelnc(i,j) = arrayr42d(i,j)
+                  if (abs(arrayr42d(i,j)-fillValue) < small) graupelnc(i,j) = spval 
+                enddo
+              enddo
+            endif
+
             ! max hourly surface precipitation rate
             if(trim(fieldname)=='pratemax') then
               !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,prate_max,arrayr42d,sm,fillValue)
@@ -2235,25 +2243,6 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
             endif
 
             if(modelname=='FV3R')then
-            !acsnow
-            if(trim(fieldname)=='accswe_land') then
-              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,accswe_land,arrayr42d,fillvalue,spval)
-              do j=jsta,jend
-                do i=ista, iend
-                  accswe_land(i,j) = arrayr42d(i,j)
-                  if(abs(arrayr42d(i,j)-fillvalue)<small) accswe_land(i,j) = spval
-                enddo
-              enddo
-            endif
-            if(trim(fieldname)=='accswe_ice') then
-              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,accswe_ice,arrayr42d,fillvalue,spval)
-              do j=jsta,jend
-                do i=ista, iend
-                  accswe_ice(i,j) = arrayr42d(i,j)
-                  if(abs(arrayr42d(i,j)-fillvalue)<small) accswe_ice(i,j) = spval
-                enddo
-              enddo
-            endif
 
             !sndepac
             if(trim(fieldname)=='snacc_land') then
@@ -3729,6 +3718,8 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
               endif
             endif
 
+            if(modelname=='FV3R') then
+
             ! model level smoke
             if(trim(fieldname)=='smoke') then
               !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,smoke,arrayr43d,spval,fillvalue)
@@ -3755,36 +3746,19 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
               enddo
             endif
 
-            if(modelname=='FV3R') then
-
-            ! model level smoke_ext
-            if(trim(fieldname)=='smoke_ext') then
-              !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,extsmoke,arrayr43d,spval,fillvalue)
-              do l=1,lm
-                do j=jsta,jend
-                  do i=ista, iend
-                    extsmoke(i,j,l)=arrayr43d(i,j,l)
-                    if(abs(arrayr43d(i,j,l)-fillvalue)<small) extsmoke(i,j,l) = spval
-                  enddo
-                enddo
-              enddo
-            endif
-
-            ! model level dust_ext
-            if(trim(fieldname)=='dust_ext') then
-              !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,extdust,arrayr43d,spval,fillvalue)
+            ! model level ext550 extinction
+            if(trim(fieldname)=='ext550') then
+              !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,ext550,arrayr43d,spval,fillvalue)
               do l=1,lm
                 do j=jsta,jend
                   do i=ista, iend
-                    extdust(i,j,l)=arrayr43d(i,j,l)
-                    if(abs(arrayr43d(i,j,l)-fillvalue)<small) extdust(i,j,l) = spval
+                    ext550(i,j,l)=arrayr43d(i,j,l)
+                    if(abs(arrayr43d(i,j,l)-fillvalue)<small) ext550(i,j,l) = spval
                   enddo
                 enddo
               enddo
             endif
 
-            endif !end FV3R
-
             ! model level coarse dust
             if(trim(fieldname)=='coarsepm') then
               !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,coarsepm,arrayr43d,spval,fillvalue)
@@ -3798,6 +3772,8 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
               enddo
             endif
 
+            endif !end FV3R
+
             ! Thompson scheme cloud ice effective radius
             if(trim(fieldname)=='cieffr') then
               !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,effri,arrayr43d,spval,fillvalue)
@@ -4276,21 +4252,25 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
 !      print *,'in post_gfs,ths=',maxval(ths(1:im,jsta:jend)), &
 !          minval(ths(1:im,jsta:jend))
 
-! compute cwm for gfdlmp
-!      if(  imp_physics == 11 ) then
+! compute cwm and max qrain in the column to be used later in precip type computation 
+        do j=jsta,jend
+          do i=ista,iend
+            qrmax(i,j)=0.
+          enddo
+        enddo
         do l=1,lm
-!$omp parallel do default(none) private(i,j) shared(l,jsta,jend,ista,iend,cwm,qqg,qqs,qqr,qqi,qqw,spval)
+!$omp parallel do default(none) private(i,j) shared(l,jsta,jend,ista,iend,cwm,qrmax,qqg,qqs,qqr,qqi,qqw,spval)
           do j=jsta,jend
             do i=ista,iend
-              if( qqg(i,j,l) /= spval) then
+              if( qqr(i,j,l) /= spval) then
                 cwm(i,j,l) = qqg(i,j,l)+qqs(i,j,l)+qqr(i,j,l)+qqi(i,j,l)+qqw(i,j,l)
+                qrmax(i,j)=max(qrmax(i,j),qqr(i,j,l))
               else
                 cwm(i,j,l) = spval
               endif
             enddo
           enddo
         enddo
-!      endif
 
 ! estimate 2m pres and convert t2m to theta
 !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,lm,pshltr,pint,tshltr,spval)
@@ -4340,31 +4320,18 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
 
       if(modelname=='FV3R') then
         ! smoke and dust extinction
-        !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,zint,taod5503d,aextc55,extsmoke,extdust,spval)
+        !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,zint,taod5503d,aextc55,ext550,spval)
         do l=1,lm
           do j=jsta,jend
             do i=ista, iend
-              if(taod5503d(i,j,l)<spval.and.aextc55(i,j,l)<spval) then
-                taod5503d(i,j,l)=extsmoke(i,j,l)+extdust(i,j,l)
+              if(ext550(i,j,l)<spval) then
+                taod5503d(i,j,l)=ext550(i,j,l)
                 aextc55(i,j,l)=taod5503d(i,j,l)/(zint(i,j,l)-zint(i,j,l+1))
               endif
             enddo
           enddo
         enddo
 
-        !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,accswe_ice,accswe_land,acsnow)
-        do j=jsta,jend
-          do i=ista, iend
-            if(accswe_land(i,j)<spval) then
-              acsnow(i,j) = accswe_land(i,j) 
-            elseif(accswe_ice(i,j)<spval) then
-              acsnow(i,j) = accswe_ice(i,j)
-            else
-              acsnow(i,j) = spval
-            endif
-          enddo
-        enddo
-
         !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,snacc_ice,snacc_land,sndepac)
         do j=jsta,jend
           do i=ista, iend
@@ -4378,10 +4345,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
           enddo
         enddo
 
-        deallocate(extsmoke)
-        deallocate(extdust)
-        deallocate(accswe_ice)
-        deallocate(accswe_land)
+        deallocate(ext550)
         deallocate(snacc_ice)
         deallocate(snacc_land)
 

From f66eaaf59fadc77710ecdeaf1e026f269214bb00 Mon Sep 17 00:00:00 2001
From: Wen Meng <wen.meng@noaa.gov>
Date: Wed, 14 Jun 2023 16:14:32 +0000
Subject: [PATCH 3/9] More clean up

---
 io/post_fv3.F90 | 22 ++--------------------
 1 file changed, 2 insertions(+), 20 deletions(-)

diff --git a/io/post_fv3.F90 b/io/post_fv3.F90
index eff4f3cd5..82c62f05a 100644
--- a/io/post_fv3.F90
+++ b/io/post_fv3.F90
@@ -693,15 +693,6 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
         allocate(snacc_land(ista:iend,jsta:jend))
       endif
 
-      !wm
-!! GFS does not have convective precip
-!!$omp parallel do default(none) private(i,j) shared(jsta,jend,cprate,ista,iend)
-!      do j=jsta,jend
-!        do i=ista,iend
-!          cprate(i,j) = 0.
-!        enddo
-!      enddo
-
 !
 ! GFS doesn not yet output soil layer thickness, assign SLDPTH to be the same as nam
       sldpth(1) = 0.10
@@ -709,18 +700,9 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
       sldpth(3) = 0.6
       sldpth(4) = 1.0
 
-! GFS does not output time averaged convective and strat cloud fraction, set acfrcv to spval, n
-! cfrcv to 1
-!                     time averaged cloud fraction, set acfrst to spval, ncfrst to 1
-!                     UNDERGROUND RUNOFF, bgroff
-!                     inst incoming sfc longwave
-!                     inst incoming sfc shortwave, rswin
-!                     inst incoming clear sky sfc shortwave, rswinc
-!                     inst outgoing sfc shortwave, rswout
-!                     snow phase change heat flux, snopcx
-! GFS does not use total momentum flux,sfcuvx
+! set ncfrcv to 1, ncfrst to 1
 !$omp parallel do default(none),private(i,j),shared(jsta,jend,spval,ista,iend), &
-!$omp& shared(acfrcv,ncfrcv,acfrst,ncfrst,bgroff,rswinc)
+!$omp& shared(ncfrcv,ncfrst)
       do j=jsta,jend
         do i=ista,iend
           ncfrcv(i,j) = 1.0

From a49d3d98544cfb245a91c896635512a5bd5f0b18 Mon Sep 17 00:00:00 2001
From: Wen Meng <wen.meng@noaa.gov>
Date: Mon, 26 Jun 2023 19:29:51 +0000
Subject: [PATCH 4/9] Update upp revision to 0544d07

---
 upp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/upp b/upp
index 7529528c9..0544d0741 160000
--- a/upp
+++ b/upp
@@ -1 +1 @@
-Subproject commit 7529528c9a5bc99188349e7d2d418ec68eaedbce
+Subproject commit 0544d074185a3fb67503ae918ef3ade2eab29a2f

From 98c3f4d921b5c973fb6e202ded71ba5aeab0fe9f Mon Sep 17 00:00:00 2001
From: Wen Meng <wen.meng@noaa.gov>
Date: Fri, 30 Jun 2023 19:08:41 +0000
Subject: [PATCH 5/9] Add initialization for some local variables.

---
 io/post_fv3.F90 | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/io/post_fv3.F90 b/io/post_fv3.F90
index 82c62f05a..c86c27630 100644
--- a/io/post_fv3.F90
+++ b/io/post_fv3.F90
@@ -691,6 +691,21 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
         allocate(ext550(ista:iend,jsta:jend,lm))
         allocate(snacc_ice(ista:iend,jsta:jend))
         allocate(snacc_land(ista:iend,jsta:jend))
+
+        do j=jsta,jend
+          do i=ista,iend
+            snacc_ice(i,j)=spval
+            snacc_land(i,j)=spval
+          end do
+        end do
+
+        do l=1,lm
+          do j=jsta,jend
+            do i=ista,iend
+              ext550(i,j,l)=spval
+            end do
+          end do
+        end do
       endif
 
 !
@@ -4309,6 +4324,9 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
               if(ext550(i,j,l)<spval) then
                 taod5503d(i,j,l)=ext550(i,j,l)
                 aextc55(i,j,l)=taod5503d(i,j,l)/(zint(i,j,l)-zint(i,j,l+1))
+              else
+                taod5503d(i,j,l)=spval
+                aextc55(i,j,l)=spval
               endif
             enddo
           enddo

From 48d6fbac0930e422c1d7d8d5b6e193002a5a894f Mon Sep 17 00:00:00 2001
From: Wen Meng <wen.meng@noaa.gov>
Date: Sat, 1 Jul 2023 01:35:30 +0000
Subject: [PATCH 6/9] Change undefined value of aextc55 as 0

---
 io/post_fv3.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/io/post_fv3.F90 b/io/post_fv3.F90
index c86c27630..ff52a330d 100644
--- a/io/post_fv3.F90
+++ b/io/post_fv3.F90
@@ -4326,7 +4326,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
                 aextc55(i,j,l)=taod5503d(i,j,l)/(zint(i,j,l)-zint(i,j,l+1))
               else
                 taod5503d(i,j,l)=spval
-                aextc55(i,j,l)=spval
+                aextc55(i,j,l)=0.
               endif
             enddo
           enddo

From 827bfcade4a926b4110ff7c67ab671005c002ad0 Mon Sep 17 00:00:00 2001
From: Wen Meng <wen.meng@noaa.gov>
Date: Thu, 6 Jul 2023 20:43:05 +0000
Subject: [PATCH 7/9] Read soill to level 9 for RUC LSM

---
 io/post_fv3.F90 | 63 +++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 63 insertions(+)

diff --git a/io/post_fv3.F90 b/io/post_fv3.F90
index ff52a330d..2b2b70564 100644
--- a/io/post_fv3.F90
+++ b/io/post_fv3.F90
@@ -1626,6 +1626,69 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
               enddo
             endif
 
+            if(nsoil==9) then
+            ! liquid volumetric soil mpisture in fraction
+            if(trim(fieldname)=='soill5') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,sh2o,arrayr42d,sm,fillValue)
+              do j=jsta,jend
+                do i=ista, iend
+                  sh2o(i,j,5) = arrayr42d(i,j)
+                  if( abs(arrayr42d(i,j)-fillValue) < small) sh2o(i,j,5) = spval
+                  if (sm(i,j) /= 0.0) sh2o(i,j,5) = spval
+                enddo
+              enddo
+            endif
+
+            ! liquid volumetric soil mpisture in fraction
+            if(trim(fieldname)=='soill6') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,sh2o,arrayr42d,sm,fillValue)
+              do j=jsta,jend
+                do i=ista, iend
+                  sh2o(i,j,6) = arrayr42d(i,j)
+                  if( abs(arrayr42d(i,j)-fillValue) < small) sh2o(i,j,6) = spval
+                  if (sm(i,j) /= 0.0) sh2o(i,j,6) = spval
+                enddo
+              enddo
+            endif
+
+            ! liquid volumetric soil mpisture in fraction
+            if(trim(fieldname)=='soill7') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,sh2o,arrayr42d,sm,fillValue)
+              do j=jsta,jend
+                do i=ista, iend
+                  sh2o(i,j,7) = arrayr42d(i,j)
+                  if( abs(arrayr42d(i,j)-fillValue) < small) sh2o(i,j,7) = spval
+                  if (sm(i,j) /= 0.0) sh2o(i,j,7) = spval
+                enddo
+              enddo
+            endif
+
+            ! liquid volumetric soil mpisture in fraction
+            if(trim(fieldname)=='soill8') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,sh2o,arrayr42d,sm,fillValue)
+              do j=jsta,jend
+                do i=ista, iend
+                  sh2o(i,j,8) = arrayr42d(i,j)
+                  if( abs(arrayr42d(i,j)-fillValue) < small) sh2o(i,j,8) = spval
+                  if (sm(i,j) /= 0.0) sh2o(i,j,8) = spval
+                enddo
+              enddo
+            endif
+
+            ! liquid volumetric soil mpisture in fraction
+            if(trim(fieldname)=='soill9') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,sh2o,arrayr42d,sm,fillValue)
+              do j=jsta,jend
+                do i=ista, iend
+                  sh2o(i,j,9) = arrayr42d(i,j)
+                  if( abs(arrayr42d(i,j)-fillValue) < small) sh2o(i,j,9) = spval
+                  if (sm(i,j) /= 0.0) sh2o(i,j,9) = spval
+                enddo
+              enddo
+            endif
+
+            endif !nsoil
+
             ! volumetric soil moisture
             if(trim(fieldname)=='soilw1') then
               !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,smc,arrayr42d,sm,fillValue)

From 60b546fd13a9c005b5f3e19157b49af88cc2576f Mon Sep 17 00:00:00 2001
From: Wen Meng <wen.meng@noaa.gov>
Date: Mon, 17 Jul 2023 01:26:01 +0000
Subject: [PATCH 8/9] Correct tke reading and update upp revision.

---
 io/post_fv3.F90 | 8 ++++++--
 upp             | 2 +-
 2 files changed, 7 insertions(+), 3 deletions(-)

diff --git a/io/post_fv3.F90 b/io/post_fv3.F90
index 2b2b70564..f8dcc80de 100644
--- a/io/post_fv3.F90
+++ b/io/post_fv3.F90
@@ -3726,13 +3726,17 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
             endif
 
             ! model level tke
-            if(trim(fieldname)=='tke') then
+            if(trim(fieldname)=='qke') then
               !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,q2,arrayr43d, fillvalue,spval)
               do l=1,lm
                 do j=jsta,jend
                   do i=ista, iend
                     q2(i,j,l)=arrayr43d(i,j,l)
-                    if(abs(arrayr43d(i,j,l)-fillvalue)<small) q2(i,j,l) = spval
+                    if(abs(arrayr43d(i,j,l)-fillvalue)<small) then
+                      q2(i,j,l) = spval
+                    else
+                      q2(i,j,l) = q2(i,j,l)/2.0
+                    endif
                   enddo
                 enddo
 !              print *,'in gfs_post, get tke=',maxval(q2(:,:,l)), minval(q2(:,:,l)),'l=',l
diff --git a/upp b/upp
index 0544d0741..83e59b5d9 160000
--- a/upp
+++ b/upp
@@ -1 +1 @@
-Subproject commit 0544d074185a3fb67503ae918ef3ade2eab29a2f
+Subproject commit 83e59b5d92ec678bcab11dd59919fac189a11a4c

From 898c2e45f56d538eda0168a00c1f33bfa129d784 Mon Sep 17 00:00:00 2001
From: Wen Meng <wen.meng@noaa.gov>
Date: Mon, 17 Jul 2023 20:08:15 +0000
Subject: [PATCH 9/9] Update upp revision to baa7751

---
 upp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/upp b/upp
index 83e59b5d9..baa77515f 160000
--- a/upp
+++ b/upp
@@ -1 +1 @@
-Subproject commit 83e59b5d92ec678bcab11dd59919fac189a11a4c
+Subproject commit baa77515f182eb6d6f2f8b57bd1e0ffb742761e8