From 0451a89ee7bb3e64a544232f024521a5e24f73be Mon Sep 17 00:00:00 2001
From: Wen Meng <wen.meng@noaa.gov>
Date: Thu, 30 Mar 2023 17:03:11 +0000
Subject: [PATCH 01/11] update upp revision to ce35e5f

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

diff --git a/upp b/upp
index 22cfb88de..ce35e5f8d 160000
--- a/upp
+++ b/upp
@@ -1 +1 @@
-Subproject commit 22cfb88dee0afa74bb753bdf85596b4ceff5e913
+Subproject commit ce35e5f8d2f44a95b71b6d142cd854f73188d074

From 5d75ad0fa220689f3237206b1e1b651591beba7e Mon Sep 17 00:00:00 2001
From: Wen Meng <wen.meng@noaa.gov>
Date: Thu, 30 Mar 2023 17:08:32 +0000
Subject: [PATCH 02/11] Update read interface for RRFS

---
 io/post_fv3.F90 | 407 ++++++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 398 insertions(+), 9 deletions(-)

diff --git a/io/post_fv3.F90 b/io/post_fv3.F90
index 33f098ab9..56555c006 100644
--- a/io/post_fv3.F90
+++ b/io/post_fv3.F90
@@ -86,7 +86,8 @@ subroutine post_run_fv3(wrt_int_state,grid_id,mype,mpicomp,lead_write, &
 
       grib      = "grib2"
       gridtype  = "A"
-      nsoil     = 4
+      !nsoil     = 4
+      nsoil     = wrt_int_state%nsoil
       nwtpg     = wrt_int_state%petcount
       jts       = wrt_int_state%out_grid_info(grid_id)%j_start     !<-- Starting J of this write task's subsection
       jte       = wrt_int_state%out_grid_info(grid_id)%j_end       !<-- Ending J of this write task's subsection
@@ -493,17 +494,19 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
 !     Jul 2019    J. Wang      Initial code
 !     Apr 2022    W. Meng      Unify set_postvars_gfs and
 !                               set_postvars_regional to set_postvars_fv3
+!     Apr 2023    W. Meng      Sync RRFS changes from off-line post
 !
 !-----------------------------------------------------------------------
 !*** set up post fields from nmint_state
 !-----------------------------------------------------------------------
 !
       use esmf
+      use vrbls4d,     only: dust, smoke, fv3dust, coarsepm
       use vrbls3d,     only: t, q, uh, vh, wh, alpint, dpres, zint, zmid, o3,  &
                              qqr, qqs, cwm, qqi, qqw, qqg, omga, cfr, pmid,    &
                              q2, rlwtt, rswtt, tcucn, tcucns, train, el_pbl,   &
                              pint, exch_h, ref_10cm, qqni, qqnr, qqnwfa,       &
-                             qqnifa
+                             qqnifa, effri, effrl, effrs, aextc55, taod5503d
       use vrbls2d,     only: f, pd, sigt4, fis, pblh, ustar, z0, ths, qs, twbs,&
                              qwbs, avgcprate, cprate, avgprec, prec, lspa, sno,&
                              cldefi, th10, q10, tshltr, pshltr, albase,        &
@@ -520,7 +523,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
                              acsnow, acsnom, sst, thz0, qz0, uz0, vz0, ptop,   &
                              htop, pbot, hbot, ptopl, pbotl, ttopl, ptopm,     &
                              pbotm, ttopm, ptoph, pboth, pblcfr, ttoph, runoff,&
-                             tecan, tetran, tedir, twa,                        &
+                             tecan, tetran, tedir, twa, sndepac,               &
                              maxtshltr, mintshltr, maxrhshltr, minrhshltr,     &
                              dzice, smcwlt, suntime, fieldcapa, htopd, hbotd,  &
                              htops, hbots, aswintoa, maxqshltr, minqshltr,     &
@@ -536,8 +539,8 @@ 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
-      use soil,        only: sldpth, sh2o, smc, stc
+                             ltg1_max, ltg2_max, ltg3_max, aodtot, ebb, hwp
+      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, &
                              ista, iend, ista_2l, iend_2u, ista_m,iend_m, &
@@ -545,7 +548,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
                              tprec, tclod, trdlw, trdsw, tsrfc, tmaxmin, theat, &
                              ardlw, ardsw, asrfc, avrain, avcnvc, iSF_SURFACE_PHYSICS,&
                              td3d, idat, sdat, ifhr, ifmin, dt, nphs, dtq2, pt_tbl, &
-                             alsl, spl, ihrst, modelname
+                             alsl, spl, ihrst, modelname, nsoil
       use params_mod,  only: erad, dtr, capa, p1000, small
       use gridspec_mod,only: latstart, latlast, lonstart, lonlast, cenlon, cenlat, &
                              dxval, dyval, truelat2, truelat1, psmapf, cenlat,     &
@@ -590,7 +593,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
+                             cw2d, cfr2d, buf, buf2
+      real,dimension(:,:,:),allocatable :: extsmoke, extdust
       character(len=80)              :: fieldname, wrtFBName, flatlon
       type(ESMF_Grid)                :: wrtGrid
       type(ESMF_Field)               :: theField
@@ -679,6 +683,12 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
       dtq2 = wrt_int_state%dtp
       nphs = 2.
       dt   = dtq2/nphs
+
+      allocate(extsmoke(ista:iend,jsta:jend,lm))
+      allocate(extdust(ista:iend,jsta:jend,lm))
+      allocate(buf(ista:iend,jsta:jend))
+      allocate(buf2(ista:iend,jsta:jend))
+
 !
 ! GFS does not have convective cloud efficiency
 !                   similated precip
@@ -799,9 +809,9 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
 !$omp& shared(smstav,sfcevp,acsnow,acsnom,qz0,uz0,vz0,maxrhshltr,minrhshltr)
       do j=jsta_2l,jend_2u
         do i=ista_2l,iend_2u
-          smstav(i,j) = spval
+          !smstav(i,j) = spval
           sfcevp(i,j) = spval
-          acsnow(i,j) = spval
+          !acsnow(i,j) = spval
           acsnom(i,j) = spval
           qz0(i,j)    = spval
           uz0(i,j)    = spval
@@ -1089,6 +1099,39 @@ 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)
+              do j=jsta,jend
+                do i=ista, iend
+                  ebb(i,j)=arrayr42d(i,j)
+                  if(abs(arrayr42d(i,j)-fillValue) < small) ebb(i,j)=spval
+                enddo
+              enddo
+            endif
+
+            ! wildfire potential
+            if(trim(fieldname)=='hwp') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,hwp,arrayr42d,fillValue,spval)
+              do j=jsta,jend
+                do i=ista, iend
+                  hwp(i,j)=arrayr42d(i,j)
+                  if(abs(arrayr42d(i,j)-fillValue) < small) hwp(i,j)=spval
+                enddo
+              enddo
+            endif
+
             ! frictional velocity
             if(trim(fieldname)=='fricv') then
               !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,ustar,arrayr42d,fillValue,spval)
@@ -1616,6 +1659,20 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
               enddo
             endif
 
+            !assign soil depths for RUC LSM, hard wire 9 soil depths here
+            !so they aren't missing.
+            if (nsoil==9) then
+              sllevel(1) = 0.0
+              sllevel(2) = 0.01
+              sllevel(3) = 0.04
+              sllevel(4) = 0.1
+              sllevel(5) = 0.3
+              sllevel(6) = 0.6
+              sllevel(7) = 1.0
+              sllevel(8) = 1.6
+              sllevel(9) = 3.0
+            endif 
+
             ! liquid volumetric soil mpisture in fraction
             if(trim(fieldname)=='soill1') then
               !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,sh2o,arrayr42d,sm,fillValue)
@@ -1712,6 +1769,65 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
               enddo
             endif
 
