Transition to Thompson Microphysics for Microwave All-sky Assimilation (

#743)

src/enkf/gridio_gfs.f90

@@ -195,9 +195,6 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, &
   sst_ind = getindex(vars2d, 'sst')
   use_full_hydro = ( ql_ind > 0 .and. qi_ind > 0  .and. &
                      qr_ind > 0 .and. qs_ind > 0 .and. qg_ind > 0 )
-  ! Currently, we do not let precipiation to affect the enkf analysis  
-  ! The following line will be removed after testing
-  use_full_hydro = .false.
 
   if (.not. isinitialized) call init_spec_vars(nlons,nlats,ntrunc,4)
 
@@ -459,9 +456,6 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, &
   end if
 
   ! cloud derivatives
-  ! Currently, we do not let precipiation to affect the enkf analysis  
-  ! The following line will be removed after testing
-  use_full_hydro = .true.
   if (.not. use_full_hydro .and. iope==0) then
      if (ql_ind > 0 .or. qi_ind > 0) then
         do k=1,nlevs

src/gsi/general_read_gfsatm.f90

@@ -189,13 +189,14 @@ subroutine general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr,
 
 end subroutine general_reload
 subroutine 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,vdflag,g_cf)
+           g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vdflag,g_ni,g_nr,g_cf)
 
 ! !USES:
 
   use kinds, only: r_kind,i_kind
   use mpimod, only: npe,mpi_comm_world,ierror,mpi_rtype
   use general_sub2grid_mod, only: sub2grid_info
+  use ncepnems_io, only: imp_physics
   implicit none
 
 ! !INPUT PARAMETERS:
@@ -211,7 +212,8 @@ subroutine general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
   real(r_kind),dimension(grd%lat2,grd%lon2),         intent(  out) :: g_ps
   real(r_kind),dimension(grd%lat2,grd%lon2),         intent(inout) :: g_z
   real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig),intent(  out) :: g_u,g_v,&
-       g_vor,g_div,g_q,g_oz,g_tv,g_ql,g_qi,g_qr,g_qs,g_qg
+       g_vor,g_div,g_q,g_oz,g_tv,g_ql,g_qi,g_qr,g_qs,g_qg,g_ni,g_nr
+
   real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig),intent(  out),optional :: g_cf
 
 
@@ -392,7 +394,25 @@ subroutine general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
                g_qg(i,j,klev)=sub(ij,k)
             enddo
          enddo
-      elseif ( iflag(k) == 15 .and. present(g_cf) ) then  
+      elseif ( iflag(k) == 15 .and. imp_physics == 8) then
+         klev=ilev(k)
+         ij=0
+         do j=1,grd%lon2
+            do i=1,grd%lat2
+               ij=ij+1
+               g_ni(i,j,klev)=sub(ij,k)
+            enddo
+         enddo
+      elseif ( iflag(k) == 16 .and. imp_physics == 8) then
+         klev=ilev(k)
+         ij=0
+         do j=1,grd%lon2
+            do i=1,grd%lat2
+               ij=ij+1
+               g_nr(i,j,klev)=sub(ij,k)
+            enddo
+         enddo
+      elseif ( iflag(k) == 17 .and. present(g_cf) ) then 
          klev=ilev(k)
          ij=0
          do j=1,grd%lon2
@@ -2828,7 +2848,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
    real(r_kind),pointer,dimension(:,:)       :: g_t2m, g_q2m
    real(r_kind),pointer,dimension(:,:,:)     :: g_vor,g_div,&
                                                 g_q,g_oz,g_tv
-   real(r_kind),pointer,dimension(:,:,:)     :: g_ql,g_qi,g_qr,g_qs,g_qg
+   real(r_kind),pointer,dimension(:,:,:)     :: g_ql,g_qi,g_qr,g_qs,g_qg,g_ni,g_nr
 
    real(r_kind),allocatable,dimension(:,:)   :: g_z
    real(r_kind),allocatable,dimension(:,:,:) :: g_u,g_v
