Skip to content

Commit

Permalink
Fix write_incr (NOAA-EMC#559)
Browse files Browse the repository at this point in the history
  • Loading branch information
@jderber-NOAA
jderber-NOAA authored and ting lei committed May 24, 2023
1 parent 60d3aa8 commit ddfe627
Show file tree
Hide file tree
Showing 4 changed files with 104 additions and 20 deletions.
4 changes: 2 additions & 2 deletions src/enkf/controlvec.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,6 @@ module controlvec
use mpeu_util, only: gettablesize, gettable, getindex
use constants, only: max_varname_length
implicit none

private

public :: read_control, write_control, controlvec_cleanup, init_controlvec
Expand Down Expand Up @@ -224,7 +223,8 @@ subroutine read_control()
call readgriddata(nanal1(nproc),nanal2(nproc),cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,nbackgrounds, &
fgfileprefixes,fgsfcfileprefixes,reducedgrid,grdin,qsat)
end if
!print *,'min/max qsat',nanal,'=',minval(qsat),maxval(qsat)
call mpi_barrier(mpi_comm_world,ierr)
write(6,*)'thinkdeb min/max qsat','=',minval(qsat),maxval(qsat)
if (use_qsatensmean) then
allocate(qsatmean(npts,nlevs,nbackgrounds))
allocate(qsat_tmp(npts))
Expand Down
89 changes: 75 additions & 14 deletions src/enkf/gridio_fv3reg.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,7 @@ module gridio
use netcdf_mod,only: nc_check
use mpimod, only: mpi_comm_world, mpi_sum, mpi_real4, mpi_real8, mpi_rtype
use mpimod, only: mpi_status_size
use mpisetup,only: mpi_wtime

implicit none

Expand Down Expand Up @@ -133,10 +134,13 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,
integer (i_kind):: ps_ind, sst_ind
integer (i_kind):: tmp_ind,ifile
logical (i_kind):: ice
real tot1,tot2
real t1,t2

!======================================================================
write (6,*)"The input fileprefix, reducedgrid are not used in the current implementation", &
fileprefixes, reducedgrid
tot1=mpi_wtime()
nlevsp1=nlevs+1
u_ind = getindex(vars3d, 'u') !< indices in the state var arrays
v_ind = getindex(vars3d, 'v') ! U and V (3D)
Expand Down Expand Up @@ -217,7 +221,9 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,

!----------------------------------------------------------------------
! Update u and v variables (same for NMM and ARW)

write(6,*)'thinkdeb readgridata 10'
call flush(6)
t1=mpi_wtime()
if (u_ind > 0) then
allocate(uworkvar3d(nx_res,ny_res+1,nlevs))
varstrname = 'u'
Expand All @@ -241,6 +247,9 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,

deallocate(uworkvar3d)
endif
t2=mpi_wtime()
write(6,*)'thinkdeb read u t ',t2-t1
t1=mpi_wtime()
if (v_ind > 0) then
allocate(vworkvar3d(nx_res+1,ny_res,nlevs))
varstrname = 'v'
Expand All @@ -264,6 +273,9 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,
deallocate(vworkvar3d)

endif
t2=mpi_wtime()
write(6,*)'thinkdeb read v t ',t2-t1
t1=mpi_wtime()
if (w_ind > 0) then
varstrname = 'W'
call fv3lamfile%get_idfn(varstrname,file_id,fv3filename)
Expand All @@ -286,8 +298,13 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,
enddo

endif
t2=mpi_wtime()
write(6,*)'thinkdeb read w t ',t2-t1
write(6,*)'thinkdeb readgridata 11'
call flush(6)

if (tv_ind > 0 .or. tsen_ind > 0) then
t1=mpi_wtime()
allocate(tsenworkvar3d(nx_res,ny_res,nlevs))
varstrname = 'T'
call fv3lamfile%get_idfn(varstrname,file_id,fv3filename)
Expand All @@ -296,6 +313,9 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,
call fv3lamfile%get_idfn(varstrname,file_id,fv3filename)
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,qworkvar3d)

t2=mpi_wtime()
write(6,*)'thinkdeb read q 1 t ',t2-t1
t1=mpi_wtime()

if (q_ind > 0) then
varstrname = 'sphum'
Expand All @@ -317,6 +337,9 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,
enddo

endif
t2=mpi_wtime()
write(6,*)'thinkdeb process q t ',t2-t1
t1=mpi_wtime()
if(tv_ind > 0) then
!$omp parallel do schedule(dynamic,1) private(k,j,i)
do k=1,nlevs
Expand All @@ -330,16 +353,18 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,
do k=1,nlevs
do j=1,ny_res
do i=1,nx_res
tvworkvar3d=workvar3d
tvworkvar3d(i,j,k)=workvar3d(i,j,k)
enddo
enddo
enddo
t2=mpi_wtime()
write(6,*)'thinkdeb tv_ind t ',t2-t1
else! tsen_id >0
!$omp parallel do schedule(dynamic,1) private(k,j,i)
do k=1,nlevs
do j=1,ny_res
do i=1,nx_res
workvar3d=tsenworkvar3d
workvar3d(i,j,k)=tsenworkvar3d(i,j,k)
enddo
enddo
enddo
Expand All @@ -365,6 +390,8 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,