+            if(nsoil==9) then
+            ! volumetric soil moisture
+            if(trim(fieldname)=='soilw5') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,smc,arrayr42d,sm,fillValue)
+              do j=jsta,jend
+                do i=ista, iend
+                  smc(i,j,5) = arrayr42d(i,j)
+                  if( abs(arrayr42d(i,j)-fillValue) < small) smc(i,j,5) = spval
+                  if (sm(i,j) /= 0.0) smc(i,j,5) = spval
+                enddo
+              enddo
+            endif
+
+            if(trim(fieldname)=='soilw6') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,smc,arrayr42d,sm,fillValue)
+              do j=jsta,jend
+                do i=ista, iend
+                  smc(i,j,6) = arrayr42d(i,j)
+                  if( abs(arrayr42d(i,j)-fillValue) < small) smc(i,j,6) = spval
+                  if (sm(i,j) /= 0.0) smc(i,j,6) = spval
+                enddo
+              enddo
+            endif
+
+            if(trim(fieldname)=='soilw7') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,smc,arrayr42d,sm,fillValue)
+              do j=jsta,jend
+                do i=ista, iend
+                  smc(i,j,7) = arrayr42d(i,j)
+                  if( abs(arrayr42d(i,j)-fillValue) < small) smc(i,j,7) = spval
+                  if (sm(i,j) /= 0.0) smc(i,j,7) = spval
+                enddo
+              enddo
+            endif
+
+            if(trim(fieldname)=='soilw8') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,smc,arrayr42d,sm,fillValue)
+              do j=jsta,jend
+                do i=ista, iend
+                  smc(i,j,8) = arrayr42d(i,j)
+                  if( abs(arrayr42d(i,j)-fillValue) < small) smc(i,j,8) = spval
+                  if (sm(i,j) /= 0.0) smc(i,j,8) = spval
+                enddo
+              enddo
+            endif
+
+            if(trim(fieldname)=='soilw9') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,smc,arrayr42d,sm,fillValue)
+              do j=jsta,jend
+                do i=ista, iend
+                  smc(i,j,9) = arrayr42d(i,j)
+                  if( abs(arrayr42d(i,j)-fillValue) < small) smc(i,j,9) = spval
+                  if (sm(i,j) /= 0.0) smc(i,j,9) = spval
+                enddo
+              enddo
+            endif
+
+            endif !nsoil
+
             ! soil temperature
             if(trim(fieldname)=='soilt1') then
               !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,stc,arrayr42d,sm,sice,fillValue)
@@ -1764,6 +1880,75 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
               enddo
             endif
 
+            if(nsoil==9) then
+
+            ! soil temperature
+            if(trim(fieldname)=='soilt5') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,stc,arrayr42d,sm,sice,fillValue)
+              do j=jsta,jend
+                do i=ista, iend
+                  stc(i,j,5) = arrayr42d(i,j)
+                  if( abs(arrayr42d(i,j)-fillValue) < small) stc(i,j,5) = spval
+                  !mask open water areas, combine with sea ice tmp
+                  if (sm(i,j) /= 0.0 .and. sice(i,j) ==0.) stc(i,j,5) = spval
+                enddo
+              enddo
+            endif
+
+            ! soil temperature
+            if(trim(fieldname)=='soilt6') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,stc,arrayr42d,sm,sice,fillValue)
+              do j=jsta,jend
+                do i=ista, iend
+                  stc(i,j,6) = arrayr42d(i,j)
+                  if( abs(arrayr42d(i,j)-fillValue) < small) stc(i,j,6) = spval
+                  !mask open water areas, combine with sea ice tmp
+                  if (sm(i,j) /= 0.0 .and. sice(i,j) ==0.) stc(i,j,6) = spval
+                enddo
+              enddo
+            endif
+
+            ! soil temperature
+            if(trim(fieldname)=='soilt7') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,stc,arrayr42d,sm,sice,fillValue)
+              do j=jsta,jend
+                do i=ista, iend
+                  stc(i,j,7) = arrayr42d(i,j)
+                  if( abs(arrayr42d(i,j)-fillValue) < small) stc(i,j,7) = spval
+                  !mask open water areas, combine with sea ice tmp
+                  if (sm(i,j) /= 0.0 .and. sice(i,j) ==0.) stc(i,j,7) = spval
+                enddo
+              enddo
+            endif
+
+            ! soil temperature
+            if(trim(fieldname)=='soilt8') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,stc,arrayr42d,sm,sice,fillValue)
+              do j=jsta,jend
+                do i=ista, iend
+                  stc(i,j,8) = arrayr42d(i,j)
+                  if( abs(arrayr42d(i,j)-fillValue) < small) stc(i,j,8) = spval
+                  !mask open water areas, combine with sea ice tmp
+                  if (sm(i,j) /= 0.0 .and. sice(i,j) ==0.) stc(i,j,8) = spval
+                enddo
+              enddo
+            endif
+
+            ! soil temperature
+            if(trim(fieldname)=='soilt9') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,stc,arrayr42d,sm,sice,fillValue)
+              do j=jsta,jend
+                do i=ista, iend
+                  stc(i,j,9) = arrayr42d(i,j)
+                  if( abs(arrayr42d(i,j)-fillValue) < small) stc(i,j,9) = spval
+                  !mask open water areas, combine with sea ice tmp
+                  if (sm(i,j) /= 0.0 .and. sice(i,j) ==0.) stc(i,j,9) = spval
+                enddo
+              enddo
+            endif
+
+            endif !nsoil
+
             ! time averaged incoming sfc longwave
             if(trim(fieldname)=='dlwrf_ave') then
               !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,alwin,arrayr42d,fillValue,spval)
@@ -2138,6 +2323,77 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
               enddo
             endif
 
+            ! wetness
+            if(trim(fieldname)=='wetness') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,smstav,arrayr42d,fillvalue,spval)
+              do j=jsta,jend
+                do i=ista, iend
+                  smstav(i,j) = arrayr42d(i,j)
+                  if(abs(arrayr42d(i,j)-fillvalue)<small) smstav(i,j) = spval
+                enddo
+              enddo
+            endif
+
+            !acsnow
+            if(trim(fieldname)=='accswe_land') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,buf,arrayr42d,fillvalue,spval)
+              do j=jsta,jend
+                do i=ista, iend
+                  buf(i,j) = arrayr42d(i,j)
+                  if(abs(arrayr42d(i,j)-fillvalue)<small) buf(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,buf2,arrayr42d,fillvalue,spval)
+              do j=jsta,jend
+                do i=ista, iend
+                  buf2(i,j) = arrayr42d(i,j)
+                  if(abs(arrayr42d(i,j)-fillvalue)<small) buf2(i,j) = spval
+                enddo
+              enddo
+            endif
+            ! !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,buf,buf2,acsnow)
+             do j=jsta,jend
+               do i=ista, iend
+                 if(buf(i,j)<spval .and. buf2(i,j)<spval) then
+                   acsnow(i,j) = buf(i,j) + buf2(i,j)
+                 else
+                   acsnow(i,j) = spval
+                 endif
+               enddo
+             enddo
+
+            !sndepac
+            if(trim(fieldname)=='snacc_land') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,buf,arrayr42d,fillvalue,spval)
+              do j=jsta,jend
+                do i=ista, iend
+                  buf(i,j) = arrayr42d(i,j)
+                  if(abs(arrayr42d(i,j)-fillvalue)<small) buf(i,j) = spval
+                enddo
+              enddo
+            endif
+            if(trim(fieldname)=='snacc_ice') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,buf2,arrayr42d,fillvalue,spval)
+              do j=jsta,jend
+                do i=ista, iend
+                  buf2(i,j) = arrayr42d(i,j)
+                  if(abs(arrayr42d(i,j)-fillvalue)<small) buf2(i,j) = spval
+                enddo
+              enddo
+            endif
+            !!$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,buf,buf2,sndepac)
+            do j=jsta,jend
+              do i=ista, iend
+                if(buf(i,j)<spval .and. buf2(i,j)<spval) then
+                  sndepac(i,j) = buf(i,j) + buf2(i,j)
+                else
+                  sndepac(i,j) = spval
+                endif
+              enddo
+            enddo
+
             ! inst  cloud top pressure
             if(trim(fieldname)=='prescnvclt') then
               !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,ptop,arrayr42d,fillValue)
@@ -3161,6 +3417,18 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
                 enddo
               endif
             else !Other MP
+              ! read cldfra_bl first
+              if(trim(fieldname)=='cldfra_bl') then
+                !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,cfr,arrayr43d,fillvalue,spval)
+                do l=1,lm
+                  do j=jsta,jend
+                    do i=ista, iend
+                      cfr(i,j,l) = arrayr43d(i,j,l)
+                      if(abs(arrayr43d(i,j,l)-fillvalue)<small) cfr(i,j,l) = spval
+                    enddo
+                  enddo
+                enddo
+              endif
               if(trim(fieldname)=='cldfra') then
                 !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,cfr,arrayr43d,fillvalue,spval)
                 do l=1,lm
