Transition to Thompson Microphysics for Microwave All-sky Assimilation (

NOAA-EMC#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)