@@ -3038,6 +3058,8 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
        call gsi_bundlegetpointer(gfs_bundle,'qs',g_qs  ,ier);istatus1=istatus1+ier
        call gsi_bundlegetpointer(gfs_bundle,'qg',g_qg  ,ier);istatus1=istatus1+ier
     !  call gsi_bundlegetpointer(gfs_bundle,'cf',g_cf  ,ier);istatus1=istatus1+ier
+       call gsi_bundlegetpointer(gfs_bundle,'ni',g_ni  ,ier);istatus1=istatus1+ier
+       call gsi_bundlegetpointer(gfs_bundle,'nr',g_nr  ,ier);istatus1=istatus1+ier
        if ( istatus1 /= 0 ) then
           if ( mype == 0 ) then
              write(6,*) 'general_read_gfsatm_allhydro_nc: ERROR'
@@ -3100,7 +3122,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
       endif
       if ( icount == icm ) 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_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
       endif
    endif
 
@@ -3135,7 +3157,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
        endif
        if ( icount == icm ) 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_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
        endif
 
        !   Thermodynamic variable:  s-->g transform, communicate to all tasks
@@ -3172,7 +3194,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
           endif
           if ( icount == icm ) 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_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
           endif
        end do
 
@@ -3229,7 +3251,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
              endif
              if ( icount == icm ) 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_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
              endif
           end do
           do k=1,nlevs
@@ -3285,7 +3307,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
              endif
              if ( icount == icm ) 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_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
              endif
 
           end do
@@ -3321,7 +3343,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
              endif
              if ( icount == icm ) 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_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
              endif
 
              icount=icount+1
@@ -3352,7 +3374,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
              endif
              if ( icount == icm ) 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_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
              endif
           end do
        endif ! if ( uvflag )
@@ -3383,7 +3405,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
           endif
           if ( icount == icm ) 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_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
           endif
        end do
 
@@ -3414,7 +3436,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
           endif
           if ( icount == icm ) 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_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
           endif
        end do
 
@@ -3445,7 +3467,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
 
           if ( icount == icm ) 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_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
 
@@ -3475,7 +3497,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
           endif
           if ( icount == icm ) 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_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
 
@@ -3505,7 +3527,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
           endif
           if ( icount == icm ) 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_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
 
@@ -3535,7 +3557,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
           endif
           if ( icount == icm ) 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_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
 
@@ -3565,13 +3587,73 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
           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_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
+
+      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
+
+   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, '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)=15
+!      iflag(icount)=17
 !      ilev(icount)=k
 !      kr = levs+1-k ! netcdf is top to bottom, need to flip
 !

src/gsi/netcdfgfs_io.f90

@@ -159,17 +159,17 @@ subroutine read_
     real(r_kind),pointer,dimension(:,:,:):: ges_qs_it  => NULL()
     real(r_kind),pointer,dimension(:,:,:):: ges_qg_it  => NULL()
     real(r_kind),pointer,dimension(:,:,:):: ges_cf_it  => NULL()
-
+    real(r_kind),pointer,dimension(:,:,:):: ges_ni_it  => NULL()
+    real(r_kind),pointer,dimension(:,:,:):: ges_nr_it  => NULL()
     type(sub2grid_info) :: grd_t
     logical regional
     logical:: l_cld_derived,zflag,inithead
-
     type(gsi_bundle) :: atm_bundle
     type(gsi_grid)   :: atm_grid
     integer(i_kind),parameter :: n2d=2
     ! integer(i_kind),parameter :: n3d=8
     integer(i_kind),parameter :: n2d_2m=4