@@ -3174,6 +3442,110 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
               endif
             endif
 
+            ! 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)
+              do l=1,lm
+                do j=jsta,jend
+                  do i=ista, iend
+                    smoke(i,j,l,1)=arrayr43d(i,j,l)
+                    if(abs(arrayr43d(i,j,l)-fillvalue)<small) smoke(i,j,l,1) = spval
+                  enddo
+                enddo
+              enddo
+            endif
+
+            ! model level dust
+            if(trim(fieldname)=='dust') then
+              !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,fv3dust,arrayr43d,spval,fillvalue)
+              do l=1,lm
+                do j=jsta,jend
+                  do i=ista, iend
+                    fv3dust(i,j,l,1)=arrayr43d(i,j,l)
+                    if(abs(arrayr43d(i,j,l)-fillvalue)<small) fv3dust(i,j,l,1) = spval
+                  enddo
+                enddo
+              enddo
+            endif
+
+            ! 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)
+              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
+                  enddo
+                enddo
+              enddo
+            endif
+
+            ! 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)
+              do l=1,lm
+                do j=jsta,jend
+                  do i=ista, iend
+                    coarsepm(i,j,l,1)=arrayr43d(i,j,l)
+                    if(abs(arrayr43d(i,j,l)-fillvalue)<small) coarsepm(i,j,l,1) = spval
+                  enddo
+                enddo
+              enddo
+            endif
+
+            ! 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)
+              do l=1,lm
+                do j=jsta,jend
+                  do i=ista, iend
+                    effri(i,j,l)=arrayr43d(i,j,l)
+                    if(abs(arrayr43d(i,j,l)-fillvalue)<small) effri(i,j,l) = spval
+                  enddo
+                enddo
+              enddo
+            endif
+
+            ! Thompson scheme cloud water effective radius
+            if(trim(fieldname)=='cleffr') then
+              !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,effrl,arrayr43d,spval,fillvalue)
+              do l=1,lm
+                do j=jsta,jend
+                  do i=ista, iend
+                    effrl(i,j,l)=arrayr43d(i,j,l)
+                    if(abs(arrayr43d(i,j,l)-fillvalue)<small) effrl(i,j,l) = spval
+                  enddo
+                enddo
+              enddo
+            endif
+
+            ! Thompson scheme cloud snow effective radius
+            if(trim(fieldname)=='cseffr') then
+              !$omp parallel do default(none) private(i,j,l) shared(lm,jsta,jend,ista,iend,effrs,arrayr43d,spval,fillvalue)
+              do l=1,lm
+                do j=jsta,jend
+                  do i=ista, iend
+                    effrs(i,j,l)=arrayr43d(i,j,l)
+                    if(abs(arrayr43d(i,j,l)-fillvalue)<small) effrs(i,j,l) = spval
+                  enddo
+                enddo
+              enddo
+            endif
+
 !3d fields
           endif
 
@@ -3384,6 +3756,23 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
         end do
       end do
 
+            ! 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
+               enddo
+             enddo
+      deallocate(extsmoke)
+      deallocate(extdust)
+      deallocate(buf)
+      deallocate(buf2)
+
 ! generate look up table for lifted parcel calculations
       thl    = 210.
       plq    = 70000.

From 487c190d72e1ff6d69fc80c4b5c492a371d74c88 Mon Sep 17 00:00:00 2001
From: Wen Meng <wen.meng@noaa.gov>
Date: Wed, 5 Apr 2023 00:11:58 +0000
Subject: [PATCH 03/11] Read MERRA2 aerosol

---
 io/post_fv3.F90           | 72 +++++++++++++++++++++++++++++++++++++--
 io/post_nems_routines.F90 |  9 +++--
 2 files changed, 75 insertions(+), 6 deletions(-)

diff --git a/io/post_fv3.F90 b/io/post_fv3.F90
index 56555c006..155f3aa1c 100644
--- a/io/post_fv3.F90
+++ b/io/post_fv3.F90
@@ -494,7 +494,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
 !     Jul 2019    J. Wang      Initial code
 !     Apr 2022    W. Meng      Unify set_postvars_gfs and
 !                               set_postvars_regional to set_postvars_fv3
-!     Apr 2023    W. Meng      Sync RRFS changes from off-line post
+!     Apr 2023    W. Meng      Sync RRFS and GFS changes from off-line post
 !
 !-----------------------------------------------------------------------
 !*** set up post fields from nmint_state
@@ -539,7 +539,9 @@ 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, aodtot, ebb, hwp,   &
+                             aod550,du_aod550,ss_aod550,su_aod550,oc_aod550,   &
+                             bc_aod550
       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, &
@@ -548,7 +550,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
                              tprec, tclod, trdlw, trdsw, tsrfc, tmaxmin, theat, &
                              ardlw, ardsw, asrfc, avrain, avcnvc, iSF_SURFACE_PHYSICS,&
                              td3d, idat, sdat, ifhr, ifmin, dt, nphs, dtq2, pt_tbl, &
-                             alsl, spl, ihrst, modelname, nsoil
+                             alsl, spl, ihrst, modelname, nsoil, rdaod
       use params_mod,  only: erad, dtr, capa, p1000, small
       use gridspec_mod,only: latstart, latlast, lonstart, lonlast, cenlon, cenlat, &
                              dxval, dyval, truelat2, truelat1, psmapf, cenlat,     &
@@ -2394,6 +2396,70 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
               enddo
             enddo
 
+            if(rdaod) then
+              ! MERRA2 aerosols 
+              if(trim(fieldname)=='aod550') then
+                !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,aod550,arrayr42d,fillValue)
+                do j=jsta,jend
+                  do i=ista, iend
+                    aod550(i,j) = arrayr42d(i,j)
+                    if(abs(arrayr42d(i,j)-fillvalue)<small) aod550(i,j) = spval
+                  enddo
+                enddo
+              endif
+
+              if(trim(fieldname)=='du_aod550') then
+                !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,du_aod550,arrayr42d,fillValue)
+                do j=jsta,jend
+                  do i=ista, iend
+                    du_aod550(i,j) = arrayr42d(i,j)
+                    if(abs(arrayr42d(i,j)-fillvalue)<small) du_aod550(i,j) = spval
+                  enddo
+                enddo
+              endif
+
+              if(trim(fieldname)=='ss_aod550') then
+                !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,ss_aod550,arrayr42d,fillValue)
+                do j=jsta,jend
+                  do i=ista, iend
+                    ss_aod550(i,j) = arrayr42d(i,j)
+                    if(abs(arrayr42d(i,j)-fillvalue)<small) ss_aod550(i,j) = spval
+                  enddo
+                enddo
+              endif
+ 
+              if(trim(fieldname)=='su_aod550') then
+                !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,su_aod550,arrayr42d,fillValue)
+                do j=jsta,jend
+                  do i=ista, iend
+                    su_aod550(i,j) = arrayr42d(i,j)
+                    if(abs(arrayr42d(i,j)-fillvalue)<small) su_aod550(i,j) = spval
+                  enddo
+                enddo
+              endif
+
+              if(trim(fieldname)=='oc_aod550') then
+                !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,oc_aod550,arrayr42d,fillValue)
+                do j=jsta,jend
+                  do i=ista, iend
+                    oc_aod550(i,j) = arrayr42d(i,j)
+                    if(abs(arrayr42d(i,j)-fillvalue)<small) oc_aod550(i,j) = spval
+                  enddo
+                enddo
+              endif
+
+              if(trim(fieldname)=='bc_aod550') then
+                !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,bc_aod550,arrayr42d,fillValue)
+                do j=jsta,jend
+                  do i=ista, iend
+                    bc_aod550(i,j) = arrayr42d(i,j)
+                    if(abs(arrayr42d(i,j)-fillvalue)<small) bc_aod550(i,j) = spval
+                  enddo
+                enddo
+              endif
+
+            endif !end rdaod
+
             ! inst  cloud top pressure
             if(trim(fieldname)=='prescnvclt') then
               !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,ptop,arrayr42d,fillValue)
