Skip to content

Commit

Permalink
Feature/enkf q2 (#568)
Browse files Browse the repository at this point in the history 
Resubmit PR for cleaning the unneeded item in the previous PR 
In correspondence to the EMC GSI Issue#566, this PR contains a quick
adding of the clipping of negative values of sphum (q) in the analysis
of FV3-LAM EnKF.
This part of codes are not tested in the current GSI regression tests,
which, hence, are not run.
The current codes are verified using local FV3-LAM case. It is found the
differences from this changes exist for sphum ( maximum values about
0.003 (units) and 0.3 K for T (sensible T). The latter is because the
sphum would be used when the analysis TV is converted to T. All
differences are on spontaneous points and values are reasonable as
expected.
Hence, the code is regarded as verified.

Fixes #566

**DUE DATE for this PR is 6/15/2023.** If this PR is not merged into
`develop` by this date, the PR will be closed and returned to the
developer.
  • Loading branch information
@TingLei-daprediction
TingLei-daprediction authored Jun 27, 2023
1 parent 8735959 commit 00cac54
Show file tree
Hide file tree
Showing 2 changed files with 22 additions and 52 deletions.
72 changes: 21 additions & 51 deletions src/enkf/gridio_fv3reg.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ module gridio
use params, only: nlevs, cliptracers, datapath, arw, nmm, datestring
use params, only: nx_res,ny_res,nlevs,ntiles,l_fv3reg_filecombined,&
fv3_io_layout_nx,fv3_io_layout_ny,nanals
use params, only: pseudo_rh, l_use_enkf_directZDA
use params, only: pseudo_rh
use mpeu_util, only: getindex
use read_fv3regional_restarts,only:read_fv3_restart_data1d,read_fv3_restart_data2d
use read_fv3regional_restarts,only:read_fv3_restart_data3d,read_fv3_restart_data4d
Expand Down Expand Up @@ -694,7 +694,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid

ps_ind = getindex(vars2d, 'ps') ! Ps (2D)


clip=tiny(clip)
!----------------------------------------------------------------------
if (nbackgrounds > 1) then
write(6,*)'gridio/writegriddata: writing multiple backgrounds not yet supported'
Expand Down Expand Up @@ -853,6 +853,8 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid
varstrname = 'sphum'
call fv3lamfile%get_idfn(varstrname,file_id,fv3filename)
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,qworkvar3d)
!enforce lower positive bound (clip) to replace negative hydrometers
if ( cliptracers ) where (qworkvar3d < clip) qworkvar3d = clip
tvworkvar3d=tsenworkvar3d*(one+fv*qworkvar3d)
tvworkvar3d=tvworkvar3d+workinc3d
if(q_ind > 0) then
Expand All @@ -866,6 +868,8 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid
enddo
enddo
qworkvar3d=qworkvar3d+workinc3d
!enforce lower positive bound (clip) to replace negative q
if ( cliptracers ) where (qworkvar3d < clip) qworkvar3d = clip
endif
tsenworkvar3d=tvworkvar3d/(one+fv*qworkvar3d)
varstrname = 'T'
Expand Down Expand Up @@ -932,10 +936,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid
enddo
enddo
workvar3d=workvar3d+workinc3d
if ( l_use_enkf_directZDA .and. cliptracers ) then ! set cliptracers to remove negative hydrometers
clip = tiny(workvar3d(1,1,1))
where (workvar3d < clip) workvar3d = clip
end if
if ( cliptracers ) where (workvar3d < clip) workvar3d = clip
call write_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d)

endif
Expand All @@ -955,10 +956,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid
enddo
enddo
workvar3d=workvar3d+workinc3d
if ( l_use_enkf_directZDA .and. cliptracers ) then ! set cliptracers to remove negative hydrometers
clip = tiny(workvar3d(1,1,1))
where (workvar3d < clip) workvar3d = clip
end if
if ( cliptracers ) where (workvar3d < clip) workvar3d = clip
call write_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d)