endif
if(allocated(tsenworkvar3d)) deallocate(tsenworkvar3d)
write(6,*)'thinkdeb readgridata 17'
call flush(6)



Expand Down Expand Up @@ -498,6 +525,8 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,

endif

write(6,*)'thinkdeb readgridata 19'
call flush(6)
if (qnr_ind > 0) then
varstrname = 'rain_nc'
call fv3lamfile%get_idfn(varstrname,file_id,fv3filename)
Expand All @@ -523,11 +552,12 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,
varstrname = 'ref_f3d'
call fv3lamfile%get_idfn(varstrname,file_id,fv3filename)
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d)
!$omp parallel do schedule(dynamic,1) private(k,j,i,nn)
do k=1,nlevs
nn = nn_tile0
!clt nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
nn=nn_tile0+(j-1)*nx_res+i
vargrid(nn,levels(dbz_ind-1)+k,nb,ne)=max(workvar3d(i,j,nlevs+1-k),0.0_r_kind)
enddo
enddo
Expand Down Expand Up @@ -582,7 +612,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,




!$omp parallel do schedule(dynamic,1) private(k,j,i)
do k=1,nlevs
do j=1,ny_res
do i=1,nx_res
Expand All @@ -609,6 +639,8 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,



write(6,*)'thinkdeb readgridata 19'
call flush(6)


if(allocated(workprsi)) deallocate(workprsi)
Expand Down Expand Up @@ -636,7 +668,9 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,

end do backgroundloop ! loop over backgrounds to read in
end do ensmemloop ! loop over ens members to read in

tot2=mpi_wtime()
write(6,*)'thinkdeb readgridata 20, tot t ',tot2-tot1
call flush(6)
return

end subroutine readgriddata
Expand Down Expand Up @@ -689,6 +723,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid
workprsi,qworkvar3d

real(r_single) :: clip
real :: t1,t2

!----------------------------------------------------------------------
! Define variables required by for extracting netcdf variable
Expand All @@ -705,7 +740,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid



write(6,*)"anlfileprefixes, fgfileprefixes are not used in the current implementation", &
write(6,*)"anlfileprefixesthinkdeb, fgfileprefixes are not used in the current implementation", &
anlfileprefixes, fgfileprefixes
write(6,*)"the no_inflate_flag is not used in the currrent implementation ",no_inflate_flag
!----------------------------------------------------------------------
Expand Down Expand Up @@ -1003,6 +1038,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid
enddo
enddo
enddo
!$omp parallel do schedule(dynamic,1) private(k,j,i)
do k=1,nlevs
do j=1,ny_res
do i=1,nx_res
Expand Down Expand Up @@ -1063,6 +1099,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid
enddo
enddo
enddo
t1=mpi_wtime()
!$omp parallel do schedule(dynamic,1) private(k,j,i)
do k=1,nlevs
do j=1,ny_res
Expand All @@ -1071,6 +1108,13 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid
enddo
enddo
enddo
t2=mpi_wtime()
write(6,*)"thinkdeb255 1 dt is ",t2-t1
t1=mpi_wtime()
workvar3d=workvar3d+workinc3d
t2=mpi_wtime()
write(6,*)"thinkdeb255 2 dt is ",t2-t1

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
Expand Down Expand Up @@ -1145,7 +1189,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d)
!$omp parallel do schedule(dynamic,1) private(k,j,i,nn)
do k=1,nlevs
nn = nn_tile0
!clt nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn_tile0+(j-1)*nx_res+i
Expand All @@ -1168,21 +1212,31 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid
call write_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d)

endif

write(6,*)'thinkdeb255 dbz_ind is ',dbz_ind
if (dbz_ind > 0) then
varstrname = 'ref_f3d'
call fv3lamfile%get_idfn(varstrname,file_id,fv3filename)
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d)
!$omp parallel do schedule(dynamic,1) private(k,j,i)
do k=1,nlevs
nn = nn_tile0
!clt nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
!clt nn=nn+1
nn=nn_tile0+(j-1)*nx_res+i
workinc3d(i,j,nlevs+1-k)=vargrid(nn,levels(dbz_ind-1)+k,nb,ne)
enddo
enddo
enddo
workvar3d=workvar3d+workinc3d
!$omp parallel do schedule(dynamic,1) private(k,j,i)
do k=1,nlevs
do j=1,ny_res
do i=1,nx_res
workvar3d(i,j,k)=workvar3d(i,j,k)+workinc3d(i,j,k)
enddo
enddo
enddo

where (workvar3d < 0.0_r_kind) workvar3d = 0.0_r_kind
call write_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d)