diff --git a/io/post_nems_routines.F90 b/io/post_nems_routines.F90
index 9edf4365a..0c9e99424 100644
--- a/io/post_nems_routines.F90
+++ b/io/post_nems_routines.F90
@@ -225,7 +225,8 @@ subroutine read_postnmlt(kpo,kth,kpv,po,th,pv,post_namelist)
       use ctlblk_mod, only : komax,fileNameD3D,lsm,lsmp1,spl,spldef,  &
                              lsmdef,ALSL,me,d3d_on,gocart_on,hyb_sigp,&
                              pthresh,novegtype,ivegsrc,icu_physics,   &
-                             isf_surface_physics,modelname,submodelname
+                             isf_surface_physics,modelname,submodelname,&
+                             rdaod
       use upp_ifi_mod, only: write_ifi_debug_files
 !
 !    revision history:
@@ -242,7 +243,7 @@ subroutine read_postnmlt(kpo,kth,kpv,po,th,pv,post_namelist)
       integer l,k
 
       namelist/nampgb/kpo,po,kth,th,kpv,pv,popascal,d3d_on,gocart_on,  &
-                      hyb_sigp,write_ifi_debug_files
+                      hyb_sigp,write_ifi_debug_files,rdaod
       namelist/model_inputs/modelname,submodelname
 !---------------------------------------------------------------------
 !
@@ -259,6 +260,7 @@ subroutine read_postnmlt(kpo,kth,kpv,po,th,pv,post_namelist)
       d3d_on      = .false.
       gocart_on   = .false.
       popascal    = .false.
+      rdaod       = .false.
 !
       if (me == 0) print *,'post_namelist=',post_namelist
 !jw post namelist is using the same file itag as standalone post
@@ -268,7 +270,8 @@ subroutine read_postnmlt(kpo,kth,kpv,po,th,pv,post_namelist)
       close (nlunit)
       if (me == 0) then
         print*,'komax,kpo,kth,th,kpv,pv,popascal= ',komax,kpo            &
-     &  ,kth,th(1:kth),kpv,pv(1:kpv),popascal,' gocart_on=',gocart_on
+     &  ,kth,th(1:kth),kpv,pv(1:kpv),popascal,' gocart_on=',gocart_on    &
+     &  ,' rdaod=',rdaod
        endif
 !
 ! set up pressure level from POSTGPVARS or DEFAULT

From 5a7928a9371576d76406937df9802bee3de8966e Mon Sep 17 00:00:00 2001
From: Wen Meng <wen.meng@noaa.gov>
Date: Sat, 8 Apr 2023 03:19:31 +0000
Subject: [PATCH 04/11] Update upp revision to dccb321;Merge aerosol/chemical
 changes from Kate Zhang.

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

diff --git a/io/post_fv3.F90 b/io/post_fv3.F90
index 155f3aa1c..e235ec451 100644
--- a/io/post_fv3.F90
+++ b/io/post_fv3.F90
@@ -506,7 +506,10 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
                              qqr, qqs, cwm, qqi, qqw, qqg, omga, cfr, pmid,    &
                              q2, rlwtt, rswtt, tcucn, tcucns, train, el_pbl,   &
                              pint, exch_h, ref_10cm, qqni, qqnr, qqnwfa,       &
-                             qqnifa, effri, effrl, effrs, aextc55, taod5503d
+                             qqnifa, effri, effrl, effrs, aextc55, taod5503d,  &
+                             duem, dusd, dudp, duwt, dusv, ssem, sssd, ssdp,   &
+                             sswt, sssv, bcem, bcsd, bcdp, bcwt, bcsv, ocem,   &
+                             ocsd, ocdp, ocwt, ocsv
       use vrbls2d,     only: f, pd, sigt4, fis, pblh, ustar, z0, ths, qs, twbs,&
                              qwbs, avgcprate, cprate, avgprec, prec, lspa, sno,&
                              cldefi, th10, q10, tshltr, pshltr, albase,        &
@@ -541,7 +544,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
                              acgraup, graup_bucket, acfrain, frzrn_bucket,     &
                              ltg1_max, ltg2_max, ltg3_max, aodtot, ebb, hwp,   &
                              aod550,du_aod550,ss_aod550,su_aod550,oc_aod550,   &
-                             bc_aod550
+                             bc_aod550,maod
       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, &
@@ -550,7 +553,9 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
                              tprec, tclod, trdlw, trdsw, tsrfc, tmaxmin, theat, &
                              ardlw, ardsw, asrfc, avrain, avcnvc, iSF_SURFACE_PHYSICS,&
                              td3d, idat, sdat, ifhr, ifmin, dt, nphs, dtq2, pt_tbl, &
-                             alsl, spl, ihrst, modelname, nsoil, rdaod
+                             alsl, spl, ihrst, modelname, nsoil, rdaod, gocart_on,  &
+                             gccpp_on, nasa_on, d2d_chem, nbin_ss, nbin_bc, nbin_oc,&
+                             nbin_du
       use params_mod,  only: erad, dtr, capa, p1000, small
       use gridspec_mod,only: latstart, latlast, lonstart, lonlast, cenlon, cenlat, &
                              dxval, dyval, truelat2, truelat1, psmapf, cenlat,     &
@@ -597,7 +602,8 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
       real,dimension(:,:),allocatable :: dummy, p2d, t2d, q2d,  qs2d,  &
                              cw2d, cfr2d, buf, buf2
       real,dimension(:,:,:),allocatable :: extsmoke, extdust
-      character(len=80)              :: fieldname, wrtFBName, flatlon
+      character(len=80)              :: fieldname, wrtFBName, flatlon, &
+                                        VarName
       type(ESMF_Grid)                :: wrtGrid
       type(ESMF_Field)               :: theField
       type(ESMF_Field), allocatable  :: fcstField(:)
