From 8a1c0e725e4a318ea72cd303c7e274f4b46ff52c Mon Sep 17 00:00:00 2001
From: emilyhcliu <36091766+emilyhcliu@users.noreply.github.com>
Date: Thu, 22 Aug 2024 14:20:36 -0400
Subject: [PATCH] Recover the capability of handling model fields from
 operation gfs.v16.3 (#785)

---
 src/gsi/general_read_gfsatm.f90 | 177 ++++++++++++++++----------------
 1 file changed, 89 insertions(+), 88 deletions(-)

diff --git a/src/gsi/general_read_gfsatm.f90 b/src/gsi/general_read_gfsatm.f90
index c34b269950..e1a5406a13 100755
--- a/src/gsi/general_read_gfsatm.f90
+++ b/src/gsi/general_read_gfsatm.f90
@@ -3588,99 +3588,100 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
           if ( icount == icm .or. k==nlevs) then
              call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
                   g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
-      endif
-   enddo ! do k=1,nlevs
-
-   do k=1,nlevs
-      icount=icount+1
-      iflag(icount)=15
-      ilev(icount)=k
-      kr = levs+1-k ! netcdf is top to bottom, need to flip
+          endif
+       enddo ! do k=1,nlevs
 
-      if (mype==mype_use(icount)) then
-         call read_vardata(filges, 'nccice', rwork3d0, nslice=kr, slicedim=3)
-         ! cloud ice water number concentration.
-         if ( diff_res ) then
-            grid_b=rwork3d0(:,:,1)
-            vector(1)=.false.
-            call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
-            call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector)
-            do kk=1,grd%itotsub
-               i=grd%ltosi_s(kk)
-               j=grd%ltosj_s(kk)
-               work(kk)=grid2(i,j,1)
-            enddo
-         else
-            grid=rwork3d0(:,:,1)
-            call general_fill_ns(grd,grid,work)
-         endif
-      endif
-      if ( icount == icm .or. k==nlevs ) then
-         call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
-              g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
-      endif
-   enddo ! do k=1,nlevs
+       ! Read fields specific to Thompson microphysics
+       if (imp_physics == 8) then
+          do k=1,nlevs
+             icount=icount+1
+             iflag(icount)=15
+             ilev(icount)=k
+             kr = levs+1-k ! netcdf is top to bottom, need to flip
 
-   do k=1,nlevs
-      icount=icount+1
-      iflag(icount)=16
-      ilev(icount)=k
-      kr = levs+1-k ! netcdf is top to bottom, need to flip
+             if (mype==mype_use(icount)) then
+                call read_vardata(filges, 'nccice', rwork3d0, nslice=kr, slicedim=3)
+                ! cloud ice water number concentration.
+                if ( diff_res ) then
+                   grid_b=rwork3d0(:,:,1)
+                   vector(1)=.false.
+                   call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
+                   call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector)
+                   do kk=1,grd%itotsub
+                      i=grd%ltosi_s(kk)
+                      j=grd%ltosj_s(kk)
+                      work(kk)=grid2(i,j,1)
+                   enddo
+                else
+                   grid=rwork3d0(:,:,1)
+                   call general_fill_ns(grd,grid,work)
+                endif
+             endif
+             if ( icount == icm .or. k==nlevs ) then
+                call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
+                     g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
+             endif
+          enddo ! do k=1,nlevs
 
-      if (mype==mype_use(icount)) then
-         call read_vardata(filges, 'nconrd', rwork3d0, nslice=kr, slicedim=3)
-         ! rain number concentration.
-         if ( diff_res ) then
-            grid_b=rwork3d0(:,:,1)
-            vector(1)=.false.
-            call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
-            call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector)
-            do kk=1,grd%itotsub
-               i=grd%ltosi_s(kk)
-               j=grd%ltosj_s(kk)
-               work(kk)=grid2(i,j,1)
-            enddo
-         else
-            grid=rwork3d0(:,:,1)
-            call general_fill_ns(grd,grid,work)
-         endif
-      endif
-      if ( icount == icm .or. k==nlevs) then
-         call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
-              g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
-          endif
-       enddo ! do k=1,nlevs
+          do k=1,nlevs
+             icount=icount+1
+             iflag(icount)=16
+             ilev(icount)=k
+             kr = levs+1-k ! netcdf is top to bottom, need to flip
 
-!   do k=1,nlevs
-!      icount=icount+1
-!      iflag(icount)=17
-!      ilev(icount)=k
-!      kr = levs+1-k ! netcdf is top to bottom, need to flip
-!
-!      if (mype==mype_use(icount)) then
-!         call read_vardata(filges, 'cld_amt', rwork3d0, nslice=kr, slicedim=3)
-!         ! Cloud amount (cloud fraction). 
-!         if ( diff_res ) then
-!            grid_b=rwork3d0(:,:,1)
-!            vector(1)=.false.
-!            call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
-!            call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector)
-!            do kk=1,grd%itotsub
-!               i=grd%ltosi_s(kk)
-!               j=grd%ltosj_s(kk)
-!               work(kk)=grid2(i,j,1)
-!            enddo
-!         else
-!            grid=rwork3d0(:,:,1)
-!            call general_fill_ns(grd,grid,work)
-!         endif
+             if (mype==mype_use(icount)) then
+                call read_vardata(filges, 'nconrd', rwork3d0, nslice=kr, slicedim=3)
+                ! rain number concentration.
+                if ( diff_res ) then
+                   grid_b=rwork3d0(:,:,1)
+                   vector(1)=.false.
+                   call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
+                   call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector)
+                   do kk=1,grd%itotsub
+                      i=grd%ltosi_s(kk)
+                      j=grd%ltosj_s(kk)
+                      work(kk)=grid2(i,j,1)
+                   enddo
+                else
+                   grid=rwork3d0(:,:,1)
+                   call general_fill_ns(grd,grid,work)
+                endif
+             endif
+             if ( icount == icm .or. k==nlevs) then
+                call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
+                     g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
+                endif
+          enddo ! do k=1,nlevs
+       endif ! imp_physics
+!       do k=1,nlevs
+!          icount=icount+1
+!          iflag(icount)=17
+!          ilev(icount)=k
+!          kr = levs+1-k ! netcdf is top to bottom, need to flip
 !
-!      endif
-!      if ( icount == icm .or. k==nlevs ) then
-!         call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
-!              g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_cf)
-!      endif
-!   enddo ! do k=1,nlevs
+!          if (mype==mype_use(icount)) then
+!             call read_vardata(filges, 'cld_amt', rwork3d0, nslice=kr, slicedim=3)
+!             ! Cloud amount (cloud fraction). 
+!             if ( diff_res ) then
+!                grid_b=rwork3d0(:,:,1)
+!                vector(1)=.false.
+!                call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
+!                call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector)
+!                do kk=1,grd%itotsub
+!                   i=grd%ltosi_s(kk)
+!                   j=grd%ltosj_s(kk)
+!                   work(kk)=grid2(i,j,1)
+!                enddo
+!             else
+!                grid=rwork3d0(:,:,1)
+!                call general_fill_ns(grd,grid,work)
+!             endif
+!          endif
+!          if ( icount == icm .or. k==nlevs ) then
+!             call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
+!                  g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_cf)
+!          endif
+!       enddo ! do k=1,nlevs
 
    else ! read_2m