endif
Expand All @@ -978,10 +976,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid
enddo
enddo
workvar3d=workvar3d+workinc3d
if ( l_use_enkf_directZDA .and. cliptracers ) then ! set cliptracers to remove negative hydrometers
clip = tiny(workvar3d(1,1,1))
where (workvar3d < clip) workvar3d = clip
end if
if ( cliptracers ) where (workvar3d < clip) workvar3d = clip
call write_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d)

endif
Expand All @@ -1001,10 +996,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid
enddo
enddo
workvar3d=workvar3d+workinc3d
if ( l_use_enkf_directZDA .and. cliptracers ) then ! set cliptracers to remove negative hydrometers
clip = tiny(workvar3d(1,1,1))
where (workvar3d < clip) workvar3d = clip
end if
if ( cliptracers ) where (workvar3d < clip) workvar3d = clip
call write_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d)

endif
Expand All @@ -1024,10 +1016,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid
enddo
enddo
workvar3d=workvar3d+workinc3d
if ( l_use_enkf_directZDA .and. cliptracers ) then ! set cliptracers to remove negative hydrometers
clip = tiny(workvar3d(1,1,1))
where (workvar3d < clip) workvar3d = clip
end if
if ( cliptracers ) where (workvar3d < clip) workvar3d = clip
call write_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d)

endif
Expand All @@ -1046,10 +1035,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid
enddo
enddo
workvar3d=workvar3d+workinc3d
if ( l_use_enkf_directZDA .and. cliptracers ) then ! set cliptracers to remove negative hydrometers
clip = tiny(workvar3d(1,1,1))
where (workvar3d < clip) workvar3d = clip
end if
if ( cliptracers ) where (workvar3d < clip) workvar3d = clip
call write_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d)

endif
Expand Down Expand Up @@ -1888,7 +1874,7 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,vargrid,no_inflat

ps_ind = getindex(vars2d, 'ps') ! Ps (2D)


clip=tiny(clip)
allocate(my_neb(4))
!----------------------------------------------------------------------
if (nbackgrounds > 1) then
Expand Down Expand Up @@ -2173,6 +2159,7 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,vargrid,no_inflat
enddo