@@ -2460,6 +2466,350 @@ 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
+                if ( K == 1) VarName='duem001'
+                if ( K == 2) VarName='duem002'
+                if ( K == 3) VarName='duem003'
+                if ( K == 4) VarName='duem004'
+                if ( K == 5) VarName='duem005'
+
+                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
+                    do i=ista, iend
+                      duem(i,j,K) = arrayr42d(i,j)
+                      if( abs(arrayr42d(i,j)-fillValue) < small) duem(i,j,K) = spval
+                    enddo
+                  enddo
+                endif
+              enddo
+
+              do K = 1, nbin_du
+                if ( K == 1) VarName='dust1sd'
+                if ( K == 2) VarName='dust2sd'
+                if ( K == 3) VarName='dust3sd'
+                if ( K == 4) VarName='dust4sd'
+                if ( K == 5) VarName='dust5sd'
+
+                if(trim(fieldname)==VarName) then
+                  !$omp parallel do default(none) private(i,j,K) shared(jsta,jend,ista,iend,spval,dusd,arrayr42d,fillvalue)
+                  do j=jsta,jend
+                    do i=ista, iend
+                      dusd(i,j,K) = arrayr42d(i,j)
+                      if( abs(arrayr42d(i,j)-fillValue) < small) dusd(i,j,K) = spval
+                    enddo
+                  enddo
+                endif
+              enddo 
+
+              do K = 1, nbin_du
+                if ( K == 1) VarName='dust1dp'
+                if ( K == 2) VarName='dust2dp'
+                if ( K == 3) VarName='dust3dp' 
+                if ( K == 4) VarName='dust4dp'
+                if ( K == 5) VarName='dust5dp' 
+                             
+                if(trim(fieldname)==VarName) then
+                  !$omp parallel do default(none) private(i,j,K) shared(jsta,jend,ista,iend,spval,dudp,arrayr42d,fillvalue)
+                  do j=jsta,jend
+                    do i=ista, iend
+                      dudp(i,j,K) = arrayr42d(i,j)
+                      if( abs(arrayr42d(i,j)-fillValue) < small) dudp(i,j,K) = spval
+                    enddo    
+                  enddo      
+                endif        
+              enddo
+
+              do K = 1, nbin_du
+                if ( K == 1) VarName='dust1wtl'
+                if ( K == 2) VarName='dust2wtl'
+                if ( K == 3) VarName='dust3wtl'
+                if ( K == 4) VarName='dust4wtl'
+                if ( K == 5) VarName='dust5wtl'
+
+                if(trim(fieldname)==VarName) then
+                  !$omp parallel do default(none) private(i,j,K) shared(jsta,jend,ista,iend,spval,duwt,arrayr42d,fillvalue)
+                  do j=jsta,jend
+                    do i=ista, iend
+                      duwt(i,j,K) = arrayr42d(i,j)
+                      if( abs(arrayr42d(i,j)-fillValue) < small) duwt(i,j,K) = spval
+                    enddo
+                  enddo
+                endif
+              enddo
+
+              do K = 1, nbin_du
+                if ( K == 1) VarName='dust1wtc'
+                if ( K == 2) VarName='dust2wtc'
+                if ( K == 3) VarName='dust3wtc'
+                if ( K == 4) VarName='dust4wtc'
+                if ( K == 5) VarName='dust5wtc'
+                             
+                if(trim(fieldname)==VarName) then
+                  !$omp parallel do default(none) private(i,j,K) shared(jsta,jend,ista,iend,spval,dusv,arrayr42d,fillvalue)
+                  do j=jsta,jend
+                    do i=ista, iend
+                      dusv(i,j,K) = arrayr42d(i,j)
+                      if( abs(arrayr42d(i,j)-fillValue) < small) dusv(i,j,K) = spval
+                    enddo    
+                  enddo      
+                endif        
+              enddo
+
+              do K = 1, nbin_ss
+                if ( K == 1) VarName='ssem001'
+                if ( K == 2) VarName='ssem002'
+                if ( K == 3) VarName='ssem003'
+                if ( K == 4) VarName='ssem004'
+                if ( K == 5) VarName='ssem005'
+
+                if(trim(fieldname)==VarName) then
+                  !$omp parallel do default(none) private(i,j,K) shared(jsta,jend,ista,iend,spval,ssem,arrayr42d,fillvalue)
+                  do j=jsta,jend
+                    do i=ista, iend
+                      ssem(i,j,K) = arrayr42d(i,j)
+                      if( abs(arrayr42d(i,j)-fillValue) < small) ssem(i,j,K) = spval
+                    enddo
+                  enddo
+                endif
+              enddo
+
+              do K = 1, nbin_ss
+                if ( K == 1) VarName='seas1sd'
+                if ( K == 2) VarName='seas2sd'
+                if ( K == 3) VarName='seas3sd'
+                if ( K == 4) VarName='seas4sd'
+                if ( K == 5) VarName='seas5sd'
+
+                if(trim(fieldname)==VarName) then
+                  !$omp parallel do default(none) private(i,j,K) shared(jsta,jend,ista,iend,spval,sssd,arrayr42d,fillvalue)
+                  do j=jsta,jend
+                    do i=ista, iend
+                      sssd(i,j,K) = arrayr42d(i,j)
+                      if( abs(arrayr42d(i,j)-fillValue) < small) sssd(i,j,K) = spval
+                    enddo
+                  enddo
+                endif
+              enddo
+
+              do K = 1, nbin_ss
+                if ( K == 1) VarName='seas1dp'
+                if ( K == 2) VarName='seas2dp'
+                if ( K == 3) VarName='seas3dp'
+                if ( K == 4) VarName='seas4dp'
+                if ( K == 5) VarName='seas5dp'
+
+                if(trim(fieldname)==VarName) then
+                  !$omp parallel do default(none) private(i,j,K) shared(jsta,jend,ista,iend,spval,ssdp,arrayr42d,fillvalue)
+                  do j=jsta,jend
+                    do i=ista, iend
+                      ssdp(i,j,K) = arrayr42d(i,j)
+                      if( abs(arrayr42d(i,j)-fillValue) < small) ssdp(i,j,K) = spval
+                    enddo
+                  enddo
+                 endif
+               enddo
+
+              do K = 1, nbin_ss
+                if ( K == 1) VarName='seas1wt'
+                if ( K == 2) VarName='seas2wt'
+                if ( K == 3) VarName='seas3wt'
+                if ( K == 4) VarName='seas4wt'
+                if ( K == 5) VarName='seas5wt'
+
+                if(trim(fieldname)==VarName) then
+                  !$omp parallel do default(none) private(i,j,K) shared(jsta,jend,ista,iend,spval,sswt,arrayr42d,fillvalue)
+                  do j=jsta,jend
+                    do i=ista, iend
+                      sswt(i,j,K) = arrayr42d(i,j)
+                      if( abs(arrayr42d(i,j)-fillValue) < small) sswt(i,j,K) = spval
+                    enddo
+                  enddo
+                 endif
+               enddo
+
+              do K = 1, nbin_ss
+                if ( K == 1) VarName='seas1wtc'
+                if ( K == 2) VarName='seas2wtc'
+                if ( K == 3) VarName='seas3wtc'
+                if ( K == 4) VarName='seas4wtc'
+                if ( K == 5) VarName='seas5wtc'
+
+                if(trim(fieldname)==VarName) then
+                  !$omp parallel do default(none) private(i,j,K) shared(jsta,jend,ista,iend,spval,sssv,arrayr42d,fillvalue)
+                  do j=jsta,jend
+                    do i=ista, iend
+                      sssv(i,j,K) = arrayr42d(i,j)
+                      if( abs(arrayr42d(i,j)-fillValue) < small) sssv(i,j,K) = spval
+                    enddo
+                  enddo
+                 endif
+               enddo
+
+              do K = 1, nbin_bc
+                if ( K == 1) VarName='bceman'
+                if ( K == 2) VarName='bcembb'
+
+                if(trim(fieldname)==VarName) then
+                  !$omp parallel do default(none) private(i,j,K) shared(jsta,jend,ista,iend,spval,bcem,arrayr42d,fillvalue)
+                  do j=jsta,jend
+                    do i=ista, iend
+                      bcem(i,j,K) = arrayr42d(i,j)
+                      if( abs(arrayr42d(i,j)-fillValue) < small) bcem(i,j,K) = spval
+                    enddo
+                  enddo
+                 endif
+               enddo
+
+              do K = 1, nbin_bc
+                if ( K == 1) VarName='bc1sd'
+                if ( K == 2) VarName='bc2sd'
+
+                if(trim(fieldname)==VarName) then
+                  !$omp parallel do default(none) private(i,j,K) shared(jsta,jend,ista,iend,spval,bcsd,arrayr42d,fillvalue)
+                  do j=jsta,jend
+                    do i=ista, iend
+                      bcsd(i,j,K) = arrayr42d(i,j)
+                      if( abs(arrayr42d(i,j)-fillValue) < small) bcsd(i,j,K) = spval
+                    enddo
+                  enddo
+                 endif
+              enddo
+
+              do K = 1, nbin_bc
+                if ( K == 1) VarName='bc1dp'
+                if ( K == 2) VarName='bc2dp'
+
+                if(trim(fieldname)==VarName) then
+                  !$omp parallel do default(none) private(i,j,K) shared(jsta,jend,ista,iend,spval,bcdp,arrayr42d,fillvalue)
+                  do j=jsta,jend
+                    do i=ista, iend
+                      bcdp(i,j,K) = arrayr42d(i,j)
+                      if( abs(arrayr42d(i,j)-fillValue) < small) bcdp(i,j,K) = spval
+                    enddo
+                  enddo
+                endif
+              enddo
+
+              do K = 1, nbin_bc
+                if ( K == 1) VarName='bc1wtl'
+                if ( K == 2) VarName='bc2wtl'
+
+                if(trim(fieldname)==VarName) then
+                  !$omp parallel do default(none) private(i,j,K) shared(jsta,jend,ista,iend,spval,bcwt,arrayr42d,fillvalue)
+                  do j=jsta,jend
+                    do i=ista, iend
+                      bcwt(i,j,K) = arrayr42d(i,j)
+                      if( abs(arrayr42d(i,j)-fillValue) < small) bcwt(i,j,K) = spval
+                    enddo
+                  enddo
+                endif
+              enddo
+
+              do K = 1, nbin_bc
+                if ( K == 1) VarName='bc1wtc'
+                if ( K == 2) VarName='bc2wtc'
+
+                if(trim(fieldname)==VarName) then
+                  !$omp parallel do default(none) private(i,j,K) shared(jsta,jend,ista,iend,spval,bcsv,arrayr42d,fillvalue)
+                  do j=jsta,jend
+                    do i=ista, iend
+                      bcsv(i,j,K) = arrayr42d(i,j)
+                      if( abs(arrayr42d(i,j)-fillValue) < small) bcsv(i,j,K) = spval
+                    enddo
+                  enddo
+                endif
+              enddo
+
+              do K = 1, nbin_oc
+                if ( K == 1) VarName='oceman'
+                if ( K == 2) VarName='ocembb'
+
+                if(trim(fieldname)==VarName) then
+                  !$omp parallel do default(none) private(i,j,K) shared(jsta,jend,ista,iend,spval,ocem,arrayr42d,fillvalue)
+                  do j=jsta,jend
+                    do i=ista, iend
+                      ocem(i,j,K) = arrayr42d(i,j)
+                      if( abs(arrayr42d(i,j)-fillValue) < small) ocem(i,j,K) = spval
+                    enddo
+                  enddo
+                endif
+              enddo
+
+              do K = 1, nbin_oc
+                if ( K == 1) VarName='oc1sd'
+                if ( K == 2) VarName='oc2sd'
+
+                if(trim(fieldname)==VarName) then
+                  !$omp parallel do default(none) private(i,j,K) shared(jsta,jend,ista,iend,spval,ocsd,arrayr42d,fillvalue)
+                  do j=jsta,jend
+                    do i=ista, iend
+                      ocsd(i,j,K) = arrayr42d(i,j)
+                      if( abs(arrayr42d(i,j)-fillValue) < small) ocsd(i,j,K) = spval
+                    enddo
+                  enddo
+                endif
+              enddo
+
+              do K = 1, nbin_oc
+                if ( K == 1) VarName='oc1dp'
+                if ( K == 2) VarName='oc2dp'
+
+                if(trim(fieldname)==VarName) then
+                  !$omp parallel do default(none) private(i,j,K) shared(jsta,jend,ista,iend,spval,ocdp,arrayr42d,fillvalue)
+                  do j=jsta,jend
+                    do i=ista, iend
+                      ocdp(i,j,K) = arrayr42d(i,j)
+                      if( abs(arrayr42d(i,j)-fillValue) < small) ocdp(i,j,K) = spval
+                    enddo
+                  enddo
+                endif
+              enddo
+
+              do K = 1, nbin_oc
+                if ( K == 1) VarName='oc1wtl'
+                if ( K == 2) VarName='oc2wtl'
+
+                if(trim(fieldname)==VarName) then
+                  !$omp parallel do default(none) private(i,j,K) shared(jsta,jend,ista,iend,spval,ocwt,arrayr42d,fillvalue)
+                  do j=jsta,jend
+                    do i=ista, iend
+                      ocwt(i,j,K) = arrayr42d(i,j)
+                      if( abs(arrayr42d(i,j)-fillValue) < small) ocwt(i,j,K) = spval
+                    enddo
+                  enddo
+                endif
+              enddo
+
+              do K = 1, nbin_oc
+                if ( K == 1) VarName='oc1wtc'
+                if ( K == 2) VarName='oc2wtc'
+
+                if(trim(fieldname)==VarName) then
+                  !$omp parallel do default(none) private(i,j,K) shared(jsta,jend,ista,iend,spval,ocsv,arrayr42d,fillvalue)
+                  do j=jsta,jend
+                    do i=ista, iend
+                      ocsv(i,j,K) = arrayr42d(i,j)
+                      if( abs(arrayr42d(i,j)-fillValue) < small) ocsv(i,j,K) = spval
+                    enddo
+                  enddo
+                endif
+              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
+                  enddo
+                endif
+
+            endif !end gocart_on
+
             ! inst  cloud top pressure
             if(trim(fieldname)=='prescnvclt') then
               !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,ptop,arrayr42d,fillValue)
