From f9519080433ea3c6bb4007bf3898f852b81a2d6a Mon Sep 17 00:00:00 2001 From: "Guoqing.Ge" Date: Thu, 22 Sep 2022 17:35:09 -0600 Subject: [PATCH 1/2] update for 2mT DA --- src/gsi/general_read_gfsatm.f90 | 338 ++++++++++++++++++++++++-------- src/gsi/gsimod.F90 | 4 +- src/gsi/jfunc.f90 | 5 + src/gsi/netcdfgfs_io.f90 | 42 +++- src/gsi/read_prepbufr.f90 | 23 ++- src/gsi/setupt.f90 | 169 +++++++++++++--- src/gsi/update_guess.f90 | 6 +- 7 files changed, 461 insertions(+), 126 deletions(-) diff --git a/src/gsi/general_read_gfsatm.f90 b/src/gsi/general_read_gfsatm.f90 index 39db75db73..790b2eb321 100755 --- a/src/gsi/general_read_gfsatm.f90 +++ b/src/gsi/general_read_gfsatm.f90 @@ -188,6 +188,72 @@ subroutine general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, return end subroutine general_reload + +! 2m reload +subroutine general_2m_reload(grd,g_t2m, g_q2m,g_ps,icount,iflag,work) +! !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 + + implicit none +! !INPUT PARAMETERS: + + type(sub2grid_info), intent(in ) :: grd + integer(i_kind), intent(inout) :: icount + integer(i_kind),dimension(npe), intent(inout) :: iflag + real(r_kind),dimension(grd%itotsub),intent(in ) :: work + +! !OUTPUT PARAMETERS: + + real(r_kind),dimension(grd%lat2,grd%lon2), intent( out) :: g_t2m, g_q2m, g_ps + +! !DESCRIPTION: 2m version og general_reload. +!------------------------------------------------------------------------- + + integer(i_kind) i,j,ij,k + real(r_kind),dimension(grd%lat2*grd%lon2,npe):: sub + + call mpi_alltoallv(work,grd%sendcounts_s,grd%sdispls_s,mpi_rtype,& + sub,grd%recvcounts_s,grd%rdispls_s,mpi_rtype,& + mpi_comm_world,ierror) + +!$omp parallel do schedule(dynamic,1) private(k,i,j,ij) + do k=1,icount + if ( iflag(k) == 1 ) then + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_t2m(i,j)=sub(ij,k) + enddo + enddo + elseif ( iflag(k) == 2 ) then + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_q2m(i,j)=sub(ij,k) + enddo + enddo + elseif ( iflag(k) == 3 ) then + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_ps(i,j)=sub(ij,k) + enddo + enddo + endif + enddo ! do k=1,icount + + icount=0 + iflag=0 + + return + +end subroutine general_2m_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) @@ -1892,7 +1958,7 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & use gsi_bundlemod, only: gsi_bundlegetpointer use module_ncio, only: Dataset, Variable, Dimension, open_dataset,& close_dataset, get_dim, read_vardata,get_idate_from_time_units - use gfsreadmod, only: general_reload + use gfsreadmod, only: general_reload, general_2m_reload implicit none @@ -1903,13 +1969,14 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & type(sub2grid_info) ,intent(in ) :: grd type(spec_vars) ,intent(in ) :: sp_a character(*) ,intent(in ) :: filename - logical ,intent(in ) :: uvflag,zflag,vordivflag,init_head + logical ,intent(in ) :: uvflag,zflag,vordivflag,init_head integer(i_kind) ,intent( out) :: iret_read type(gsi_bundle) ,intent(inout) :: gfs_bundle real(r_kind),pointer,dimension(:,:) :: ptr2d real(r_kind),pointer,dimension(:,:,:) :: ptr3d real(r_kind),pointer,dimension(:,:) :: g_ps + real(r_kind),pointer,dimension(:,:) :: g_t2m, g_q2m real(r_kind),pointer,dimension(:,:,:) :: g_vor,g_div,& g_cwmr,g_q,g_oz,g_tv @@ -1944,8 +2011,7 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & logical,dimension(1) :: vector type(Dataset) :: atmges type(Dimension) :: ncdim - - + logical :: read_2m, read_z !****************************************************************************** ! Initialize variables used below @@ -1959,11 +2025,26 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & mype_use=-1 icount=0 procuse=.false. + + if (filename(1:3) == 'sfc') then + read_2m = .true. + read_z = .false. + if ( mype == 0 ) write(6,* ) & + trim(my_name), ': reading 2m variables from ', trim(filename) + else + read_2m = .false. + read_z = zflag + if ( mype == 0 ) write(6,* ) & + trim(my_name), ': reading atmos variables from ', trim(filename) + endif + + if ( mype == 0 ) procuse = .true. do i=1,npe if ( grd%recvcounts_s(i-1) > 0 ) then icount = icount+1 mype_use(icount)=i-1 + ! mype_use lists the pes being used, procuse true if this proc being used if ( i-1 == mype ) procuse=.true. endif enddo @@ -1997,6 +2078,7 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & ! get dimension sizes ncdim = get_dim(atmges, 'grid_xt'); lonb = ncdim%len ncdim = get_dim(atmges, 'grid_yt'); latb = ncdim%len + if (.not. read_2m) & ncdim = get_dim(atmges, 'pfull'); levs = ncdim%len ! get time information @@ -2030,12 +2112,14 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & trim(my_name),grd%nlon,lonb !call stop2(101) endif + if (.not. read_2m) then if ( levs /= grd%nsig ) then if ( mype == 0 ) write(6, & '(a,'': inconsistent spatial dimension nsig = '',i4,tr1,''levs = '',i4)') & trim(my_name),grd%nsig,levs call stop2(101) endif + endif allocate( spec_vor(sp_a%nc), spec_div(sp_a%nc) ) allocate( grid(grd%nlon,nlatm2), grid_v(grd%nlon,nlatm2) ) @@ -2073,57 +2157,73 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & endif ! if ( procuse ) ! Get pointer to relevant variables (this should be made flexible and general) - iredundant=0 - call gsi_bundlegetpointer(gfs_bundle,'sf',g_div ,ier) - if ( ier == 0 ) iredundant = iredundant + 1 - call gsi_bundlegetpointer(gfs_bundle,'div',g_div ,ier) - if ( ier == 0 ) iredundant = iredundant + 1 - if ( iredundant==2 ) then - if ( mype == 0 ) then - write(6,*) 'general_read_gfsatm_nems: ERROR' - write(6,*) 'cannot handle having both sf and div' - write(6,*) 'Aborting ... ' - endif - call stop2(999) - endif - iredundant=0 - call gsi_bundlegetpointer(gfs_bundle,'vp',g_vor ,ier) - if ( ier == 0 ) iredundant = iredundant + 1 - call gsi_bundlegetpointer(gfs_bundle,'vor',g_vor ,ier) - if ( ier == 0 ) iredundant = iredundant + 1 - if ( iredundant==2 ) then - if ( mype == 0 ) then - write(6,*) 'general_read_gfsatm_nems: ERROR' - write(6,*) 'cannot handle having both vp and vor' - write(6,*) 'Aborting ... ' - endif - call stop2(999) - endif - iredundant=0 - call gsi_bundlegetpointer(gfs_bundle,'t' ,g_tv ,ier) - if ( ier == 0 ) iredundant = iredundant + 1 - call gsi_bundlegetpointer(gfs_bundle,'tv',g_tv ,ier) - if ( ier == 0 ) iredundant = iredundant + 1 - if ( iredundant==2 ) then - if ( mype == 0 ) then - write(6,*) 'general_read_gfsatm_nems: ERROR' - write(6,*) 'cannot handle having both t and tv' - write(6,*) 'Aborting ... ' - endif - call stop2(999) - endif - istatus=0 - call gsi_bundlegetpointer(gfs_bundle,'ps',g_ps ,ier);istatus=istatus+ier - call gsi_bundlegetpointer(gfs_bundle,'q' ,g_q ,ier);istatus=istatus+ier - call gsi_bundlegetpointer(gfs_bundle,'oz',g_oz ,ier);istatus=istatus+ier - call gsi_bundlegetpointer(gfs_bundle,'cw',g_cwmr,ier);istatus=istatus+ier - if ( istatus /= 0 ) then - if ( mype == 0 ) then - write(6,*) 'general_read_gfsatm_nems: ERROR' - write(6,*) 'Missing some of the required fields' - write(6,*) 'Aborting ... ' - endif - call stop2(999) + if (.not. read_2m) then + iredundant=0 + call gsi_bundlegetpointer(gfs_bundle,'sf',g_div ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + call gsi_bundlegetpointer(gfs_bundle,'div',g_div ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + if ( iredundant==2 ) then + if ( mype == 0 ) then + write(6,*) 'general_read_gfsatm_nems: ERROR' + write(6,*) 'cannot handle having both sf and div' + write(6,*) 'Aborting ... ' + endif + call stop2(999) + endif + iredundant=0 + call gsi_bundlegetpointer(gfs_bundle,'vp',g_vor ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + call gsi_bundlegetpointer(gfs_bundle,'vor',g_vor ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + if ( iredundant==2 ) then + if ( mype == 0 ) then + write(6,*) 'general_read_gfsatm_nems: ERROR' + write(6,*) 'cannot handle having both vp and vor' + write(6,*) 'Aborting ... ' + endif + call stop2(999) + endif + iredundant=0 + call gsi_bundlegetpointer(gfs_bundle,'t' ,g_tv ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + call gsi_bundlegetpointer(gfs_bundle,'tv',g_tv ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + if ( iredundant==2 ) then + if ( mype == 0 ) then + write(6,*) 'general_read_gfsatm_nems: ERROR' + write(6,*) 'cannot handle having both t and tv' + write(6,*) 'Aborting ... ' + endif + call stop2(999) + endif + + istatus=0 + call gsi_bundlegetpointer(gfs_bundle,'ps',g_ps ,ier);istatus=istatus+ier + call gsi_bundlegetpointer(gfs_bundle,'q' ,g_q ,ier);istatus=istatus+ier + call gsi_bundlegetpointer(gfs_bundle,'oz',g_oz ,ier);istatus=istatus+ier + call gsi_bundlegetpointer(gfs_bundle,'cw',g_cwmr,ier);istatus=istatus+ier + if ( istatus /= 0 ) then + if ( mype == 0 ) then + write(6,*) 'general_read_gfsatm_nems: ERROR' + write(6,*) 'Missing some of the required fields' + write(6,*) 'Aborting ... ' + endif + call stop2(999) + endif + else ! read 2m vars + istatus=0 + call gsi_bundlegetpointer(gfs_bundle,'t2m',g_t2m ,ier);istatus=istatus+ier + call gsi_bundlegetpointer(gfs_bundle,'q2m',g_q2m ,ier);istatus=istatus+ier + call gsi_bundlegetpointer(gfs_bundle,'ps',g_ps ,ier);istatus=istatus+ier + if ( istatus /= 0 ) then + if ( mype == 0 ) then + write(6,*) 'general_read_gfsatm_nems: ERROR' + write(6,*) 'Missing 2m required variabless' + write(6,*) 'Aborting ... ' + endif + call stop2(999) + endif endif allocate(g_u(grd%lat2,grd%lon2,grd%nsig),g_v(grd%lat2,grd%lon2,grd%nsig)) allocate(g_z(grd%lat2,grd%lon2)) @@ -2135,8 +2235,8 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & ! Once on the grid, fields need to be scattered from the full domain to ! sub-domains. - ! Only read Terrain when zflag is true. - if ( zflag ) then + ! Only read Terrain when read_z is true. + if ( read_z ) then icount=icount+1 iflag(icount)=1 @@ -2161,12 +2261,14 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & call general_fill_ns(grd,grid,work) endif endif - if ( icount == icm ) then + if ( read_2m .or. ( (.not. read_2m ) .and. ( icount == icm)) ) then call general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, & icount,iflag,ilev,work,uvflag,vordivflag) endif endif + if (.not. read_2m) then + icount=icount+1 iflag(icount)=2 ilev(icount)=1 @@ -2508,6 +2610,82 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & endif enddo ! do k=1,nlevs + elseif (read_2m) then ! read_2m + icount=1 + iflag(icount)=1 + + ! 2m temperature from sfc file + if (mype==mype_use(icount)) then + ! read t2m + call read_vardata(atmges, 'tmp2m', rwork2d) + + if ( diff_res ) then + vector(1)=.false. + grid_b=rwork2d + 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=rwork2d + call general_fill_ns(grd,grid,work) + endif + endif + !if ( icount == icm ) then + ! call general_2m_reload(grd,g_t2m, g_q2m, icount,iflag,work) + !endif + + icount=icount + 1 + iflag(icount)=2 + + ! 2m temperature from sfc file + if (mype==mype_use(icount)) then + ! read ps + call read_vardata(atmges, 'spfh2m', rwork2d) + if ( diff_res ) then + vector(1)=.false. + grid_b=rwork2d + 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=rwork2d + call general_fill_ns(grd,grid,work) + endif + endif + if (mype==mype_use(icount)) then + ! read ps + call read_vardata(atmges, 'pressfc', rwork2d) + rwork2d = r0_001*rwork2d ! convert Pa to cb + if ( diff_res ) then + vector(1)=.false. + grid_b=rwork2d + 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=rwork2d + call general_fill_ns(grd,grid,work) + endif + endif + + ! not using all procs. doesn't trigger. todo: fix task allocation + !if ( icount == icm ) then + call general_2m_reload(grd,g_t2m, g_q2m, g_ps, icount,iflag,work) + !endif + + endif ! read_2m if ( procuse ) then if ( diff_res) deallocate(grid_b,grid_b2,grid_c,grid_c2,grid2) @@ -2530,27 +2708,29 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & !enddo ! Load u->div and v->vor slot when uv are used instead - if ( uvflag ) then - call gsi_bundlegetpointer(gfs_bundle,'u' ,ptr3d,ier) - if ( ier == 0 ) then - ptr3d=g_u - call gsi_bundlegetpointer(gfs_bundle,'v' ,ptr3d,ier) - if ( ier == 0 ) ptr3d=g_v - else ! in this case, overload: return u/v in sf/vp slot - call gsi_bundlegetpointer(gfs_bundle,'sf' ,ptr3d,ier) - if ( ier == 0 ) then - ptr3d=g_u - call gsi_bundlegetpointer(gfs_bundle,'vp' ,ptr3d,ier) - if ( ier == 0 ) ptr3d=g_v - endif - endif - else ! in this case, overload: return u/v in sf/vp slot - call gsi_bundlegetpointer(gfs_bundle,'sf' ,ptr3d,ier) - if ( ier == 0 ) ptr3d=g_u - call gsi_bundlegetpointer(gfs_bundle,'vp' ,ptr3d,ier) - if ( ier == 0 ) ptr3d=g_v - endif - if (zflag) then + if ( .not. read_2m ) then + if ( uvflag ) then + call gsi_bundlegetpointer(gfs_bundle,'u' ,ptr3d,ier) + if ( ier == 0 ) then + ptr3d=g_u + call gsi_bundlegetpointer(gfs_bundle,'v' ,ptr3d,ier) + if ( ier == 0 ) ptr3d=g_v + else ! in this case, overload: return u/v in sf/vp slot + call gsi_bundlegetpointer(gfs_bundle,'sf' ,ptr3d,ier) + if ( ier == 0 ) then + ptr3d=g_u + call gsi_bundlegetpointer(gfs_bundle,'vp' ,ptr3d,ier) + if ( ier == 0 ) ptr3d=g_v + endif + endif + else ! in this case, overload: return u/v in sf/vp slot + call gsi_bundlegetpointer(gfs_bundle,'sf' ,ptr3d,ier) + if ( ier == 0 ) ptr3d=g_u + call gsi_bundlegetpointer(gfs_bundle,'vp' ,ptr3d,ier) + if ( ier == 0 ) ptr3d=g_v + endif + endif ! read_2m + if (read_z) then call gsi_bundlegetpointer(gfs_bundle,'z' ,ptr2d,ier) if ( ier == 0 ) ptr2d=g_z endif diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index ff636bd70f..b5019c7dfa 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -99,7 +99,7 @@ module gsimod factv,factl,factp,factg,factw10m,facthowv,factcldch,niter,niter_no_qc,biascor,& init_jfunc,qoption,cwoption,switch_on_derivatives,tendsflag,jiterstart,jiterend,R_option,& bcoption,diurnalbc,print_diag_pcg,tsensible,diag_precon,step_start,pseudo_q2,& - clip_supersaturation,cnvw_option + clip_supersaturation,cnvw_option,hofx_2m_sfcfile use state_vectors, only: init_anasv,final_anasv use control_vectors, only: init_anacv,final_anacv,nrf,nvars,nrf_3d,cvars3d,cvars2d,& nrf_var,lcalc_gfdl_cfrac,incvars_to_zero,incvars_zero_strat,incvars_efold @@ -1042,7 +1042,7 @@ module gsimod ! l_foreaft_thin - separate TDR fore/aft scan for thinning namelist/obs_input/dmesh,time_window_max,time_window_rad, & - ext_sonde,l_foreaft_thin + ext_sonde,l_foreaft_thin,hofx_2m_sfcfile ! SINGLEOB_TEST (one observation test case setup): ! maginnov - magnitude of innovation for one ob diff --git a/src/gsi/jfunc.f90 b/src/gsi/jfunc.f90 index d65b4c8fd3..1f5924102e 100644 --- a/src/gsi/jfunc.f90 +++ b/src/gsi/jfunc.f90 @@ -136,10 +136,12 @@ module jfunc public :: pseudo_q2 public :: varq public :: cnvw_option + public :: hofx_2m_sfcfile logical first,last,switch_on_derivatives,tendsflag,print_diag_pcg,tsensible,diag_precon logical clip_supersaturation,R_option logical pseudo_q2,limitqobs + logical hofx_2m_sfcfile logical cnvw_option integer(i_kind) iout_iter,miter,iguess,nclen,qoption,cwoption integer(i_kind) jiter,jiterstart,jiterend,iter @@ -249,6 +251,9 @@ subroutine init_jfunc ! option for including convective clouds in the all-sky assimilation cnvw_option=.false. +! to calculate hofx for 2m obs, interpolate from gfs sfc file. + hofx_2m_sfcfile=.false. + return end subroutine init_jfunc diff --git a/src/gsi/netcdfgfs_io.f90 b/src/gsi/netcdfgfs_io.f90 index 41eb05c228..013cbcb4f1 100644 --- a/src/gsi/netcdfgfs_io.f90 +++ b/src/gsi/netcdfgfs_io.f90 @@ -130,6 +130,7 @@ subroutine read_ use mpimod, only: npe,mype use cloud_efr_mod, only: cloud_calc_gfs,set_cloud_lower_bound use gridmod, only: fv3_full_hydro + use jfunc, only: hofx_2m_sfcfile implicit none @@ -141,6 +142,8 @@ subroutine read_ real(r_kind),pointer,dimension(:,: ):: ges_ps_it =>NULL() real(r_kind),pointer,dimension(:,: ):: ges_z_it =>NULL() + real(r_kind),pointer,dimension(:,: ):: ges_t2m_it =>NULL() + real(r_kind),pointer,dimension(:,: ):: ges_q2m_it =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_u_it =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_v_it =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_div_it =>NULL() @@ -163,13 +166,11 @@ subroutine read_ type(gsi_bundle) :: atm_bundle type(gsi_grid) :: atm_grid integer(i_kind),parameter :: n2d=2 + integer(i_kind),parameter :: n2d_2m=4 ! integer(i_kind),parameter :: n3d=8 integer(i_kind),parameter :: n3d=14 character(len=4), parameter :: vars2d(n2d) = (/ 'z ', 'ps ' /) - ! character(len=4), parameter :: vars3d(n3d) = (/ 'u ', 'v ', & - ! 'vor ', 'div ', & - ! 'tv ', 'q ', & - ! 'cw ', 'oz ' /) + character(len=4), parameter :: vars2d_with2m(n2d_2m) = (/ 'z ', 'ps ','t2m ','q2m ' /) character(len=4), parameter :: vars3d(n3d) = (/ 'u ', 'v ', & 'vor ', 'div ', & 'tv ', 'q ', & @@ -189,8 +190,11 @@ subroutine read_ ! Allocate bundle used for reading members call gsi_gridcreate(atm_grid,lat2,lon2,nsig) - - call gsi_bundlecreate(atm_bundle,atm_grid,'aux-atm-read',istatus,names2d=vars2d,names3d=vars3d) + if (hofx_2m_sfcfile) then + call gsi_bundlecreate(atm_bundle,atm_grid,'aux-atm-read',istatus,names2d=vars2d_with2m,names3d=vars3d) + 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) @@ -198,16 +202,26 @@ subroutine read_ do it=1,nfldsig - write(filename,'(''sigf'',i2.2)') ifilesig(it) - ! Read background fields into bundle - if (fv3_full_hydro) then + if (hofx_2m_sfcfile) then ! current read_sfc routines called from different part of + ! of the code, can't easily read into the met-bundle + ! wrote a new routine here + write(filename,'(''sfcf'',i2.2)') ifilesig(it) + call general_read_gfsatm_nc(grd_t,sp_a,filename,.true.,.true.,.true.,& + atm_bundle,.true.,istatus) + + write(filename,'(''sigf'',i2.2)') ifilesig(it) + call general_read_gfsatm_nc(grd_t,sp_a,filename,.true.,.true.,.true.,& + atm_bundle,.true.,istatus) + elseif (fv3_full_hydro) then if (mype==0) write(*,*) 'calling general_read_gfsatm_allhydro_nc', it + write(filename,'(''sigf'',i2.2)') ifilesig(it) call general_read_gfsatm_allhydro_nc(grd_t,sp_a,filename,.true.,.true.,.true.,& atm_bundle,istatus) ! this loads cloud and precip if (mype==0) write(*,*) 'done with general_read_gfsatm_allhydro_nc', it else if (mype==0) write(*,*) 'calling general_read_gfsatm_nc' + write(filename,'(''sigf'',i2.2)') ifilesig(it) call general_read_gfsatm_nc(grd_t,sp_a,filename,.true.,.true.,.true.,& atm_bundle,.true.,istatus) if (mype==0) write(*,*) 'done with general_read_gfsatm_nc' @@ -273,6 +287,16 @@ subroutine set_guess_ call gsi_bundlegetpointer (gsi_metguess_bundle(it),'z' ,ges_z_it ,istatus) if(istatus==0) ges_z_it = ptr2d endif + call gsi_bundlegetpointer (atm_bundle,'t2m',ptr2d,istatus) + if (istatus==0) then + call gsi_bundlegetpointer (gsi_metguess_bundle(it),'t2m' ,ges_t2m_it ,istatus) + if(istatus==0) ges_t2m_it = ptr2d + endif + call gsi_bundlegetpointer (atm_bundle,'q2m',ptr2d,istatus) + if (istatus==0) then + call gsi_bundlegetpointer (gsi_metguess_bundle(it),'q2m' ,ges_q2m_it ,istatus) + if(istatus==0) ges_q2m_it = ptr2d + endif call gsi_bundlegetpointer (atm_bundle,'u',ptr3d,istatus) if (istatus==0) then call gsi_bundlegetpointer (gsi_metguess_bundle(it),'u' ,ges_u_it ,istatus) diff --git a/src/gsi/read_prepbufr.f90 b/src/gsi/read_prepbufr.f90 index d2cb503926..21d73b41ce 100644 --- a/src/gsi/read_prepbufr.f90 +++ b/src/gsi/read_prepbufr.f90 @@ -212,7 +212,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& use hilbertcurve,only: init_hilbertcurve, accum_hilbertcurve, & apply_hilbertcurve,destroy_hilbertcurve use ndfdgrids,only: init_ndfdgrid,destroy_ndfdgrid,relocsfcob,adjust_error - use jfunc, only: tsensible + use jfunc, only: tsensible,hofx_2m_sfcfile use deter_sfc_mod, only: deter_sfc_type,deter_sfc2 use gsi_nstcouplermod, only: nst_gsi,nstinfo use gsi_nstcouplermod, only: gsi_nstcoupler_deter @@ -1617,7 +1617,9 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& ! If temperature ob, extract information regarding virtual ! versus sensible temperature if(tob) then - if (.not. twodvar_regional .or. .not.tsensible) then + ! use tvirtual if tsensible flag not set, and not in either 2Dregional or global_2m DA mode + if ( (.not. tsensible) .and. .not. (twodvar_regional .or. hofx_2m_sfcfile) ) then + do k=1,levs tvflg(k)=one ! initialize as sensible do j=1,20 @@ -1914,6 +1916,12 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& ! Missing Values ==> Cycling! In this case for howv only. #ww3 if (howvob .and. owave(1,k) > r0_1_bmiss) cycle LOOP_K_LEVS +! CSD - temporary hack ( move to prepbufr pre-processing later) +! Over-ride QM=9 for sfc obs + if (sfctype .and. hofx_2m_sfcfile ) then + if (tob .and. qm == 9 ) qm = 2 ! 2=not checked + if (qob .and. qm == 9 ) qm = 2 ! 2=not checked + endif ! Set usage variable usage = zero if(icuse(nc) <= 0)usage=100._r_kind @@ -2068,10 +2076,13 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& ! Temperature if(tob) then ppb=obsdat(1,k) + errout = one if (aircraftobs .and. aircraft_t_bc .and. acft_profl_file) then call errormod_aircraft(pqm,tqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm,hdr3) else - call errormod(pqm,tqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm) + if (.not. (sfctype .and. hofx_2m_sfcfile)) then + call errormod(pqm,tqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm) + endif end if toe=obserr(3,k)*errout qtflg=tvflg(k) @@ -2326,7 +2337,11 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& end if qobcon=obsdat(2,k)*convert tdry=r999 - if (tqm(k)0) hofx_2m_sfcfile=.false. + thead => obsLL(:) save_jacobian = conv_diagsave .and. jiter==jiterstart .and. lobsdiag_forenkf @@ -432,8 +444,16 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav hr_offset=min_offset/60.0_r_kind dup=one do k=1,nobs + ikx=nint(data(ikxx,k)) + itype=ictype(ikx) + sfctype=(itype>179.and.itype<190).or.(itype>=192.and.itype<=199) + ! landsfctype below is used to restrict hofx_2m_sfcfile to land obs only. + ! GDAS assmilates 180 and 182 over ocean. Would probably be better to read h(x) from + ! the surface file for these, but have left as is (default is LML) for zero-diff + ! changes. + landsfctype =( itype==181 .or. itype==183 .or. itype==187 ) do l=k+1,nobs - if (twodvar_regional) then + if (twodvar_regional .or. (hofx_2m_sfcfile .and. sfctype) ) then duplogic=data(ilat,k) == data(ilat,l) .and. & data(ilon,k) == data(ilon,l) .and. & data(ier,k) < r1000 .and. data(ier,l) < r1000 .and. & @@ -464,8 +484,16 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav end do end do -! Run a buddy-check - if (twodvar_regional .and. buddycheck_t) call buddy_check_t(is,data,luse,mype,nele,nobs,muse,buddyuse) +! Run a buddy-check (does not work for hofx_2m_sfcfile, +! pre-existing routine crashes with an MPI problem ) +! Note: current params have buddy radius of 108 km, max diff of 8 K. +! gross error check removes O-F > 7., so this is probably removing +! most obs that fail the buddy check already +! GSD uses the buddy check to expand gross error check for obs that pass +! the buddy check. Not sure if we want to do this globally. Turn off for now. + !if ( (twodvar_regional .or. hofx_2m_sfcfile) .and. buddycheck_t) & + if ( (twodvar_regional) .and. buddycheck_t) & + call buddy_check_t(is,data,luse,mype,nele,nobs,muse,buddyuse) ! If requested, save select data for output to diagnostic file if(conv_diagsave)then @@ -654,17 +682,24 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav call tintrp2a1(ges_lnprsl,prsltmp,dlat,dlon,dtime,hrdifsig,& nsig,mype,nfldsig) - drpx=zero - if(sfctype .and. .not.twodvar_regional) then - drpx=abs(one-((one/exp(dpres-log(psges))))**rd_over_cp)*t0c - end if + if ( hofx_2m_sfcfile .and. landsfctype) then + drpx = zero + dpres = one ! put obs at surface + else + drpx=zero + if(sfctype .and. .not.twodvar_regional) then + drpx=abs(one-((one/exp(dpres-log(psges))))**rd_over_cp)*t0c + end if -! Put obs pressure in correct units to get grid coord. number - call grdcrd1(dpres,prsltmp(1),nsig,-1) + ! Put obs pressure in correct units to get grid coord. number + call grdcrd1(dpres,prsltmp(1),nsig,-1) + endif ! Implementation of forward model ---------- - if(sfctype.and.sfcmodel) then + if(sfctype .and. sfcmodel) then +! SCENARIO 1: If obs is sfctype, and sfcmodel is requested +! outdated (delete?) tgges=data(iskint,i) roges=data(isfcr,i) @@ -693,8 +728,48 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav prsltmp2(2), tvtmp(2), qtmp(2), hsges(1), roges, msges, & f10ges,u10ges,v10ges, t2ges, q2ges, regime, iqtflg) tges = t2ges - + elseif (landsfctype .and. hofx_2m_sfcfile ) then +! SCENARIO 2: obs is sfctype, and hofx_2m_sfcfile scheme is on. +! 2m forecast has been read from the sfc guess files + + ! mask: 0 - sea, 1 - land, 2-ice (val + 3 = interpolated from at least one different land cover + ! for now, use only pure land (interpolation from mixed grid cells is pretty sketchy, but + ! may result in lack of obs near coasts) + if (int(data(idomsfc,i)) .NE. 1 ) muse(i) = .false. + + call tintrp2a11(ges_t2m,tges2m,dlat,dlon,dtime,hrdifsig,& + mype,nfldsig) + +! correct obs to model terrain height - equivalent to gsd_terrain_match, which uses lapse rate +! calculated from model levels. For now, use standard lapse rate, then look +! at updating once the number of model levels has been finalised (and perhaps make dependent on the number of model +! levels) + + delta_z = data(izz,i) - data(istnelv,i) + tob = tob + delta_z*T_lapse + !update the station elevation + data(istnelv,i) = data(izz,i) + + if(save_jacobian) then + t_ind = getindex(svars2d, 't2m') + if (t_ind < 0) then + print *, 'Error: no variable t2m in state vector.Exiting.' + call stop2(1300) + endif + dhx_dx%st_ind(1) = sum(levels(1:ns3d)) + t_ind + dhx_dx%end_ind(1) = sum(levels(1:ns3d)) + t_ind + dhx_dx%val(1) = one + dhx_dx%val(2) = zero ! in this case, there is no vertical interp + ! and nnz (=dim(dhx_dx%val)) should be one, + ! but nnz is a file attribute, so need to use + ! same value as for vertical profile obs. Get + ! around this by setting val(2) to zero. + endif else +! SCENARIO 3: obs is sfctype, and neither sfcmodel nor hofx_2m_sfcfile is chosen +! .or. obs is not sfctype +! ! SCENARIO 3a: obs is a virtual temp. + write(6,*) ' CSD going into iqtflg', iqtflg, dpres if(iqtflg)then ! Interpolate guess tv to observation location and time call tintrp31(ges_tv,tges,dlat,dlon,dpres,dtime, & @@ -716,6 +791,8 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav dhx_dx%val(1) = one - delz ! weight for iz's level dhx_dx%val(2) = delz ! weight for iz+1's level endif +! ! SCENARIO 3b: obs is a sensibl temp. +! ! this was the default for T2m else ! Interpolate guess tsen to observation location and time call tintrp31(ges_tsen,tges,dlat,dlon,dpres,dtime, & @@ -738,7 +815,10 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav dhx_dx%val(2) = delz ! weight for iz+1's level endif end if + write(6,*) ' CSD coming out iqtflg', tges + +! SCENARIO 4: obs is sfctype, and i_use_2mt4b flag is on (turns on LAM sfc DA) if(i_use_2mt4b>0 .and. sfctype) then if(i_coastline==1 .or. i_coastline==3) then @@ -773,17 +853,23 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav call grdcrd1(sfcchk,prsltmp(1),nsig,-1) ! Check to see if observations is above the top of the model (regional mode) - if(sfctype)then + if(sfctype .and. .not. hofx_2m_sfcfile )then if(abs(dpres)>four) drpx=1.0e10_r_kind pres_diff=prest-r10*psges if (twodvar_regional .and. abs(pres_diff)>=r1000) drpx=1.0e10_r_kind end if - rlow=max(sfcchk-dpres,zero) -! linear variation of observation ramp [between grid points 1(~3mb) and 15(~45mb) below the surface] - if(l_sfcobserror_ramp_t) then - ramp=min(max(((rlow-1.0_r_kind)/(15.0_r_kind-1.0_r_kind)),0.0_r_kind),1.0_r_kind) + + if (.not. (hofx_2m_sfcfile .and. sfctype) ) then + rlow=max(sfcchk-dpres,zero) ! sfcchk = k of sfc, dpres k of obs + ! linear variation of observation ramp [between grid points 1(~3mb) and 15(~45mb) below the surface] + if(l_sfcobserror_ramp_t) then + ramp=min(max(((rlow-1.0_r_kind)/(15.0_r_kind-1.0_r_kind)),0.0_r_kind),1.0_r_kind) + else + ramp=rlow + endif else - ramp=rlow + rlow = zero + ramp = zero endif rhgh=max(zero,dpres-rsigp-r0_001) @@ -795,11 +881,28 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav if(rlow/=zero) awork(2) = awork(2) + one if(rhgh/=zero) awork(3) = awork(3) + one end if - - ratio_errors=error/(data(ier,i)+drpx+1.0e6_r_kind*rhgh+r8*ramp) + + ! inflate error for uncertainty in the terrain adjustment + lapse_error = 0. + if ( hofx_2m_sfcfile .and. sfctype) then + if (abs(delta_z)300 m. do not assim. + ! inflate obs error to account for error in lapse_rate + ! also include some representativity error here (assuming + ! delta_z ~ heterogeneity) + lapse_error = abs(lapse_error_frac*T_lapse*delta_z) + else + muse(i)=.false. + endif + endif + +! note: for hofx_2m_sfcfile, three middle terms in denominator will be zero +! (elevation mistmatch between obs and model dealt with elsewhere) + ratio_errors=error/(data(ier,i)+drpx+1.0e6_r_kind*rhgh+r8*ramp + lapse_error) + + ! Compute innovation - if(i_use_2mt4b>0 .and. sfctype) then + if(sfctype .and. (( i_use_2mt4b>0) .or. hofx_2m_sfcfile ) ) then ddiff = tob-tges2m else ddiff = tob-tges @@ -1411,7 +1514,7 @@ subroutine init_vars_ write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus call stop2(999) endif - if(i_use_2mt4b>0) then + if(i_use_2mt4b>0 .or. hofx_2m_sfcfile) then ! get t2m ... varname='t2m' call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus) @@ -1676,12 +1779,12 @@ subroutine contents_netcdf_diag_(odiag) call nc_diag_metadata("Observation_Subtype", icsubtype(ikx) ) call nc_diag_metadata("Latitude", sngl(data(ilate,i)) ) call nc_diag_metadata("Longitude", sngl(data(ilone,i)) ) - call nc_diag_metadata("Station_Elevation", sngl(data(istnelv,i)) ) + call nc_diag_metadata("Station_Elevation", sngl(data(istnelv,i)) ) ! this is the obs elevation, after being interpolated to the model (=model height) call nc_diag_metadata("Pressure", sngl(prest) ) - call nc_diag_metadata("Height", sngl(data(iobshgt,i)) ) + call nc_diag_metadata("Height", sngl(data(iobshgt,i)) ) ! this is the obs height (= stn elevation, before being corrected) call nc_diag_metadata("Time", sngl(dtime-time_offset)) call nc_diag_metadata("Prep_QC_Mark", sngl(data(iqc,i)) ) - call nc_diag_metadata("Setup_QC_Mark", sngl(data(iqt,i)) ) + call nc_diag_metadata("Setup_QC_Mark", sngl(data(iqt,i)) ) ! this is the virtual temp flag call nc_diag_metadata("Prep_Use_Flag", sngl(data(iuse,i)) ) if(muse(i)) then call nc_diag_metadata("Analysis_Use_Flag", sngl(one) ) @@ -1693,7 +1796,12 @@ subroutine contents_netcdf_diag_(odiag) call nc_diag_metadata("Errinv_Input", sngl(errinv_input) ) call nc_diag_metadata("Errinv_Adjust", sngl(errinv_adjst) ) call nc_diag_metadata("Errinv_Final", sngl(errinv_final) ) - call nc_diag_metadata("Observation", sngl(data(itob,i)) ) + if (hofx_2m_sfcfile ) then + call nc_diag_metadata("Observation", sngl(tob) ) + call nc_diag_metadata("Observation_Before_Elev_Correction", sngl(data(itob,i)) ) + else + call nc_diag_metadata("Observation", sngl(data(itob,i)) ) + endif call nc_diag_metadata("Obs_Minus_Forecast_adjusted", sngl(ddiff) ) call nc_diag_metadata("Obs_Minus_Forecast_unadjusted", sngl(tob-tges) ) call nc_diag_metadata("Forecast_unadjusted", sngl(tges)) @@ -1735,9 +1843,12 @@ subroutine contents_netcdf_diag_(odiag) call nc_diag_data2d("ObsDiagSave_obssen", odiag%obssen ) endif - if (twodvar_regional .or. l_obsprvdiag) then + if (twodvar_regional .or. l_obsprvdiag .or. hofx_2m_sfcfile ) then call nc_diag_metadata("Dominant_Sfc_Type", data(idomsfc,i) ) +! this is the model height interpolated to the obs location in read_prepbufr call nc_diag_metadata("Model_Terrain", data(izz,i) ) + endif + if (twodvar_regional .or. l_obsprvdiag) then r_prvstg = data(iprvd,i) call nc_diag_metadata("Provider_Name", c_prvstg ) r_sprvstg = data(isprvd,i) diff --git a/src/gsi/update_guess.f90 b/src/gsi/update_guess.f90 index a90e7a19d6..ce811ae90e 100644 --- a/src/gsi/update_guess.f90 +++ b/src/gsi/update_guess.f90 @@ -113,7 +113,7 @@ subroutine update_guess(sval,sbias) use mpimod, only: mype use constants, only: zero,one,fv,max_varname_length,qmin,qcmin,tgmin,& r100,one_tenth,tiny_r_kind - use jfunc, only: iout_iter,bcoption,tsensible,clip_supersaturation,superfact + use jfunc, only: iout_iter,bcoption,tsensible,clip_supersaturation,superfact,hofx_2m_sfcfile use gridmod, only: lat2,lon2,nsig,& regional,twodvar_regional,regional_ozone,& l_reg_update_hydro_delz @@ -454,7 +454,7 @@ subroutine update_guess(sval,sbias) endif call gsd_update_soil_tq(tinc_1st,is_t,qinc_1st,is_q,it) endif ! l_gsd_soilTQ_nudge - if (i_use_2mt4b > 0 .and. is_t>0) then + if ( (i_use_2mt4b > 0.or. hofx_2m_sfcfile) .and. is_t>0) then do j=1,lon2 do i=1,lat2 tinc_1st(i,j)=p_tv(i,j,1) @@ -462,7 +462,7 @@ subroutine update_guess(sval,sbias) end do call gsd_update_t2m(tinc_1st,it) endif ! l_gsd_t2m_adjust - if (i_use_2mq4b > 0 .and. is_q>0) then + if ( (i_use_2mq4b > 0.or. hofx_2m_sfcfile) .and. is_q>0) then do j=1,lon2 do i=1,lat2 qinc_1st(i,j)=p_q(i,j,1) From 9473171843bf389893db5e3f3087a43957caaa5e Mon Sep 17 00:00:00 2001 From: "Guoqing.Ge" Date: Tue, 25 Oct 2022 20:42:20 -0600 Subject: [PATCH 2/2] change tqm and obserr for surface T obs --- src/gsi/read_prepbufr.f90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/gsi/read_prepbufr.f90 b/src/gsi/read_prepbufr.f90 index 21d73b41ce..6d0ab1b365 100644 --- a/src/gsi/read_prepbufr.f90 +++ b/src/gsi/read_prepbufr.f90 @@ -1919,6 +1919,10 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& ! CSD - temporary hack ( move to prepbufr pre-processing later) ! Over-ride QM=9 for sfc obs if (sfctype .and. hofx_2m_sfcfile ) then + if (tob) then!hardwire 187 obs error to 1.2 for now + tqm(k)=2 + obserr(3,k)=1.2 + endif if (tob .and. qm == 9 ) qm = 2 ! 2=not checked if (qob .and. qm == 9 ) qm = 2 ! 2=not checked endif