-    integer(i_kind),parameter :: n3d=14
+    integer(i_kind),parameter :: n3d=16
     character(len=4), parameter :: vars2d(n2d) = (/ 'z   ', 'ps  ' /)
     character(len=4), parameter :: vars2d_with2m(n2d_2m) = (/ 'z   ', 'ps  ','t2m ','q2m ' /)
     ! character(len=4), parameter :: vars3d(n3d) = (/ 'u   ', 'v   ', &
@@ -182,13 +182,16 @@ subroutine read_
                                                     'cw  ', 'oz  ', &
                                                     'ql  ', 'qi  ', &
                                                     'qr  ', 'qs  ', &
-                                                    'qg  ', 'cf  ' /)
+                                                    'qg  ', 'ni  ', &
+                                                    'nr  ', 'cf  ' /)
+
     real(r_kind),pointer,dimension(:,:):: ptr2d   =>NULL()
     real(r_kind),pointer,dimension(:,:,:):: ptr3d =>NULL()
 
     regional=.false.
     inner_vars=1
-    num_fields=min(14*grd_a%nsig+2,npe)
+
+    num_fields=min(n3d*grd_a%nsig+2,npe)
 !  Create temporary communication information fore read routines
     call general_sub2grid_create_info(grd_t,inner_vars,grd_a%nlat,grd_a%nlon, &
           grd_a%nsig,num_fields,regional)
@@ -200,6 +203,7 @@ subroutine read_
     else
         call gsi_bundlecreate(atm_bundle,atm_grid,'aux-atm-read',istatus,names2d=vars2d,names3d=vars3d)
     endif
+
     if(istatus/=0) then
       write(6,*) myname_,': trouble creating atm_bundle'
       call stop2(999)
@@ -248,6 +252,8 @@ subroutine read_
                    if (associated(ges_qr_it)) ges_qr_it(i,j,k) = max(qcmin,ges_qr_it(i,j,k))
                    if (associated(ges_qs_it)) ges_qs_it(i,j,k) = max(qcmin,ges_qs_it(i,j,k))
                    if (associated(ges_qg_it)) ges_qg_it(i,j,k) = max(qcmin,ges_qg_it(i,j,k))
+                   if (associated(ges_ni_it)) ges_ni_it(i,j,k) = max(qcmin,ges_ni_it(i,j,k))
+                   if (associated(ges_nr_it)) ges_nr_it(i,j,k) = max(qcmin,ges_nr_it(i,j,k))
                    if (associated(ges_cf_it)) ges_cf_it(i,j,k) = min(max(zero,ges_cf_it(i,j,k)),one)
                 enddo
              enddo
@@ -256,6 +262,8 @@ subroutine read_
           l_cld_derived = associated(ges_cwmr_it).and.&
                           associated(ges_q_it)   .and.&
                           associated(ges_ql_it)  .and.&
+                          associated(ges_ni_it)  .and.&
+                          associated(ges_nr_it)  .and.&
                           associated(ges_qi_it)  .and.& 
                           associated(ges_tv_it)
 !         call set_cloud_lower_bound(ges_cwmr_it)
@@ -343,6 +351,16 @@ subroutine set_guess_
        call gsi_bundlegetpointer (gsi_metguess_bundle(it),'ql',ges_ql_it,istatus)
        if(istatus==0) ges_ql_it = ptr3d
     endif
+    call gsi_bundlegetpointer (atm_bundle,'ni',ptr3d,istatus)
+    if (istatus==0) then
+       call gsi_bundlegetpointer (gsi_metguess_bundle(it),'ni',ges_ni_it,istatus)
+       if(istatus==0) ges_ni_it = ptr3d
+    endif
+    call gsi_bundlegetpointer (atm_bundle,'nr',ptr3d,istatus)
+    if (istatus==0) then
+       call gsi_bundlegetpointer (gsi_metguess_bundle(it),'nr',ges_nr_it,istatus)
+       if(istatus==0) ges_nr_it = ptr3d
+    endif
     call gsi_bundlegetpointer (atm_bundle,'qi',ptr3d,istatus)
     if (istatus==0) then
        call gsi_bundlegetpointer (gsi_metguess_bundle(it),'qi',ges_qi_it,istatus)