diff --git a/io/post_nems_routines.F90 b/io/post_nems_routines.F90
index 0c9e99424..337139b1d 100644
--- a/io/post_nems_routines.F90
+++ b/io/post_nems_routines.F90
@@ -226,7 +226,7 @@ subroutine read_postnmlt(kpo,kth,kpv,po,th,pv,post_namelist)
                              lsmdef,ALSL,me,d3d_on,gocart_on,hyb_sigp,&
                              pthresh,novegtype,ivegsrc,icu_physics,   &
                              isf_surface_physics,modelname,submodelname,&
-                             rdaod
+                             rdaod,d2d_chem,nasa_on,gccpp_on
       use upp_ifi_mod, only: write_ifi_debug_files
 !
 !    revision history:
@@ -243,7 +243,8 @@ subroutine read_postnmlt(kpo,kth,kpv,po,th,pv,post_namelist)
       integer l,k
 
       namelist/nampgb/kpo,po,kth,th,kpv,pv,popascal,d3d_on,gocart_on,  &
-                      hyb_sigp,write_ifi_debug_files,rdaod
+                      hyb_sigp,write_ifi_debug_files,rdaod,nasa_on,gccpp_on, &
+                      d2d_chem
       namelist/model_inputs/modelname,submodelname
 !---------------------------------------------------------------------
 !
@@ -261,6 +262,9 @@ subroutine read_postnmlt(kpo,kth,kpv,po,th,pv,post_namelist)
       gocart_on   = .false.
       popascal    = .false.
       rdaod       = .false.
+      nasa_on     = .false.
+      gccpp_on    = .false.
+      d2d_chem    = .false.
 !
       if (me == 0) print *,'post_namelist=',post_namelist
 !jw post namelist is using the same file itag as standalone post
diff --git a/upp b/upp
index ce35e5f8d..dccb32176 160000
--- a/upp
+++ b/upp
@@ -1 +1 @@
-Subproject commit ce35e5f8d2f44a95b71b6d142cd854f73188d074
+Subproject commit dccb32176930676ef2a258eb65571ab4e3f7e7a4

From f535eb632ce1a20388a07f063404ccd0ccb1048e Mon Sep 17 00:00:00 2001
From: Wen Meng <wen.meng@noaa.gov>
Date: Wed, 12 Apr 2023 02:13:36 +0000
Subject: [PATCH 05/11] Incorporate changes of reading chemical fields from
 Kate Zhang.

---
 io/post_fv3.F90 | 481 +++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 474 insertions(+), 7 deletions(-)