Expand Down Expand Up @@ -1613,7 +1667,14 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, &
enddo
enddo
enddo
tvworkvar3d=workvar3d
!$omp parallel do schedule(dynamic,1) private(k,j,i)
do k=1,nlevs
do j=1,ny_res
do i=1,nx_res
tvworkvar3d(i,j,k)=workvar3d(i,j,k)
enddo
enddo
enddo
else! tsen_id >0
!$omp parallel do schedule(dynamic,1) private(k,j,i)
do k=1,nlevs
Expand Down
30 changes: 26 additions & 4 deletions src/enkf/read_fv3reg_restarts.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,8 @@ module read_fv3regional_restarts
use netcdf, only: nf90_inq_dimid,nf90_inq_varid
use netcdf, only: nf90_nowrite,nf90_inquire,nf90_inquire_dimension
use netcdf_mod, only: nc_check
use mpisetup, only: mpi_wtime
implicit none
public read_fv3_restart_data1d,read_fv3_restart_data2d
public read_fv3_restart_data3d,read_fv3_restart_data4d

Expand Down Expand Up @@ -45,15 +47,29 @@ subroutine read_fv3_restart_data2d(varname,filename,file_id,data_arr)
end subroutine read_fv3_restart_data2d

subroutine read_fv3_restart_data3d(varname,filename,file_id,data_arr)
use omp_lib
real(r_single), intent(inout), dimension(:,:,:) :: data_arr
character(len=24),parameter :: myname_ = 'read_fv3_restart_data3d'
character(len=*), intent(in) :: varname
character(len=*), intent(in) :: filename
integer(i_kind), intent(in) :: file_id
real(r_single) :: rtmp


integer(i_kind) :: var_id
integer(i_kind)::i,j,k,it
integer(i_kind):: ilow,iup, jlow,jup,klow,kup,klev


integer numThreads
real t1,t2

!$omp parallel
numThreads = omp_get_num_threads()
!$omp end parallel
write(6,*)'thinkdeb number of threads is ',numThreads


ilow=lbound(data_arr,1)
iup=ubound(data_arr,1)
jlow=lbound(data_arr,2)
Expand All @@ -62,8 +78,9 @@ subroutine read_fv3_restart_data3d(varname,filename,file_id,data_arr)
kup=ubound(data_arr,3)
call nc_check( nf90_inq_varid(file_id,trim(adjustl(varname)),var_id),&
myname_,'inq_varid '//trim(adjustl(varname))//' '//trim(filename) )
call nc_check( nf90_get_var(file_id,var_id,data_tmp),&
call nc_check( nf90_get_var(file_id,var_id,data_arr),&
myname_,'get_var '//trim(adjustl(varname))//' '//trim(filename) )
t1=mpi_wtime()
!$omp parallel do schedule(dynamic,1) private(k,j,i,klev,rtmp)
do k=klow, klow+ ceiling((kup-klow+1)/2.0)
klev=kup-k+klow
Expand All @@ -75,9 +92,14 @@ subroutine read_fv3_restart_data3d(varname,filename,file_id,data_arr)
enddo
enddo
enddo
t2=mpi_wtime()
write(6,*)'thinkdeb openmp time is ',t2-t1

! data_arr=data_arr(:,:, &
! ubound(data_arr,3):lbound(data_arr,3):-1)
t1=mpi_wtime()
data_arr=data_arr(:,:, &
ubound(data_arr,3):lbound(data_arr,3):-1)
t2=mpi_wtime()
write(6,*)'thinkdeb no-openmp time is ',t2-t1
end subroutine read_fv3_restart_data3d

subroutine read_fv3_restart_data4d(varname,filename,file_id,data_arr)
Expand All @@ -100,7 +122,7 @@ subroutine read_fv3_restart_data4d(varname,filename,file_id,data_arr)
tup=ubound(data_arr,4)
call nc_check( nf90_inq_varid(file_id,trim(adjustl(varname)),var_id),&
myname_,'inq_varid '//trim(adjustl(varname))//' '//trim(filename) )
call nc_check( nf90_get_var(file_id,var_id,data_tmp),&
call nc_check( nf90_get_var(file_id,var_id,data_arr),&
myname_,'get_var '//trim(adjustl(varname))//' '//trim(filename) )
do it=tlow, tup !normally, this dimnsions is very small or even 1,
! hence, the openmp parallelization is started from next
Expand Down
1 change: 1 addition & 0 deletions src/enkf/write_fv3reg_restarts.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ module write_fv3regional_restarts
use netcdf, only: nf90_inq_dimid,nf90_inq_varid
use netcdf, only: nf90_nowrite,nf90_inquire,nf90_inquire_dimension
use netcdf_mod, only: nc_check
implicit none
public write_fv3_restart_data1d,write_fv3_restart_data2d
public write_fv3_restart_data3d,write_fv3_restart_data4d

Expand Down

0 comments on commit ddfe627

Please sign in to comment.