if(iope ==0 ) then
if ( cliptracers ) where (qworkvar3d < clip) qworkvar3d = clip
tvworkvar3d=tsenworkvar3d*(one+fv*qworkvar3d)
tvworkvar3d=tvworkvar3d+workinc3d
if(q_ind > 0) then
Expand All @@ -2186,6 +2173,7 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,vargrid,no_inflat
enddo
enddo
qworkvar3d=qworkvar3d+workinc3d
if ( cliptracers ) where (qworkvar3d < clip) qworkvar3d = clip
endif
tsenworkvar3d=tvworkvar3d/(one+fv*qworkvar3d)
endif
Expand Down Expand Up @@ -2281,10 +2269,7 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,vargrid,no_inflat
enddo
enddo
workvar3d=workvar3d+workinc3d
if ( l_use_enkf_directZDA .and. cliptracers ) then ! set cliptracers to remove negative hydrometers
clip = tiny(workvar3d(1,1,1))
where (workvar3d < clip) workvar3d = clip
end if
if ( cliptracers ) where (workvar3d < clip) workvar3d = clip
endif
do k=1,nlevs
call mpi_scatterv(workvar3d(:,:,k),recvcounts2d,displs2d,mpi_real4,&
Expand Down Expand Up @@ -2315,10 +2300,7 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,vargrid,no_inflat
enddo
enddo
workvar3d=workvar3d+workinc3d
if ( l_use_enkf_directZDA .and. cliptracers ) then ! set cliptracers to remove negative hydrometers
clip = tiny(workvar3d(1,1,1))
where (workvar3d < clip) workvar3d = clip
end if
if ( cliptracers ) where (workvar3d < clip) workvar3d = clip
endif
do k=1,nlevs
call mpi_scatterv(workvar3d(:,:,k),recvcounts2d,displs2d,mpi_real4,&
Expand Down Expand Up @@ -2349,10 +2331,7 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,vargrid,no_inflat
enddo
enddo
workvar3d=workvar3d+workinc3d
if ( l_use_enkf_directZDA .and. cliptracers ) then ! set cliptracers to remove negative hydrometers
clip = tiny(workvar3d(1,1,1))
where (workvar3d < clip) workvar3d = clip
end if
if ( cliptracers ) where (workvar3d < clip) workvar3d = clip
endif
do k=1,nlevs
call mpi_scatterv(workvar3d(:,:,k),recvcounts2d,displs2d,mpi_real4,&
Expand Down Expand Up @@ -2383,10 +2362,7 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,vargrid,no_inflat
enddo
enddo
workvar3d=workvar3d+workinc3d
if ( l_use_enkf_directZDA .and. cliptracers ) then ! set cliptracers to remove negative hydrometers
clip = tiny(workvar3d(1,1,1))
where (workvar3d < clip) workvar3d = clip
end if
if ( cliptracers ) where (workvar3d < clip) workvar3d = clip
endif
do k=1,nlevs
call mpi_scatterv(workvar3d(:,:,k),recvcounts2d,displs2d,mpi_real4,&
Expand Down Expand Up @@ -2417,10 +2393,7 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,vargrid,no_inflat
enddo
enddo
workvar3d=workvar3d+workinc3d
if ( l_use_enkf_directZDA .and. cliptracers ) then ! set cliptracers to remove negative hydrometers
clip = tiny(workvar3d(1,1,1))
where (workvar3d < clip) workvar3d = clip
end if
if ( cliptracers ) where (workvar3d < clip) workvar3d = clip
endif
do k=1,nlevs
call mpi_scatterv(workvar3d(:,:,k),recvcounts2d,displs2d,mpi_real4,&
Expand Down Expand Up @@ -2450,10 +2423,7 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,vargrid,no_inflat
enddo
enddo
workvar3d=workvar3d+workinc3d
if ( l_use_enkf_directZDA .and. cliptracers ) then ! set cliptracers to remove negative hydrometers
clip = tiny(workvar3d(1,1,1))
where (workvar3d < clip) workvar3d = clip
end if
if ( cliptracers ) where (workvar3d < clip) workvar3d = clip
endif
do k=1,nlevs
call mpi_scatterv(workvar3d(:,:,k),recvcounts2d,displs2d,mpi_real4,&
Expand Down
2 changes: 1 addition & 1 deletion src/enkf/innovstats.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -213,7 +213,7 @@ subroutine print_innovstats(obfit,obsprd)
call printstats(' all gps',sumgps_nh,biasq_nh,sumgps_spread_nh,sumgps_oberr_nh,nobsgps_nh,&
sumgps_sh,biasgps_sh,sumgps_spread_sh,sumgps_oberr_sh,nobsgps_sh,&
sumgps_tr,biasgps_tr,sumgps_spread_tr,sumgps_oberr_tr,nobsgps_tr)
call printstats(' all dbz',sumdbz_nh,biasq_nh,sumdbz_spread_nh,sumdbz_oberr_nh,nobsdbz_nh,&
call printstats(' all dbz',sumdbz_nh,biasdbz_nh,sumdbz_spread_nh,sumdbz_oberr_nh,nobsdbz_nh,&
sumdbz_sh,biasdbz_sh,sumdbz_spread_sh,sumdbz_oberr_sh,nobsdbz_sh,&
sumdbz_tr,biasdbz_tr,sumdbz_spread_tr,sumdbz_oberr_tr,nobsdbz_tr)
call printstats(' all rw',sumrw_nh,biasq_nh,sumrw_spread_nh,sumrw_oberr_nh,nobsrw_nh,&
Expand Down

0 comments on commit 00cac54

Please sign in to comment.