Skip to content

Commit

Permalink
Fix a bug with q when tsen is used as control variable in fv3 regiona…
Browse files Browse the repository at this point in the history 
…l EnKF analysis (NOAA-EMC#6)

* For fv3-lam enkf, fix a bug with q when tsen is used as control variable

* For fv3-lam enkf, fix a bug with q when tsen is used as control variable

Co-authored-by: TingLei-NOAA <Ting.Lei@noaa.gov>
  • Loading branch information
@BinLiu-NOAA @TingLei-NOAA
BinLiu-NOAA and TingLei-NOAA authored Mar 23, 2021
1 parent 92dd834 commit 6d370e8
Showing 1 changed file with 103 additions and 105 deletions.
208 changes: 103 additions & 105 deletions src/enkf/gridio_fv3reg.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,19 +5,15 @@ module gridio
!$$$ Module documentation block
!
! This module contains various routines to ingest and update
! variables from Weather Research and Forecasting (WRF) model Advanced
! Research WRF (ARW) and Non-hydrostatic Mesoscale Model (NMM) dynamical
! cores which are required by the Ensemble Kalman Filter (ENKF) currently
! designed for operations within the National Centers for Environmental
! Prediction (NCEP) Global Forecasting System (GFS)
! variables from FV3-LAM warmstart files
!
! prgmmr: Winterbottom org: ESRL/PSD1 date: 2011-11-30
!
! program history log:
!
! 2011-11-30 Winterbottom - Initial version.
!
! 2019-01- Ting -- modified for fv3sar
! 2019-01- Ting -- modified for fv3-lam
! attributes:
! language: f95
!
Expand Down Expand Up @@ -191,17 +187,9 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,f
deallocate(vworkvar3d)

endif

if (tv_ind > 0.or.tsen_ind) then
allocate(tsenworkvar3d(nx_res,ny_res,nlevs))
varstrname = 'T'
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,tsenworkvar3d)
varstrname = 'sphum'
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,qworkvar3d)


if (q_ind > 0) then
varstrname = 'sphum'
if (q_ind > 0) then
varstrname = 'sphum'
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,qworkvar3d)
do k=1,nlevs
nn = nn_tile0
do j=1,ny_res
Expand All @@ -217,18 +205,27 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,f
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo

endif
if(tv_ind > 0) then
endif
if (tv_ind > 0.or.tsen_ind) then
allocate(tsenworkvar3d(nx_res,ny_res,nlevs))
varstrname = 'T'
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,tsenworkvar3d)
if(tv_ind > 0) then
if(.not. (q_ind > 0)) then
varstrname = 'sphum'
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,qworkvar3d)
endif
do k=1,nlevs
do j=1,ny_res
do i=1,nx_res
workvar3d(i,j,k)=tsenworkvar3d(i,j,k)*(one+fv*qworkvar3d(i,j,k))
enddo
enddo
enddo
tvworkvar3d=workvar3d
tvworkvar3d=workvar3d
else! tsen_id >0
workvar3d=tsenworkvar3d
tvworkvar3d=workvar3d*(one+fv*qworkvar3d(i,j,k))
endif
tmp_ind=max(tv_ind,tsen_ind) !then can't be both >0
do k=1,nlevs
Expand All @@ -242,7 +239,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,f
enddo
do k = levels(tmp_ind-1)+1, levels(tmp_ind)
if (nproc .eq. 0) then
write(6,*) 'READFVregional : t ', &
write(6,*) 'READFVregional : tv or tsen ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
endif
enddo
Expand Down Expand Up @@ -361,12 +358,12 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,f
end subroutine readgriddata

!========================================================================
! readgriddata_nmm.f90: read FV3-Lam state or control vector
! readgriddata: read fv3-lam state or control vector
!-------------------------------------------------------------------------


!========================================================================
! writegriddata.f90: write FV3-LAM analysis
! writegriddata: write fv3-lam
!-------------------------------------------------------------------------

subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid,no_inflate_flag)
Expand Down Expand Up @@ -404,7 +401,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid
real(r_single), dimension(:,:), allocatable ::pswork
real(r_single), dimension(:,:,:), allocatable ::workvar3d,workinc3d,workinc3d2,uworkvar3d,&
vworkvar3d,tvworkvar3d,tsenworkvar3d,&
workprsi,qworkvar3d
workprsi,qworkvar3d,qbgworkvar3d

!----------------------------------------------------------------------
! Define variables required by for extracting netcdf variable
Expand Down Expand Up @@ -438,7 +435,6 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid

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


!----------------------------------------------------------------------
if (nbackgrounds > 1) then
write(6,*)'gridio/writegriddata: writing multiple backgrounds not yet supported'
Expand All @@ -452,6 +448,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid
allocate(workinc3d(nx_res,ny_res,nlevs),workinc3d2(nx_res,ny_res,nlevsp1))
allocate(workvar3d(nx_res,ny_res,nlevs))
allocate(qworkvar3d(nx_res,ny_res,nlevs))
allocate(qbgworkvar3d(nx_res,ny_res,nlevs))
allocate(tvworkvar3d(nx_res,ny_res,nlevs))


Expand Down Expand Up @@ -503,96 +500,101 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid
allocate(vworkvar3d(nx_res+1,ny_res,nlevs))