diff --git a/io/post_fv3.F90 b/io/post_fv3.F90
index e235ec451..e98b30ed2 100644
--- a/io/post_fv3.F90
+++ b/io/post_fv3.F90
@@ -40,7 +40,7 @@ subroutine post_run_fv3(wrt_int_state,grid_id,mype,mpicomp,lead_write, &
 !-----------------------------------------------------------------------
 !
       use ctlblk_mod, only : komax,ifhr,ifmin,modelname,datapd,fld_info, &
-                             npset,grib,gocart_on,jsta,  &
+                             npset,grib,jsta,  &
                              jend,ista,iend, im, nsoil, filenameflat,numx
       use gridspec_mod, only : maptype, gridtype,latstart,latlast,       &
                                lonstart,lonlast
@@ -501,7 +501,8 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
 !-----------------------------------------------------------------------
 !
       use esmf
-      use vrbls4d,     only: dust, smoke, fv3dust, coarsepm
+      use vrbls4d,     only: dust, smoke, fv3dust, coarsepm, SALT, SUSO, SOOT, &
+                             WASO,no3,nh4, PP25, PP10
       use vrbls3d,     only: t, q, uh, vh, wh, alpint, dpres, zint, zmid, o3,  &
                              qqr, qqs, cwm, qqi, qqw, qqg, omga, cfr, pmid,    &
                              q2, rlwtt, rswtt, tcucn, tcucns, train, el_pbl,   &
@@ -509,7 +510,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
                              qqnifa, effri, effrl, effrs, aextc55, taod5503d,  &
                              duem, dusd, dudp, duwt, dusv, ssem, sssd, ssdp,   &
                              sswt, sssv, bcem, bcsd, bcdp, bcwt, bcsv, ocem,   &
-                             ocsd, ocdp, ocwt, ocsv
+                             ocsd, ocdp, ocwt, ocsv, rhomid
       use vrbls2d,     only: f, pd, sigt4, fis, pblh, ustar, z0, ths, qs, twbs,&
                              qwbs, avgcprate, cprate, avgprec, prec, lspa, sno,&
                              cldefi, th10, q10, tshltr, pshltr, albase,        &
@@ -544,19 +545,22 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
                              acgraup, graup_bucket, acfrain, frzrn_bucket,     &
                              ltg1_max, ltg2_max, ltg3_max, aodtot, ebb, hwp,   &
                              aod550,du_aod550,ss_aod550,su_aod550,oc_aod550,   &
-                             bc_aod550,maod
+                             bc_aod550,maod,                                   &
+                             dustpm10, dustcb, bccb, occb, sulfcb, sscb,       &
+                             dustallcb, ssallcb, dustpm, sspm, pp25cb, pp10cb, &
+                             no3cb, nh4cb, dusmass, ducmass, dusmass25,ducmass25
       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, &
-                             ista, iend, ista_2l, iend_2u, ista_m,iend_m, &
+                             ista, iend, ista_2l, iend_2u, ista_m,iend_m,qmin, &
                              lsm, pt, imp_physics, spval, mpi_comm_comp, gdsdegr,  &
                              tprec, tclod, trdlw, trdsw, tsrfc, tmaxmin, theat, &
                              ardlw, ardsw, asrfc, avrain, avcnvc, iSF_SURFACE_PHYSICS,&
                              td3d, idat, sdat, ifhr, ifmin, dt, nphs, dtq2, pt_tbl, &
                              alsl, spl, ihrst, modelname, nsoil, rdaod, gocart_on,  &
                              gccpp_on, nasa_on, d2d_chem, nbin_ss, nbin_bc, nbin_oc,&
-                             nbin_du
-      use params_mod,  only: erad, dtr, capa, p1000, small
+                             nbin_du,nbin_su, nbin_no3, nbin_nh4
+      use params_mod,  only: erad, dtr, capa, p1000, small,h1, d608, pi, rd
       use gridspec_mod,only: latstart, latlast, lonstart, lonlast, cenlon, cenlat, &
                              dxval, dyval, truelat2, truelat1, psmapf, cenlat,     &
                              lonstartv, lonlastv, cenlonv, latstartv, latlastv,    &
@@ -3962,6 +3966,469 @@ 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
+
+        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)=='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)=='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)=='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
+
+        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
+
+        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)=='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)=='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
+
+        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 ! 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(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
+
+        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
+
+              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)
+           endif !gocart_on or gccpp_on
+
+           if (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 !end gocart_on, nasa_on
+
 !3d fields
           endif
 

From 2e46deaa25f6a6eb987282d0c70515c6e86de64d Mon Sep 17 00:00:00 2001
From: Wen Meng <wen.meng@noaa.gov>
Date: Wed, 12 Apr 2023 19:43:28 +0000
Subject: [PATCH 06/11] Update based on Dusan's comment

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

diff --git a/io/post_fv3.F90 b/io/post_fv3.F90
index e98b30ed2..1e6bb4266 100644
--- a/io/post_fv3.F90
+++ b/io/post_fv3.F90
@@ -86,7 +86,6 @@ subroutine post_run_fv3(wrt_int_state,grid_id,mype,mpicomp,lead_write, &
 
       grib      = "grib2"
       gridtype  = "A"
-      !nsoil     = 4
       nsoil     = wrt_int_state%nsoil
       nwtpg     = wrt_int_state%petcount
       jts       = wrt_int_state%out_grid_info(grid_id)%j_start     !<-- Starting J of this write task's subsection
@@ -818,12 +817,10 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
 !                   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(smstav,sfcevp,acsnow,acsnom,qz0,uz0,vz0,maxrhshltr,minrhshltr)
+!$omp& shared(sfcevp,acsnom,qz0,uz0,vz0,maxrhshltr,minrhshltr)
       do j=jsta_2l,jend_2u
         do i=ista_2l,iend_2u
-          !smstav(i,j) = spval
           sfcevp(i,j) = spval
-          !acsnow(i,j) = spval
           acsnom(i,j) = spval
           qz0(i,j)    = spval
           uz0(i,j)    = spval
@@ -2472,11 +2469,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
 
             if ((gocart_on .or. gccpp_on) .and. d2d_chem) then
               do K = 1, nbin_du
-                if ( K == 1) VarName='duem001'
-                if ( K == 2) VarName='duem002'
-                if ( K == 3) VarName='duem003'
-                if ( K == 4) VarName='duem004'
-                if ( K == 5) VarName='duem005'
+                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)

From 2951068995663ebd95b2ca1d2f5e653648ffcf2b Mon Sep 17 00:00:00 2001
From: Wen Meng <wen.meng@noaa.gov>
Date: Thu, 13 Apr 2023 14:19:34 +0000
Subject: [PATCH 07/11] Limit new local variable allocation for RRFS.

---
 io/post_fv3.F90 | 21 +++++++++++++++++----
 1 file changed, 17 insertions(+), 4 deletions(-)

diff --git a/io/post_fv3.F90 b/io/post_fv3.F90
index 1e6bb4266..f8864e305 100644
--- a/io/post_fv3.F90
+++ b/io/post_fv3.F90
@@ -695,10 +695,13 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
       nphs = 2.
       dt   = dtq2/nphs
 
-      allocate(extsmoke(ista:iend,jsta:jend,lm))
-      allocate(extdust(ista:iend,jsta:jend,lm))
-      allocate(buf(ista:iend,jsta:jend))
-      allocate(buf2(ista:iend,jsta:jend))
+      !Allocate for regional models only
+      if(modelname=='FV3R') then
+        allocate(extsmoke(ista:iend,jsta:jend,lm))
+        allocate(extdust(ista:iend,jsta:jend,lm))
+        allocate(buf(ista:iend,jsta:jend))
+        allocate(buf2(ista:iend,jsta:jend))
+      endif
 
 !
 ! GFS does not have convective cloud efficiency
@@ -2343,6 +2346,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
               enddo
             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,buf,arrayr42d,fillvalue,spval)
@@ -2372,6 +2376,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
                  endif
                enddo
              enddo
+             endif !FV3R
 
             !sndepac
             if(trim(fieldname)=='snacc_land') then
@@ -3881,6 +3886,8 @@ 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)
@@ -3907,6 +3914,8 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
               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)
@@ -4632,6 +4641,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
         end do
       end do
 
+      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
@@ -4644,11 +4654,14 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
                  enddo
                enddo
              enddo
+
       deallocate(extsmoke)
       deallocate(extdust)
       deallocate(buf)
       deallocate(buf2)
 
+      endif !end FV3R
+
 ! generate look up table for lifted parcel calculations
       thl    = 210.
       plq    = 70000.

From 78cfb352b3b6955f29e7975fb4fbdecf13adf348 Mon Sep 17 00:00:00 2001
From: Wen Meng <wen.meng@noaa.gov>
Date: Thu, 13 Apr 2023 17:22:29 +0000
Subject: [PATCH 08/11] Modify accsnow and sndepac calculation location based
 on Dusan's comment.

---
 io/post_fv3.F90 | 100 ++++++++++++++++++++++++++----------------------
 1 file changed, 54 insertions(+), 46 deletions(-)

diff --git a/io/post_fv3.F90 b/io/post_fv3.F90
index f8864e305..672ff3388 100644
--- a/io/post_fv3.F90
+++ b/io/post_fv3.F90
@@ -603,7 +603,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, buf, buf2
+                             cw2d, cfr2d, accswe_ice, accswe_land, &
+                             snacc_land, snacc_ice
       real,dimension(:,:,:),allocatable :: extsmoke, extdust
       character(len=80)              :: fieldname, wrtFBName, flatlon, &
                                         VarName
@@ -699,8 +700,10 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
       if(modelname=='FV3R') then
         allocate(extsmoke(ista:iend,jsta:jend,lm))
         allocate(extdust(ista:iend,jsta:jend,lm))
-        allocate(buf(ista:iend,jsta:jend))
-        allocate(buf2(ista:iend,jsta:jend))
+        allocate(accswe_ice(ista:iend,jsta:jend))
+        allocate(accswe_land(ista:iend,jsta:jend))
+        allocate(snacc_ice(ista:iend,jsta:jend))
+        allocate(snacc_land(ista:iend,jsta:jend))
       endif
 
 !
@@ -2349,64 +2352,45 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
             if(modelname=='FV3R')then
             !acsnow
             if(trim(fieldname)=='accswe_land') then
