Skip to content

Commit

Permalink
Merge pull request NOAA-EMC#9 from guoqing-noaa/bugfix
Browse files Browse the repository at this point in the history 
bug fix on the latest merge of EMC develop and 2mDA
  • Loading branch information
@jswhit
jswhit authored Oct 26, 2022
2 parents 5f08c29 + 2046b38 commit 0f00543
Show file tree
Hide file tree
Showing 4 changed files with 205 additions and 88 deletions.
269 changes: 191 additions & 78 deletions src/gsi/general_read_gfsatm.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1959,7 +1959,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

Expand All @@ -1977,6 +1977,7 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, &
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

Expand Down Expand Up @@ -2011,8 +2012,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
Expand All @@ -2026,6 +2026,20 @@ 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
Expand Down Expand Up @@ -2064,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
Expand Down Expand Up @@ -2097,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) )
Expand Down Expand Up @@ -2140,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))
Expand All @@ -2202,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
Expand All @@ -2228,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
Expand Down Expand Up @@ -2575,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)
Expand All @@ -2597,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
Expand Down
7 changes: 3 additions & 4 deletions src/gsi/gsd_update_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -516,9 +516,8 @@ subroutine gsd_update_t2m(tinc,it)
!$$$
use kinds, only: r_kind,i_kind
use jfunc, only: tsensible
use constants, only: zero,one,fv,rd_over_cp_mass,one_tenth
use gridmod, only: lat2,lon2,aeta1_ll,pt_ll,aeta2_ll
use guess_grids, only:ges_prsi
use constants, only: zero,one,fv,one_tenth
use gridmod, only: lat2,lon2
use constants, only: half

implicit none
Expand Down Expand Up @@ -555,7 +554,7 @@ subroutine gsd_update_t2m(tinc,it)
if(ihaveq/=0) cycle
dth2=tinc(i,j)/(one+fv*ges_q(i,j,1))
endif
! do not need to convert sensible temperature to potential temperature

ges_t2m(i,j) = ges_t2m(i,j) + dth2
end do
end do
Expand Down
6 changes: 5 additions & 1 deletion src/gsi/read_prepbufr.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1915,9 +1915,13 @@ 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)
! 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
Expand Down
Loading

0 comments on commit 0f00543

Please sign in to comment.