call read_fv3_restart_data3d(varstrname,fv3filename,file_id,vworkvar3d)
do k=1,nlevs
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
workinc3d(i,j,k)=vargrid(nn,levels(v_ind-1)+k,nb,ne)
enddo
enddo
enddo
vworkvar3d(1:nx_res,:,:)=vworkvar3d(1:nx_res,:,:)+workinc3d
vworkvar3d(nx_res+1,:,:)=vworkvar3d(nx_res,:,:)
do k=1,nlevs
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
workinc3d(i,j,k)=vargrid(nn,levels(v_ind-1)+k,nb,ne)
enddo
enddo
enddo
vworkvar3d(1:nx_res,:,:)=vworkvar3d(1:nx_res,:,:)+workinc3d
vworkvar3d(nx_res+1,:,:)=vworkvar3d(nx_res,:,:)
call write_fv3_restart_data3d(varstrname,fv3filename,file_id,vworkvar3d)

deallocate(vworkvar3d)
endif
if (tv_ind > 0.or.tsen_ind>0 ) then

varstrname = 'T'
if(tsen_ind>0) then
do k=1,nlevs
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
workinc3d(i,j,k)=vargrid(nn,levels(tsen_ind-1)+k,nb,ne)
enddo
enddo
enddo
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d)
workvar3d=workvar3d+workinc3d
call write_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d)
else ! tv_ind >0
do k=1,nlevs
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
workinc3d(i,j,k)=vargrid(nn,levels(tv_ind-1)+k,nb,ne)
enddo
enddo
enddo

varstrname = 'T'
allocate(tsenworkvar3d(nx_res,ny_res,nlevs))
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,tsenworkvar3d)
varstrname = 'sphum'
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,qworkvar3d)
tvworkvar3d=tsenworkvar3d*(one+fv*qworkvar3d)
tvworkvar3d=tvworkvar3d+workinc3d
if(q_ind > 0) then
do k=1,nlevs
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
workinc3d(i,j,k)=vargrid(nn,levels(q_ind-1)+k,nb,ne)
enddo
enddo
enddo
qworkvar3d=qworkvar3d+workinc3d
endif
tsenworkvar3d=tvworkvar3d/(one+fv*qworkvar3d)
varstrname = 'T'
call write_fv3_restart_data3d(varstrname,fv3filename,file_id,tsenworkvar3d)
do k=1,nlevs
if (nproc .eq. 0) &
write(6,*) 'WRITEregional : T ', &
& k, minval(tsenworkvar3d(:,:,k)), maxval(tsenworkvar3d(:,:,k))
enddo




if(q_ind>0) then
if(q_ind>0) then

varstrname='sphum'
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,qbgworkvar3d)
do k=1,nlevs
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
workinc3d(i,j,k)=vargrid(nn,levels(q_ind-1)+k,nb,ne)
enddo
enddo
enddo
qworkvar3d=qbgworkvar3d+workinc3d

call write_fv3_restart_data3d(varstrname,fv3filename,file_id,qworkvar3d)
do k=1,nlevs
if (nproc .eq. 0) &
write(6,*) 'WRITEregional : sphum ', &
& k, minval(qworkvar3d(:,:,k)), maxval(qworkvar3d(:,:,k))
enddo
endif


end if



if (tv_ind > 0.or.tsen_ind>0 ) then

varstrname = 'T'
if(tsen_ind>0) then
do k=1,nlevs
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
workinc3d(i,j,k)=vargrid(nn,levels(tsen_ind-1)+k,nb,ne)
enddo
enddo
enddo
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d)
workvar3d=workvar3d+workinc3d
call write_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d)
else ! tv_ind >0
do k=1,nlevs
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
workinc3d(i,j,k)=vargrid(nn,levels(tv_ind-1)+k,nb,ne)
enddo
enddo
enddo

varstrname = 'T'
allocate(tsenworkvar3d(nx_res,ny_res,nlevs))
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,tsenworkvar3d)
if(.not. (q_ind > 0) ) then
varstrname = 'sphum'
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,qbgworkvar3d)
endif
tvworkvar3d=tsenworkvar3d*(one+fv*qbgworkvar3d)
tvworkvar3d=tvworkvar3d+workinc3d
if(.not. ( q_ind > 0)) then
tsenworkvar3d=tvworkvar3d/(one+fv*qbgworkvar3d)
else
tsenworkvar3d=tvworkvar3d/(one+fv*qworkvar3d)
endif
varstrname = 'T'
call write_fv3_restart_data3d(varstrname,fv3filename,file_id,tsenworkvar3d)
do k=1,nlevs
if (nproc .eq. 0) &
write(6,*) 'WRITEregional : T ', &
& k, minval(tsenworkvar3d(:,:,k)), maxval(tsenworkvar3d(:,:,k))
enddo

deallocate(tsenworkvar3d)
endif
deallocate(tsenworkvar3d)
endif !if tsens else tv

endif

if (oz_ind > 0) then
varstrname = 'o3mr'

Expand Down Expand Up @@ -621,9 +623,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid
workprsi(:,:,i)=workvar3d(:,:,i)*0.01_r_kind+workprsi(:,:,i+1)
enddo



nn = nn_tile0
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
Expand Down Expand Up @@ -657,16 +657,14 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid
myname_,'close '//trim(filename) )


!----------------------------------------------------------------------
! update time stamp is to be considered NSTART_HOUR in NMM (HWRF) restart file.
!======================================================================
end do ! tiles
if(allocated(workinc3d)) deallocate(workinc3d)
if(allocated(workinc3d2)) deallocate(workinc3d2)
if(allocated(workprsi)) deallocate(workprsi)
if(allocated(pswork)) deallocate(pswork)
if(allocated(tvworkvar3d)) deallocate(tvworkvar3d)
if(allocated(qworkvar3d)) deallocate(qworkvar3d)
if(allocated(qbgworkvar3d)) deallocate(qbgworkvar3d)



Expand Down

0 comments on commit 6d370e8

Please sign in to comment.