-              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,buf,arrayr42d,fillvalue,spval)
+              !$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
-                  buf(i,j) = arrayr42d(i,j)
-                  if(abs(arrayr42d(i,j)-fillvalue)<small) buf(i,j) = spval
+                  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,buf2,arrayr42d,fillvalue,spval)
+              !$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
-                  buf2(i,j) = arrayr42d(i,j)
-                  if(abs(arrayr42d(i,j)-fillvalue)<small) buf2(i,j) = spval
+                  accswe_ice(i,j) = arrayr42d(i,j)
+                  if(abs(arrayr42d(i,j)-fillvalue)<small) accswe_ice(i,j) = spval
                 enddo
               enddo
             endif
-            ! !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,buf,buf2,acsnow)
-             do j=jsta,jend
-               do i=ista, iend
-                 if(buf(i,j)<spval .and. buf2(i,j)<spval) then
-                   acsnow(i,j) = buf(i,j) + buf2(i,j)
-                 else
-                   acsnow(i,j) = spval
-                 endif
-               enddo
-             enddo
-             endif !FV3R
 
             !sndepac
             if(trim(fieldname)=='snacc_land') then
-              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,buf,arrayr42d,fillvalue,spval)
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,snacc_land,arrayr42d,fillvalue,spval)
               do j=jsta,jend
                 do i=ista, iend
-                  buf(i,j) = arrayr42d(i,j)
-                  if(abs(arrayr42d(i,j)-fillvalue)<small) buf(i,j) = spval
+                  snacc_land(i,j) = arrayr42d(i,j)
+                  if(abs(arrayr42d(i,j)-fillvalue)<small) snacc_land(i,j) = spval
                 enddo
               enddo
             endif
             if(trim(fieldname)=='snacc_ice') then
-              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,buf2,arrayr42d,fillvalue,spval)
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,snacc_ice,arrayr42d,fillvalue,spval)
               do j=jsta,jend
                 do i=ista, iend
-                  buf2(i,j) = arrayr42d(i,j)
-                  if(abs(arrayr42d(i,j)-fillvalue)<small) buf2(i,j) = spval
+                  snacc_ice(i,j) = arrayr42d(i,j)
+                  if(abs(arrayr42d(i,j)-fillvalue)<small) snacc_ice(i,j) = spval
                 enddo
               enddo
             endif
-            !!$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,buf,buf2,sndepac)
-            do j=jsta,jend
-              do i=ista, iend
-                if(buf(i,j)<spval .and. buf2(i,j)<spval) then
-                  sndepac(i,j) = buf(i,j) + buf2(i,j)
-                else
-                  sndepac(i,j) = spval
-                endif
-              enddo
-            enddo
+
+            endif !FV3R
 
             if(rdaod) then
               ! MERRA2 aerosols 
@@ -4642,23 +4626,47 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
       end do
 
       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
+        ! 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
+            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_ice(i,j)<spval .and. accswe_land(i,j)<spval) then
+                   acsnow(i,j) = accswe_ice(i,j) + accswe_land(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_ice(i,j)<spval .and. snacc_land(i,j)<spval) then
+                  sndepac(i,j) = snacc_ice(i,j) + snacc_land(i,j)
+                else
+                  sndepac(i,j) = spval
+                endif
+              enddo
+            enddo
+
       deallocate(extsmoke)
       deallocate(extdust)
       deallocate(buf)
       deallocate(buf2)
+      deallocate(buf3)
+      deallocate(buf4)
 
       endif !end FV3R
 

From 113dac23e9b0a1f7eb6744c36b243dc124076ee5 Mon Sep 17 00:00:00 2001
From: Wen Meng <wen.meng@noaa.gov>
Date: Thu, 13 Apr 2023 18:06:30 +0000
Subject: [PATCH 09/11] Bug fix

---
 io/post_fv3.F90 | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/io/post_fv3.F90 b/io/post_fv3.F90
index 672ff3388..2f2409328 100644
--- a/io/post_fv3.F90
+++ b/io/post_fv3.F90
@@ -4663,10 +4663,10 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
 
       deallocate(extsmoke)
       deallocate(extdust)
-      deallocate(buf)
-      deallocate(buf2)
-      deallocate(buf3)
-      deallocate(buf4)
+      deallocate(accswe_ice)
+      deallocate(accswe_land)
+      deallocate(snacc_ice)
+      deallocate(snacc_land)
 
       endif !end FV3R
 

From 246acf25b8b3b34f2dadca7d3066f243439d22d5 Mon Sep 17 00:00:00 2001
From: Wen Meng <wen.meng@noaa.gov>
Date: Thu, 13 Apr 2023 19:44:15 +0000
Subject: [PATCH 10/11] 1)pick up minor change from Jun Wang; 2)Changes based
 on Jun and Dusan's   comments.

---
 io/module_wrt_grid_comp.F90 |  2 +-
 io/post_fv3.F90             | 76 +++++++++++++++++++------------------
 2 files changed, 41 insertions(+), 37 deletions(-)

diff --git a/io/module_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90
index 8d6c44b17..6a3b8ee01 100644
--- a/io/module_wrt_grid_comp.F90
+++ b/io/module_wrt_grid_comp.F90
@@ -2295,7 +2295,7 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc)
 !** write out log file
 !
     if (mype == lead_write_task) then
-      open(newunit=nolog,file='atm.logf'//trim(cfhour),form='FORMATTED')
+      open(newunit=nolog,file='log.atm.f'//trim(cfhour),form='FORMATTED') 
         write(nolog,100)nfhour,idate(1:6)
 100     format(' completed fv3gfs fhour=',f10.3,2x,6(i4,2x))
       close(nolog)
diff --git a/io/post_fv3.F90 b/io/post_fv3.F90
index 2f2409328..9734817ed 100644
--- a/io/post_fv3.F90
+++ b/io/post_fv3.F90
@@ -4353,8 +4353,8 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
         do j=jsta,jend
           do i=ista,iend
 
-              TV = T(I,J,L) * (H1+D608*MAX(Q(I,J,L),QMIN))
-              RHOMID(I,J,L) = PMID(I,J,L) / (RD*TV)
+            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)
@@ -4378,36 +4378,35 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
             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)
-           endif !gocart_on or gccpp_on
-
-           if (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)
+             !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
@@ -4415,6 +4414,7 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
 
         endif !end gocart_on, nasa_on
 
+
 !3d fields
           endif
 
@@ -4642,8 +4642,10 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
              !$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_ice(i,j)<spval .and. accswe_land(i,j)<spval) then
-                   acsnow(i,j) = accswe_ice(i,j) + accswe_land(i,j)
+                 if(accswe_ice(i,j)<spval) then
+                   acsnow(i,j) = accswe_ice(i,j) 
+                 elseif(accswe_land(i,j)<spval) then
+                   acsnow(i,j) = accswe_land(i,j)
                  else
                    acsnow(i,j) = spval
                  endif
@@ -4653,8 +4655,10 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
             !$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_ice(i,j)<spval .and. snacc_land(i,j)<spval) then
-                  sndepac(i,j) = snacc_ice(i,j) + snacc_land(i,j)
+                if(snacc_ice(i,j)<spval) then
+                  sndepac(i,j) = snacc_ice(i,j) 
+                elseif(snacc_land(i,j)<spval) then
+                  sndepac(i,j) = snacc_land(i,j)
                 else
                   sndepac(i,j) = spval
                 endif

From 5b930585699ff75b691fe0c72b32a50dfe7158a2 Mon Sep 17 00:00:00 2001
From: Wen Meng <wen.meng@noaa.gov>
Date: Thu, 13 Apr 2023 20:02:51 +0000
Subject: [PATCH 11/11] modify acsnow and sndepac calculation based on Jun's
 comment

---
 io/post_fv3.F90 | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/io/post_fv3.F90 b/io/post_fv3.F90
index 9734817ed..7cc6ab45e 100644
--- a/io/post_fv3.F90
+++ b/io/post_fv3.F90
@@ -4642,10 +4642,10 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
              !$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_ice(i,j)<spval) then
-                   acsnow(i,j) = accswe_ice(i,j) 
-                 elseif(accswe_land(i,j)<spval) then
-                   acsnow(i,j) = accswe_land(i,j)
+                 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
@@ -4655,10 +4655,10 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
             !$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_ice(i,j)<spval) then
-                  sndepac(i,j) = snacc_ice(i,j) 
-                elseif(snacc_land(i,j)<spval) then
-                  sndepac(i,j) = snacc_land(i,j)
+                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