diff --git a/src/gsi/anberror.f90 b/src/gsi/anberror.f90 index aef1adaeae..1acc9040c9 100644 --- a/src/gsi/anberror.f90 +++ b/src/gsi/anberror.f90 @@ -284,7 +284,7 @@ subroutine create_anberror_vars(mype) !$$$ end documentation block use fgrid2agrid_mod, only: create_fgrid2agrid - use jfunc, only: nrclen,nclen,diag_precon + use jfunc, only: nrclen,nclen use berror, only: varprd,vprecond,bnf=>nf,bnr=>nr use gridmod, only: nlat,nlon implicit none @@ -292,7 +292,7 @@ subroutine create_anberror_vars(mype) integer(i_kind),intent(in ) :: mype allocate(varprd(max(1,nrclen))) - if(diag_precon)allocate(vprecond(nclen)) + allocate(vprecond(nclen)) allocate(an_amp(max_ngauss,nvars)) an_amp=one/three @@ -388,14 +388,13 @@ subroutine destroy_anberror_vars ! !$$$ use fgrid2agrid_mod, only: destroy_fgrid2agrid - use jfunc, only: diag_precon use berror, only: vprecond implicit none deallocate(an_amp) deallocate(afact0) deallocate(qvar3d) - if(diag_precon)deallocate(vprecond) + deallocate(vprecond) call destroy_fgrid2agrid(pf2aP1) call destroy_fgrid2agrid(pf2aP2) @@ -431,7 +430,7 @@ subroutine create_anberror_vars_reg(mype) ! !$$$ use fgrid2agrid_mod, only: create_fgrid2agrid - use jfunc, only: nrclen,nclen,diag_precon + use jfunc, only: nrclen,nclen use berror, only: varprd,vprecond use gridmod, only: nlat,nlon,istart,jstart use general_commvars_mod, only: s2g_raf @@ -442,7 +441,7 @@ subroutine create_anberror_vars_reg(mype) logical regional allocate(varprd(max(1,nrclen))) - if(diag_precon)allocate(vprecond(nclen)) + allocate(vprecond(nclen)) allocate(an_amp(max_ngauss,nvars)) an_amp=one/three @@ -654,7 +653,6 @@ subroutine destroy_anberror_vars_reg !$$$ end documentation block use fgrid2agrid_mod, only: destroy_fgrid2agrid - use jfunc, only: diag_precon use berror, only: vprecond use general_sub2grid_mod, only: general_sub2grid_destroy_info implicit none @@ -662,7 +660,7 @@ subroutine destroy_anberror_vars_reg deallocate(an_amp) deallocate(afact0) deallocate(qvar3d) - if(diag_precon)deallocate(vprecond) + deallocate(vprecond) call destroy_fgrid2agrid(pf2aP1) call general_sub2grid_destroy_info(s2g_rff) !write(6,'(" FOR TEST ONLY--REMOVE THIS MESSAGE BEFORE FINAL COMMIT--SUCCESSFUL CALL TO general_sub2grid_destroy_info to remove s2g_rff in destroy_anberror_vars_reg")') diff --git a/src/gsi/anbkerror.f90 b/src/gsi/anbkerror.f90 index adcd769ddf..999a00f1ab 100644 --- a/src/gsi/anbkerror.f90 +++ b/src/gsi/anbkerror.f90 @@ -1,4 +1,4 @@ -subroutine anbkerror(gradx,grady) +subroutine anbkerror(grady) !$$$ subprogram documentation block ! . . . . ! subprogram: anbkerror apply anisotropic background error covariance @@ -33,10 +33,10 @@ subroutine anbkerror(gradx,grady) ! 2015-07-02 pondeca - update slab mode option to work with any number of control variables ! ! input argument list: -! gradx - input field +! grady - input field ! ! output -! grady - background structure * gradx +! grady - background structure * grady ! ! attributes: ! language: f90 @@ -57,7 +57,6 @@ subroutine anbkerror(gradx,grady) implicit none ! Declare passed variables - type(control_vector),intent(inout) :: gradx type(control_vector),intent(inout) :: grady ! Declare local variables @@ -82,9 +81,6 @@ subroutine anbkerror(gradx,grady) ! Initialize timer call timer_ini('anbkerror') -! Put things in grady first since operations change input variables - grady=gradx - ! Since each internal vector [step(jj)] of grad has the same structure, pointers ! are the same independent of the subwindow jj call gsi_bundlegetpointer (grady%step(1),myvnames,ipnts,istatus) diff --git a/src/gsi/berror.f90 b/src/gsi/berror.f90 index ba501c6262..71bd4f1ea8 100644 --- a/src/gsi/berror.f90 +++ b/src/gsi/berror.f90 @@ -249,7 +249,7 @@ subroutine create_berror_vars !$$$ use balmod, only: llmin,llmax use gridmod, only: nlat,nlon,lat2,lon2,nsig,nnnn1o - use jfunc, only: nrclen,nclen,diag_precon + use jfunc, only: nrclen,nclen use constants, only: zero,one implicit none @@ -293,7 +293,7 @@ subroutine create_berror_vars dssvs = zero endif allocate(varprd(nrclen)) - if(diag_precon)allocate(vprecond(nclen)) + allocate(vprecond(nclen)) allocate(inaxs(nf,nlon/8),inxrs(nlon/8,mr:nr) ) allocate(slw(ny*nx,nnnn1o),& @@ -330,7 +330,6 @@ subroutine destroy_berror_vars ! machine: ibm RS/6000 SP ! !$$$ - use jfunc, only: diag_precon implicit none if(allocated(table)) deallocate(table) deallocate(wtaxs) @@ -342,7 +341,7 @@ subroutine destroy_berror_vars if(allocated(alv)) deallocate(alv) if(allocated(dssv)) deallocate(dssv) if(allocated(dssvs)) deallocate(dssvs) - if(diag_precon)deallocate(vprecond) + deallocate(vprecond) deallocate(slw,slw1,slw2) deallocate(ii,jj,ii1,jj1,ii2,jj2) @@ -566,7 +565,7 @@ subroutine pcinfo use kinds, only: r_kind,i_kind use radinfo, only: ostats,rstats,varA,jpch_rad,npred,newpc4pred use aircraftinfo, only: aircraft_t_bc_pof,aircraft_t_bc,ntail,npredt,ostats_t,rstats_t,varA_t - use jfunc, only: nclen,nrclen,diag_precon,step_start,ntclen + use jfunc, only: nclen,nrclen,step_start,ntclen,diag_precon use constants, only: zero,one implicit none @@ -580,10 +579,10 @@ subroutine pcinfo ! Only diagonal elements are considered ! set a coeff. factor for variances of control variables - if(diag_precon)then - lfact=step_start - vprecond=lfact + lfact=step_start + vprecond=lfact + if(diag_precon)then if(newpc4pred)then ! for radiance bias predictor coeff. nclen1=nclen-nrclen @@ -591,6 +590,7 @@ subroutine pcinfo do i=1,jpch_rad do j=1,npred ii=ii+1 +! if (ostats(i)>zero) vprecond(nclen1+ii)=vprecond(nclen1+ii)/(one+rstats(j,i)*varprd(ii)) if (ostats(i)>zero) vprecond(nclen1+ii)=one/(one+rstats(j,i)*varprd(ii)) if (ostats(i)>20.0_r_kind) then if (rstats(j,i)>zero) then @@ -615,6 +615,7 @@ subroutine pcinfo if (aircraft_t_bc_pof) obs_count = ostats_t(j,i) if (aircraft_t_bc) obs_count = ostats_t(1,i) +! if (obs_count>zero) vprecond(nclen1+ii)=vprecond(nclen1+ii)/(one+rstats_t(j,i)*varprd(jj)) if (obs_count>zero) vprecond(nclen1+ii)=one/(one+rstats_t(j,i)*varprd(jj)) if (obs_count>3.0_r_kind) then varA_t(j,i)=one/(one/varprd(jj)+rstats_t(j,i)) @@ -917,7 +918,7 @@ subroutine create_berror_vars_reg use constants, only: zero use balmod, only: llmin,llmax use gridmod, only: nlat,nlon,nsig,nnnn1o,lat2,lon2 - use jfunc, only: nrclen,nclen,diag_precon + use jfunc, only: nrclen,nclen implicit none nx=nlon @@ -940,7 +941,7 @@ subroutine create_berror_vars_reg endif allocate(varprd(max(1,nrclen) ) ) - if(diag_precon)allocate(vprecond(nclen)) + allocate(vprecond(nclen)) allocate(slw(ny*nx,nnnn1o) ) allocate(ii(ny,nx,3,nnnn1o),jj(ny,nx,3,nnnn1o) ) @@ -973,7 +974,6 @@ subroutine destroy_berror_vars_reg ! machine: ibm RS/6000 SP ! !$$$ - use jfunc, only:diag_precon implicit none deallocate(be,qvar3d) @@ -984,7 +984,7 @@ subroutine destroy_berror_vars_reg deallocate(ii,jj) deallocate(slw) deallocate(varprd) - if(diag_precon)deallocate(vprecond) + deallocate(vprecond) return end subroutine destroy_berror_vars_reg diff --git a/src/gsi/bicg.f90 b/src/gsi/bicg.f90 index 4fc7dbe863..6eb2f78905 100644 --- a/src/gsi/bicg.f90 +++ b/src/gsi/bicg.f90 @@ -31,8 +31,7 @@ subroutine bicg() use kinds, only: r_kind,i_kind,r_quad use gsi_4dvar, only: l4dvar, & ladtest, lgrtest, lanczosave, ltcost, nwrvecs -use jfunc, only: jiter,miter,niter,xhatsave,yhatsave,jiterstart, & - diag_precon +use jfunc, only: jiter,miter,niter,xhatsave,yhatsave,jiterstart use constants, only: zero,tiny_r_kind use mpimod, only: mype use obs_sensitivity, only: lobsensmin, lobsensfc, lobsensincr, & @@ -99,23 +98,24 @@ subroutine bicg() if(LMPCGL) then call pcgprecond(gradx,grady) else - call bkerror(gradx,grady) + grady=gradx + call bkerror(grady) ! If hybrid ensemble run, then multiply ensemble control variable a_en ! by its localization correlation if(l_hyb_ens) then if(aniso_a_en) then - ! call anbkerror_a_en(gradx,grady) ! not available yet + ! call anbkerror_a_en(grady) ! not available yet write(6,*)' ANBKERROR_A_EN not written yet, program stops' call stop2 (999) else - call bkerror_a_en(gradx,grady) + call bkerror_a_en(grady) end if end if ! Add potential additional preconditioner - if(diag_precon) call precond(grady) + call precond(grady) endif zg0=dot_product(gradx,grady,r_quad) @@ -188,24 +188,25 @@ subroutine bicg() else ! not sensitivity run - call bkerror(gradf,grads) + grads=gradf + call bkerror(grads) ! If hybrid ensemble run, then multiply ensemble control variable a_en ! by its localization correlation if(l_hyb_ens) then if(aniso_a_en) then -! call anbkerror_a_en(gradf,grads) ! not available yet +! call anbkerror_a_en(grads) ! not available yet write(6,*)' ANBKERROR_A_EN not written yet, program stops' stop else - call bkerror_a_en(gradf,grads) + call bkerror_a_en(grads) end if end if ! Add potential additional preconditioner - if(diag_precon) call precond(grads) + call precond(grads) ! Update xhatsave do ii=1,xhat%lencv diff --git a/src/gsi/bicglanczos.F90 b/src/gsi/bicglanczos.F90 index ea0d89ddd0..1914b0214d 100755 --- a/src/gsi/bicglanczos.F90 +++ b/src/gsi/bicglanczos.F90 @@ -66,7 +66,7 @@ module bicglanczos use gsi_bundlemod, only: assignment(=) use mpimod , only : mpi_comm_world use mpimod, only: mype -use jfunc , only : iter, jiter, diag_precon,step_start +use jfunc , only : iter, jiter use gsi_4dvar, only : nwrvecs,l4dvar,lanczosave use gsi_4dvar, only : nsubwin, nobs_bins use hybrid_ensemble_parameters,only : l_hyb_ens,aniso_a_en @@ -246,7 +246,7 @@ subroutine pcglanczos(xhat,yhat,pcost,gradx,grady,preduc,kmaxit,lsavevecs) call allocate_cv(dirx) call allocate_cv(diry) if(nprt>=1.and.ltcost_) call allocate_cv(gradf) -if(diag_precon) call allocate_cv(dirw) +call allocate_cv(dirw) !--- 'zeta' is an upper bound on the relative error of the gradient. @@ -256,7 +256,7 @@ subroutine pcglanczos(xhat,yhat,pcost,gradx,grady,preduc,kmaxit,lsavevecs) ilen=xhat%lencv -if(diag_precon) dirw=zero +dirw=zero !$omp parallel do do jj=1,ilen @@ -266,12 +266,10 @@ subroutine pcglanczos(xhat,yhat,pcost,gradx,grady,preduc,kmaxit,lsavevecs) end do !$omp end parallel do -if(diag_precon) then - do jj=1,ilen - dirw%values(jj)=diry%values(jj) - end do - call precond(diry) -end if +do jj=1,ilen + dirw%values(jj)=diry%values(jj) +end do +call precond(diry) if(LMPCGL) then dirx=zero @@ -284,7 +282,7 @@ subroutine pcglanczos(xhat,yhat,pcost,gradx,grady,preduc,kmaxit,lsavevecs) zg0=dot_product(gradx,grady,r_quad) if(zg0=1.and.ltcost_) call deallocate_cv(gradf) call deallocate_cv(diry) call deallocate_cv(dirx) @@ -370,24 +368,25 @@ subroutine pcglanczos(xhat,yhat,pcost,gradx,grady,preduc,kmaxit,lsavevecs) if(LMPCGL) then call pcgprecond(gradx,grady) else - call bkerror(gradx,grady) + grady=gradx + call bkerror(grady) ! If hybrid ensemble run, then multiply ensemble control variable a_en ! by its localization correlation if(l_hyb_ens) then if(aniso_a_en) then - ! call anbkerror_a_en(gradx,grady) ! not available yet + ! call anbkerror_a_en(grady) ! not available yet write(6,*)' ANBKERROR_A_EN not written yet, program stops' stop else - call bkerror_a_en(gradx,grady) + call bkerror_a_en(grady) end if end if endif ! Add potential additional preconditioner - if(diag_precon) call precond(grady) + call precond(grady) ! Second re-orthogonalization @@ -416,21 +415,17 @@ subroutine pcglanczos(xhat,yhat,pcost,gradx,grady,preduc,kmaxit,lsavevecs) endif ! Update search direction - if(diag_precon) then - do jj=1,ilen - diry%values(jj)=dirw%values(jj) - enddo - end if + do jj=1,ilen + diry%values(jj)=dirw%values(jj) + enddo do jj=1,ilen dirx%values(jj)=-grady%values(jj)+beta(iter)*dirx%values(jj) diry%values(jj)=-gradx%values(jj)+beta(iter)*diry%values(jj) end do - if(diag_precon) then - do jj=1,ilen - dirw%values(jj)=diry%values(jj) - end do - call precond(diry) - end if + do jj=1,ilen + dirw%values(jj)=diry%values(jj) + end do + call precond(diry) ! Diagnostics if(zgk < zero) then @@ -682,7 +677,7 @@ subroutine pcglanczos(xhat,yhat,pcost,gradx,grady,preduc,kmaxit,lsavevecs) call deallocate_cv(dirx) call deallocate_cv(diry) if(nprt>=1.and.ltcost_) call deallocate_cv(gradf) -if(diag_precon) call deallocate_cv(dirw) +call deallocate_cv(dirw) call inquire_cv @@ -905,18 +900,19 @@ subroutine pcgprecond(xcvx,ycvx) enddo !Apply B -call bkerror(xcvx,ycvx) +ycvx=xcvx +call bkerror(ycvx) ! If hybrid ensemble run, then multiply ensemble control variable a_en ! by its localization correlation if(l_hyb_ens) then if(aniso_a_en) then -! call anbkerror_a_en(xcvx,ycvx) ! not available yet +! call anbkerror_a_en(ycvx) ! not available yet write(6,*)' ANBKERROR_A_EN not written yet, program stops' call stop2(999) else - call bkerror_a_en(xcvx,ycvx) + call bkerror_a_en(ycvx) end if end if diff --git a/src/gsi/bkerror.f90 b/src/gsi/bkerror.f90 index 261a4c66ca..b3a0140691 100644 --- a/src/gsi/bkerror.f90 +++ b/src/gsi/bkerror.f90 @@ -1,4 +1,4 @@ -subroutine bkerror(gradx,grady) +subroutine bkerror(grady) !$$$ subprogram documentation block ! . . . . @@ -41,10 +41,10 @@ subroutine bkerror(gradx,grady) ! this is a hybrid ensemble run. ! ! input argument list: -! gradx - input field +! grady - input field ! ! output -! grady - background structure * gradx +! grady - background structure * grady ! ! attributes: ! language: f90 @@ -55,7 +55,6 @@ subroutine bkerror(gradx,grady) use berror, only: varprd,fpsproj,fut2ps use balmod, only: balance,tbalance use gsi_4dvar, only: nsubwin, lsqrtb - use gridmod, only: nlat,nlon,periodic use jfunc, only: nsclen,npclen,ntclen use jfunc, only: set_sqrt_2dsize use constants, only: zero @@ -63,15 +62,11 @@ subroutine bkerror(gradx,grady) use control_vectors, only: mvars,nrf,nrf_var,nrf_3d use timermod, only: timer_ini,timer_fnl use gsi_bundlemod, only: gsi_bundlegetpointer,gsi_bundlemerge,gsi_bundle,gsi_bundledup,gsi_bundledestroy - use general_sub2grid_mod, only: general_sub2grid,general_grid2sub -! use general_commvars_mod, only: s2g_raf - use general_commvars_mod, only: s2g_cv use hybrid_ensemble_isotropic, only: sqrt_beta_s_mult use hybrid_ensemble_parameters, only: l_hyb_ens implicit none ! Declare passed variables - type(control_vector),intent(inout) :: gradx type(control_vector),intent(inout) :: grady ! Declare local variables @@ -79,7 +74,6 @@ subroutine bkerror(gradx,grady) integer(i_kind) i_t,i_p,i_st,i_vp integer(i_kind) ipnts(4),istatus ! integer(i_kind) nval_lenz,ndim2d - real(r_kind),dimension(nlat*nlon*s2g_cv%nlevs_alloc)::workcv real(r_kind),pointer,dimension(:,:,:):: p_t =>NULL() real(r_kind),pointer,dimension(:,:,:):: p_st =>NULL() real(r_kind),pointer,dimension(:,:,:):: p_vp =>NULL() @@ -100,23 +94,6 @@ subroutine bkerror(gradx,grady) ! Initialize timer call timer_ini('bkerror') -! If dealing with periodic (sub)domain, gather full domain grids, -! account for periodicity, and redistribute to subdomains. This -! only needs to be done when running with a single mpi task and -! then only for array gradx. - if (periodic) then - do ii=1,nsubwin - call general_sub2grid(s2g_cv,gradx%step(ii)%values,workcv) - call general_grid2sub(s2g_cv,workcv,gradx%step(ii)%values) - end do - endif - -! Put things in grady first since operations change input variables - grady=gradx - -! if ensemble run, multiply by sqrt_beta_s - if(l_hyb_ens) call sqrt_beta_s_mult(grady) - ! Only need to get pointer for ii=1 - all other are the same call gsi_bundlegetpointer ( grady%step(1), (/'t ','sf','vp','ps'/), & ipnts, istatus ) @@ -126,6 +103,9 @@ subroutine bkerror(gradx,grady) i_p = ipnts(4) dobal = i_t>0.and.i_p>0.and.i_st>0.and.i_vp>0 +! if ensemble run, multiply by sqrt_beta_s + if(l_hyb_ens) call sqrt_beta_s_mult(grady) + ! Loop on control steps do ii=1,nsubwin @@ -186,7 +166,12 @@ subroutine bkerror(gradx,grady) call stop2(999) endif +! if ensemble run, multiply by sqrt_beta_s + if(l_hyb_ens) call sqrt_beta_s_mult(grady) + + ! Take care of background error for bias correction terms + do i=1,nsclen grady%predr(i)=grady%predr(i)*varprd(i) end do @@ -199,9 +184,6 @@ subroutine bkerror(gradx,grady) end do end if -! if ensemble run, multiply by sqrt_beta_s - if(l_hyb_ens) call sqrt_beta_s_mult(grady) - ! Finalize timer call timer_fnl('bkerror') diff --git a/src/gsi/bkgvar.f90 b/src/gsi/bkgvar.f90 index f2f5e449a5..705212eeb5 100644 --- a/src/gsi/bkgvar.f90 +++ b/src/gsi/bkgvar.f90 @@ -83,7 +83,7 @@ subroutine bkgvar(cvec,iflg) real(r_kind),pointer,dimension(:,:) ::ptrsst=>NULL() real(r_kind),pointer,dimension(:,:) ::ptrstl=>NULL() real(r_kind),pointer,dimension(:,:) ::ptrsti=>NULL() - real(r_kind),dimension(lat2,lon2) :: sst,stl,sti + real(r_kind),dimension(lat2,lon2) :: sst ! Multipy by variances !$omp parallel do schedule(dynamic,1) private(n,k,i,j,ptr3d,istatus) @@ -105,15 +105,17 @@ subroutine bkgvar(cvec,iflg) call gsi_bundlegetpointer(cvec,'sst',ptrsst,istatus) call gsi_bundlegetpointer(cvec,'stl',ptrstl,istatus) call gsi_bundlegetpointer(cvec,'sti',ptrsti,istatus) - sst=zero - stl=zero - sti=zero + if(iflg == 0)then + if(i_sst > 0)sst=zero + if(i_stl > 0)ptrstl=zero + if(i_sti > 0)ptrsti=zero + end if ! Surface fields ! !$omp parallel do schedule(dynamic,1) private(n,i,j,ptr2d,istatus) do n=1,cvec%n2d - call gsi_bundlegetpointer(cvec,cvec%r2(n)%shortname,ptr2d,istatus) if(n/=i_sst.and.n/=i_stl.and.n/=i_sti) then + call gsi_bundlegetpointer(cvec,cvec%r2(n)%shortname,ptr2d,istatus) do i=1,lon2 do j=1,lat2 ptr2d(j,i)=ptr2d(j,i)*dssvs(j,i,n) @@ -130,13 +132,13 @@ subroutine bkgvar(cvec,iflg) elseif(n==i_stl) then do i=1,lon2 do j=1,lat2 - if(isli2(j,i) == 1) stl(j,i)=ptrsst(j,i)*dssvs(j,i,n) + if(isli2(j,i) == 1) ptrstl(j,i)=ptrsst(j,i)*dssvs(j,i,n) end do end do elseif(n==i_sti) then do i=1,lon2 do j=1,lat2 - if(isli2(j,i) == 2) sti(j,i)=ptrsst(j,i)*dssvs(j,i,n) + if(isli2(j,i) == 2) ptrsti(j,i)=ptrsst(j,i)*dssvs(j,i,n) end do end do end if @@ -144,19 +146,19 @@ subroutine bkgvar(cvec,iflg) if(n==i_sst) then do i=1,lon2 do j=1,lat2 - if(isli2(j,i)/=1.and.isli2(j,i)/=2) sst(j,i)=ptrsst(j,i)*dssvs(j,i,n) + if(isli2(j,i)/=1.and.isli2(j,i)/=2) ptrsst(j,i)=ptrsst(j,i)*dssvs(j,i,n) end do end do elseif(n==i_stl) then do i=1,lon2 do j=1,lat2 - if(isli2(j,i) == 1) sst(j,i)=ptrstl(j,i)*dssvs(j,i,n) + if(isli2(j,i) == 1) ptrsst(j,i)=ptrstl(j,i)*dssvs(j,i,n) end do end do elseif(n==i_sti) then do i=1,lon2 do j=1,lat2 - if(isli2(j,i) == 2) sst(j,i)=ptrsti(j,i)*dssvs(j,i,n) + if(isli2(j,i) == 2) ptrsst(j,i)=ptrsti(j,i)*dssvs(j,i,n) end do end do end if @@ -164,15 +166,8 @@ subroutine bkgvar(cvec,iflg) end if end do + if(iflg == 0 .and. i_sst>0) ptrsst=sst - if(iflg==0) then - if(i_sst>0) ptrsst=sst - if(i_stl>0) ptrstl=stl - if(i_sti>0) ptrsti=sti - else - if(i_sst>0) ptrsst=sst -! ignore contents of ptrstl,ptrsti - end if return end subroutine bkgvar diff --git a/src/gsi/bkgvar_rewgt.f90 b/src/gsi/bkgvar_rewgt.f90 index 5c1174dc6e..ca82882af3 100644 --- a/src/gsi/bkgvar_rewgt.f90 +++ b/src/gsi/bkgvar_rewgt.f90 @@ -222,15 +222,15 @@ subroutine bkgvar_rewgt(sfvar,vpvar,tvar,psvar,mype) do k=1,nsig do j=1,lon2 do i=1,lat2 - delpsi(i,j,k)=sqrt( delpsi(i,j,k)**two ) - delchi(i,j,k)=sqrt( delchi(i,j,k)**two ) - deltv (i,j,k)=sqrt( deltv (i,j,k)**two ) + delpsi(i,j,k)=abs( delpsi(i,j,k) ) + delchi(i,j,k)=abs( delchi(i,j,k) ) + deltv (i,j,k)=abs( deltv (i,j,k) ) end do end do end do do j=1,lon2 do i=1,lat2 - delps(i,j)=sqrt( delps(i,j)**two ) + delps(i,j)=abs( delps(i,j) ) end do end do diff --git a/src/gsi/control2state.f90 b/src/gsi/control2state.f90 index c652d1574d..d32e37b13a 100644 --- a/src/gsi/control2state.f90 +++ b/src/gsi/control2state.f90 @@ -237,7 +237,7 @@ subroutine control2state(xhat,sval,bval) call general_grid2sub(s2g_cv,hwork,wbundle%values) end if -!$omp parallel sections private(istatus,ii,ic,id,sv_u,sv_v,sv_prse,sv_q,sv_tsen,uland,vland,uwter,vwter) +!$omp parallel sections private(istatus,ii,ic,id,uland,vland,uwter,vwter) !$omp section @@ -289,6 +289,8 @@ subroutine control2state(xhat,sval,bval) call gsi_bundlegetpointer (wbundle,'t' ,cv_t, istatus) call gsi_bundlegetpointer (wbundle,'q' ,cv_rh ,istatus) +! Copy other variables + call gsi_bundlegetvar ( wbundle, 't' , sv_tv, istatus ) ! Get 3d pressure if(do_getprs_tl) call getprs_tl(cv_ps,cv_t,sv_prse) @@ -298,8 +300,6 @@ subroutine control2state(xhat,sval,bval) ! Calculate sensible temperature if(do_tv_to_tsen) call tv_to_tsen(cv_t,sv_q,sv_tsen) -! Copy other variables - call gsi_bundlegetvar ( wbundle, 't' , sv_tv, istatus ) if (do_cw_to_hydro .and. .not.do_cw_to_hydro_hwrf) then ! Case when cloud-vars do not map one-to-one (cv-to-sv) @@ -323,10 +323,10 @@ subroutine control2state(xhat,sval,bval) endif enddo end if - -!$omp section call gsi_bundlegetpointer (sval(jj),'ps' ,sv_ps, istatus) call gsi_bundlegetvar ( wbundle, 'ps' , sv_ps, istatus ) + +!$omp section call gsi_bundlegetpointer (sval(jj),'sst' ,sv_sst, istatus) call gsi_bundlegetvar ( wbundle, 'sst', sv_sst, istatus ) call gsi_bundlegetpointer (sval(jj),'oz' ,sv_oz , istatus_oz) diff --git a/src/gsi/convinfo.f90 b/src/gsi/convinfo.f90 index 9b19a29b06..f6ba29c2e1 100755 --- a/src/gsi/convinfo.f90 +++ b/src/gsi/convinfo.f90 @@ -250,9 +250,11 @@ subroutine convinfo_read end if if(cflg == '!')cycle read(crecord,*)ictypet,icsubtypet,icuset - if (mype==0 .and. icuset < use_limit) write(6, *) & - 'line ignored in convinfo due to use flag ',cflg,iotype,ictypet,icsubtypet,icuset - if(icuset < use_limit)cycle + if(icuset < use_limit)then + if (mype==0) write(6, *) 'line ignored in convinfo due to use flag ',& + cflg,iotype,ictypet,icsubtypet,icuset + cycle + end if nc=nc+1 ioctype(nc)=iotype !otype type isub iuse twindow numgrp ngroup nmiter gross ermax ermin var_b var_pg ithin rmesh pmesh npred pmot ptime diff --git a/src/gsi/cplr_gfs_ensmod.f90 b/src/gsi/cplr_gfs_ensmod.f90 index d9614b83aa..4b98d904b6 100644 --- a/src/gsi/cplr_gfs_ensmod.f90 +++ b/src/gsi/cplr_gfs_ensmod.f90 @@ -104,26 +104,16 @@ subroutine get_gfs_Nens(this,grd,members,ntindex,atm_bundle,iret) ! Declare internal variables character(len=*),parameter :: myname_='get_user_ens_gfs' - real(r_single),allocatable,dimension(:,:,:,:) :: en_loc3 - integer(i_kind) :: m_cvars2d(nc2d),m_cvars3d(nc3d) integer(i_kind) :: n - real(r_kind),allocatable,dimension(:) :: clons,slons associate( this => this ) ! eliminates warning for unused dummy argument needed for binding end associate if ( (use_gfs_nemsio .or. use_gfs_ncio) .and. ens_fast_read ) then - allocate(en_loc3(grd_ens%lat2,grd_ens%lon2,nc2d+nc3d*grd_ens%nsig,members)) - allocate(clons(grd_ens%nlon),slons(grd_ens%nlon)) - call get_user_ens_gfs_fastread_(ntindex,en_loc3,m_cvars2d,m_cvars3d, & - grd_ens%lat2,grd_ens%lon2,grd_ens%nsig, & - nc2d,nc3d,members,iret,clons,slons) - do n=1,members - call move2bundle_(grd,en_loc3(:,:,:,n),atm_bundle(n), & - m_cvars2d,m_cvars3d,iret,clons,slons) - end do - deallocate(en_loc3,clons,slons) + call get_user_ens_gfs_fastread_(ntindex,atm_bundle, & + grd_ens%lat2,grd_ens%lon2, & + nc2d,nc3d,iret,grd) else do n = 1,members call get_gfs_ens(this,grd,n,ntindex,atm_bundle(n),iret) @@ -134,8 +124,8 @@ subroutine get_gfs_Nens(this,grd,members,ntindex,atm_bundle,iret) end subroutine get_gfs_Nens -subroutine get_user_ens_gfs_fastread_(ntindex,en_loc3,m_cvars2d,m_cvars3d, & - lat2in,lon2in,nsigin,nc2din,nc3din,n_ensin,iret,clons,slons) +subroutine get_user_ens_gfs_fastread_(ntindex,atm_bundle, & + lat2in,lon2in,nc2din,nc3din,iret,grd) !$$$ subprogram documentation block ! . . . . ! subprogram: get_user_ens_gfs_fastread_ @@ -178,8 +168,9 @@ subroutine get_user_ens_gfs_fastread_(ntindex,en_loc3,m_cvars2d,m_cvars3d, & use mpimod, only: mpi_comm_world,ierror,mpi_real8,mpi_integer4,mpi_max use kinds, only: i_kind,r_single,r_kind use constants, only: zero - use general_sub2grid_mod, only: sub2grid_info + use general_sub2grid_mod, only: sub2grid_info,general_sub2grid_destroy_info use gsi_4dvar, only: ens_fhrlevs + use gsi_bundlemod, only: gsi_bundle use hybrid_ensemble_parameters, only: n_ens,grd_ens use hybrid_ensemble_parameters, only: ensemble_path use control_vectors, only: nc2d,nc3d @@ -192,11 +183,10 @@ subroutine get_user_ens_gfs_fastread_(ntindex,en_loc3,m_cvars2d,m_cvars3d, & ! Declare passed variables integer(i_kind), intent(in ) :: ntindex - real(r_single), intent(inout) :: en_loc3(lat2in,lon2in,nc2din+nc3din*nsigin,n_ensin) - integer(i_kind), intent(inout) :: m_cvars2d(nc2din),m_cvars3d(nc3din) - integer(i_kind), intent(in ) :: lat2in,lon2in,nsigin,nc2din,nc3din,n_ensin + integer(i_kind), intent(in ) :: lat2in,lon2in,nc2din,nc3din integer(i_kind), intent( out) :: iret - real(r_kind), intent(inout) :: clons(grd_ens%nlon),slons(grd_ens%nlon) + type(sub2grid_info), intent(in ) :: grd + type(gsi_bundle), intent(inout) :: atm_bundle(:) ! Declare internal variables @@ -214,27 +204,21 @@ subroutine get_user_ens_gfs_fastread_(ntindex,en_loc3,m_cvars2d,m_cvars3d, & integer(i_kind) :: nlon,nlat,nsig type(genex_info) :: s_a2b real(r_single),allocatable,dimension(:,:,:,:) :: en_full,en_loc + real(r_kind),allocatable,dimension(:,:,:) :: en_loc3 integer(i_kind),allocatable,dimension(:) :: m_cvars2dw,m_cvars3dw - integer(i_kind) base_pe,base_pe0 + integer(i_kind) :: m_cvars2d(nc2d),m_cvars3d(nc3d) + type(sub2grid_info) :: grd3d - iret = 0 nlat=grd_ens%nlat nlon=grd_ens%nlon nsig=grd_ens%nsig - ! write out contents of cvars2d, cvars3d - - !if (mype == 0 ) then - ! write(6,*) ' in get_user_ens_fastread_,cvars2d=',(trim(cvars2d(i)),i=1,2) - ! write(6,*) ' in get_user_ens_fastread_,cvars3d=',(trim(cvars3d(i)),i=1,6) - !endif - ! set up partition of available processors for parallel read if ( n_ens > npe ) & call die(myname_, ': ***ERROR*** CANNOT READ ENSEMBLE n_ens > npe, increase npe >= n_ens', 99) - call ens_io_partition_(n_ens,io_pe,n_io_pe_s,n_io_pe_e,n_io_pe_em,i_ens) + call ens_io_partition_(n_ens,ntindex,io_pe,n_io_pe_s,n_io_pe_e,n_io_pe_em,i_ens) ! setup communicator for scatter to subdomains: @@ -253,9 +237,9 @@ subroutine get_user_ens_gfs_fastread_(ntindex,en_loc3,m_cvars2d,m_cvars3d, & n2d=nc3d*grd_ens%nsig+nc2d ias=1 ; iae=0 ; jas=1 ; jae=0 ; kas=1 ; kae=0 ; mas=1 ; mae=0 if(mype==io_pe) then - ias=1 ; iae=nlat - jas=1 ; jae=nlon - kas=1 ; kae=n2d + iae=nlat + jae=nlon + kae=n2d mas=n_io_pe_s ; mae=n_io_pe_em endif iasm=ias ; iaem=iae ; jasm=jas ; jaem=jae ; kasm=kas ; kaem=kae ; masm=mas ; maem=mae @@ -273,52 +257,44 @@ subroutine get_user_ens_gfs_fastread_(ntindex,en_loc3,m_cvars2d,m_cvars3d, & iasm,iaem,jasm,jaem,kasm,kaem,masm,maem, & ibsm,ibem,jbsm,jbem,kbsm,kbem,mbsm,mbem) -!! read ensembles - - allocate(en_full(iasm:iaemz,jasm:jaemz,kasm:kaemz,masm:maemz)) - write(filename,22) trim(adjustl(ensemble_path)),ens_fhrlevs(ntindex),mas 22 format(a,'sigf',i2.2,'_ens_mem',i3.3) - if (cnvw_option) then - write(filenamesfc,23) trim(adjustl(ensemble_path)),ens_fhrlevs(ntindex),mas -23 format(a,'sfcf',i2.2,'_ens_mem',i3.3) - end if - allocate(m_cvars2dw(nc2din),m_cvars3dw(nc3din)) m_cvars2dw=-999 m_cvars3dw=-999 + allocate(en_full(iasm:iaemz,jasm:jaemz,kasm:kaemz,masm:maemz)) + +!! read ensembles + if ( mas == mae ) then if ( use_gfs_nemsio ) then if (cnvw_option) then + write(filenamesfc,23) trim(adjustl(ensemble_path)),ens_fhrlevs(ntindex),mas +23 format(a,'sfcf',i2.2,'_ens_mem',i3.3) call parallel_read_nemsio_state_(en_full,m_cvars2dw,m_cvars3dw,nlon,nlat,nsig, & ias,jas,mas, & iasm,iaemz,jasm,jaemz,kasm,kaemz,masm,maemz, & - filename,.true.,clons,slons,filenamesfc) + filename,.true.,filenamesfc) else call parallel_read_nemsio_state_(en_full,m_cvars2dw,m_cvars3dw,nlon,nlat,nsig, & ias,jas,mas, & iasm,iaemz,jasm,jaemz,kasm,kaemz,masm,maemz, & - filename,.true.,clons,slons) + filename,.true.) end if else call parallel_read_gfsnc_state_(en_full,m_cvars2dw,m_cvars3dw,nlon,nlat,nsig, & ias,jas,mas, & iasm,iaemz,jasm,jaemz,kasm,kaemz,masm,maemz, & - filename,.true.,clons,slons) + filename) end if end if - base_pe0=-999 - if ( mas == 1 .and. mae == 1 ) base_pe0=mype - - call mpi_allreduce(base_pe0,base_pe,1,mpi_integer4,mpi_max,mpi_comm_world,ierror) - call mpi_bcast(clons,grd_ens%nlon,mpi_real8,base_pe,mpi_comm_world,ierror) - call mpi_bcast(slons,grd_ens%nlon,mpi_real8,base_pe,mpi_comm_world,ierror) call mpi_allreduce(m_cvars2dw,m_cvars2d,nc2d,mpi_integer4,mpi_max,mpi_comm_world,ierror) call mpi_allreduce(m_cvars3dw,m_cvars3d,nc3d,mpi_integer4,mpi_max,mpi_comm_world,ierror) + deallocate(m_cvars2dw,m_cvars3dw) ! scatter to subdomains: allocate(en_loc(ibsm:ibemz,jbsm:jbemz,kbsm:kbemz,mbsm:mbemz)) @@ -329,9 +305,12 @@ subroutine get_user_ens_gfs_fastread_(ntindex,en_loc3,m_cvars2d,m_cvars3d, & deallocate(en_full) call genex_destroy_info(s_a2b) ! check on actual routine name -! transfer en_loc to en_loc3 +! transfer en_loc to en_loc3 then to atm_bundle + + allocate(en_loc3(lat2in,lon2in,nc2d+nc3d*nsig)) -! Look to thread here OMP + iret = 0 + call create_grd23d_(grd3d,nc2d+nc3d*grd%nsig) do n=1,n_ens do k=1,nc2d+nc3d*nsig jj=0 @@ -340,15 +319,19 @@ subroutine get_user_ens_gfs_fastread_(ntindex,en_loc3,m_cvars2d,m_cvars3d, & ii=0 do i=ibsm,ibem ii=ii+1 - en_loc3(ii,jj,k,n)=en_loc(i,j,k,n) + en_loc3(ii,jj,k)=en_loc(i,j,k,n) enddo enddo enddo + call move2bundle_(grd3d,en_loc3,atm_bundle(n),m_cvars2d,m_cvars3d,iret) enddo + call general_sub2grid_destroy_info(grd3d,grd) + deallocate(en_loc,en_loc3) + end subroutine get_user_ens_gfs_fastread_ -subroutine move2bundle_(grd,en_loc3,atm_bundle,m_cvars2d,m_cvars3d,iret,clons,slons) +subroutine move2bundle_(grd3d,en_loc3,atm_bundle,m_cvars2d,m_cvars3d,iret) !$$$ subprogram documentation block ! . . . . @@ -380,10 +363,10 @@ subroutine move2bundle_(grd,en_loc3,atm_bundle,m_cvars2d,m_cvars3d,iret,clons,sl use kinds, only: i_kind,r_kind,r_single use constants, only: zero,one,two,fv - use general_sub2grid_mod, only: sub2grid_info,general_sub2grid_destroy_info + use general_sub2grid_mod, only: sub2grid_info use hybrid_ensemble_parameters, only: en_perts use gsi_bundlemod, only: gsi_bundle - use gsi_bundlemod, only: gsi_bundlegetpointer,gsi_bundleputvar + use gsi_bundlemod, only: gsi_bundlegetpointer use gsi_bundlemod, only : assignment(=) use control_vectors, only: cvars2d,cvars3d,nc2d,nc3d use mpeu_util, only: getindex @@ -391,35 +374,28 @@ subroutine move2bundle_(grd,en_loc3,atm_bundle,m_cvars2d,m_cvars3d,iret,clons,sl implicit none ! Declare passed variables - type(sub2grid_info), intent(in ) :: grd - real(r_single), intent(in ) :: en_loc3(grd%lat2,grd%lon2,nc2d+nc3d*grd%nsig) + type(sub2grid_info), intent(in ) :: grd3d + real(r_kind), intent(inout) :: en_loc3(grd3d%lat2,grd3d%lon2,nc2d+nc3d*grd3d%nsig) type(gsi_bundle), intent(inout) :: atm_bundle integer(i_kind), intent(in ) :: m_cvars2d(nc2d),m_cvars3d(nc3d) integer(i_kind), intent( out) :: iret - real(r_kind), intent(in ) :: clons(grd%nlon),slons(grd%nlon) ! Declare internal variables character(len=*),parameter :: myname_='move2bundle_' character(len=70) :: filename integer(i_kind) :: ierr - integer(i_kind) :: im,jm,km,m,k + integer(i_kind) :: km,m integer(i_kind) :: icw,iql,iqi,iqr,iqs,iqg real(r_kind),pointer,dimension(:,:) :: ps !real(r_kind),pointer,dimension(:,:) :: sst real(r_kind),pointer,dimension(:,:,:) :: u,v,tv,q,oz,cwmr real(r_kind),pointer,dimension(:,:,:) :: qlmr,qimr,qrmr,qsmr,qgmr - real(r_single),allocatable,dimension(:,:) :: scr2 - real(r_single),allocatable,dimension(:,:,:) :: scr3 - type(sub2grid_info) :: grd2d,grd3d real(r_kind),parameter :: r0_001 = 0.001_r_kind - im = en_perts(1,1)%grid%im - jm = en_perts(1,1)%grid%jm - km = en_perts(1,1)%grid%km - allocate(scr2(im,jm)) - allocate(scr3(im,jm,km)) +!--- now update halo values of all variables using general_sub2grid + call update_halos_(grd3d,en_loc3) ! Check hydrometeors in control variables icw=getindex(cvars3d,'cw') @@ -457,81 +433,44 @@ subroutine move2bundle_(grd,en_loc3,atm_bundle,m_cvars2d,m_cvars3d,iret,clons,sl return endif + do m=1,nc2d - scr2(:,:)=en_loc3(:,:,m_cvars2d(m)) - if(trim(cvars2d(m))=='ps') ps=scr2 - ! if(trim(cvars2d(m))=='sst') sst=scr2 ! no sst for now +! convert ps from Pa to cb + if(trim(cvars2d(m))=='ps') ps=r0_001*en_loc3(:,:,m_cvars2d(m)) +! if(trim(cvars2d(m))=='sst') sst=en_loc3(:,:,m_cvars2d(m)) !no sst for now enddo + + km = en_perts(1,1)%grid%km +!$omp parallel do schedule(dynamic,1) private(m) do m=1,nc3d - do k=1,km - scr3(:,:,k)=en_loc3(:,:,m_cvars3d(m)+k-1) - enddo - if(trim(cvars3d(m))=='sf') u = scr3 - if(trim(cvars3d(m))=='vp') v = scr3 - if(trim(cvars3d(m))=='t') tv = scr3 - if(trim(cvars3d(m))=='q') q = scr3 - if(trim(cvars3d(m))=='oz') oz = scr3 - if(trim(cvars3d(m))=='cw') cwmr = scr3 - if(trim(cvars3d(m))=='ql') qlmr = scr3 - if(trim(cvars3d(m))=='qi') qimr = scr3 - if(trim(cvars3d(m))=='qr') qrmr = scr3 - if(trim(cvars3d(m))=='qs') qsmr = scr3 - if(trim(cvars3d(m))=='qg') qgmr = scr3 + if(trim(cvars3d(m))=='sf')then + u = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km) + else if(trim(cvars3d(m))=='vp') then + v = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km) + else if(trim(cvars3d(m))=='t') then + tv = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km) + else if(trim(cvars3d(m))=='q') then + q = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km) + else if(trim(cvars3d(m))=='oz') then + oz = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km) + else if(trim(cvars3d(m))=='cw') then + cwmr = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km) + else if(trim(cvars3d(m))=='ql') then + qlmr = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km) + else if(trim(cvars3d(m))=='qi') then + qimr = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km) + else if(trim(cvars3d(m))=='qr') then + qrmr = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km) + else if(trim(cvars3d(m))=='qs') then + qsmr = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km) + else if(trim(cvars3d(m))=='qg') then + qgmr = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km) + end if enddo -! convert ps from Pa to cb - ps=r0_001*ps ! convert t to virtual temperature tv=tv*(one+fv*q) -!--- now update pole values of atm_bundle using general_sub2grid (so halos also -! automatically updated. - - call create_grd23d_(grd2d,1) - call create_grd23d_(grd3d,grd%nsig) - - call update_scalar_poles_(grd2d,ps) - call update_vector_poles_(grd3d,u,v,clons,slons) - call update_scalar_poles_(grd3d,tv) - call update_scalar_poles_(grd3d,q) - call update_scalar_poles_(grd3d,oz) - if (icw>0) call update_scalar_poles_(grd3d,cwmr) - if (iql>0) call update_scalar_poles_(grd3d,qlmr) - if (iqi>0) call update_scalar_poles_(grd3d,qimr) - if (iqr>0) call update_scalar_poles_(grd3d,qrmr) - if (iqs>0) call update_scalar_poles_(grd3d,qsmr) - if (iqg>0) call update_scalar_poles_(grd3d,qgmr) - - call gsi_bundleputvar(atm_bundle,'ps',ps, ierr); iret = ierr - !call gsi_bundleputvar(atm_bundle,'sst',sst,ierr); iret = ierr + iret ! no sst for now - call gsi_bundleputvar(atm_bundle,'sf',u , ierr); iret = ierr + iret - call gsi_bundleputvar(atm_bundle,'vp',v , ierr); iret = ierr + iret - call gsi_bundleputvar(atm_bundle,'t' ,tv, ierr); iret = ierr + iret - call gsi_bundleputvar(atm_bundle,'q' ,q , ierr); iret = ierr + iret - call gsi_bundleputvar(atm_bundle,'oz',oz, ierr); iret = ierr + iret - if (icw>0) call gsi_bundleputvar(atm_bundle,'cw',cwmr,ierr); iret = ierr + iret - if (iql>0) call gsi_bundleputvar(atm_bundle,'ql',qlmr,ierr); iret = ierr + iret - if (iqi>0) call gsi_bundleputvar(atm_bundle,'qi',qimr,ierr); iret = ierr + iret - if (iqr>0) call gsi_bundleputvar(atm_bundle,'qr',qrmr,ierr); iret = ierr + iret - if (iqs>0) call gsi_bundleputvar(atm_bundle,'qs',qsmr,ierr); iret = ierr + iret - if (iqg>0) call gsi_bundleputvar(atm_bundle,'qg',qgmr,ierr); iret = ierr + iret - if ( iret /= 0 ) then - if ( mype == 0 ) then - write(6,'(A)') trim(myname_) // ': ERROR!' - write(6,'(A)') trim(myname_) // ': For now, GFS needs to put all MetFields: ps,u,v,(sf,vp)tv,q,oz,cw' - write(6,'(A)') trim(myname_) // ': but some have not been found. Aborting ... ' - write(6,'(A)') trim(myname_) // ': WARNING!' - write(6,'(3A,I5)') trim(myname_) // ': Trouble reading ensemble file : ', trim(filename), ', IRET = ', iret - endif - return - endif - - call general_sub2grid_destroy_info(grd2d,grd) - call general_sub2grid_destroy_info(grd3d,grd) - - if ( allocated(scr2) ) deallocate(scr2) - if ( allocated(scr3) ) deallocate(scr3) - return end subroutine move2bundle_ @@ -559,7 +498,7 @@ subroutine create_grd23d_(grd23d,nvert) end subroutine create_grd23d_ -subroutine update_scalar_poles_(grd,s) +subroutine update_halos_(grd,s) use kinds, only: i_kind,r_kind use general_sub2grid_mod, only: sub2grid_info,general_sub2grid,general_grid2sub @@ -598,9 +537,6 @@ subroutine update_scalar_poles_(grd,s) enddo call general_sub2grid(grd,sloc,work) - do k=kbegin_loc,kend_loc - call fillpoles_s_(work(1,:,:,k),nlon,nlat) - enddo call general_grid2sub(grd,work,sloc) ii=0 do k=1,nvert @@ -614,93 +550,9 @@ subroutine update_scalar_poles_(grd,s) deallocate(sloc,work) -end subroutine update_scalar_poles_ - -subroutine update_vector_poles_(grd,u,v,clons,slons) +end subroutine update_halos_ - use kinds, only: i_kind,r_kind - use constants, only: zero - use general_sub2grid_mod, only: sub2grid_info,general_sub2grid,general_grid2sub - - implicit none - - ! Declare local parameters - - ! Declare passed variables - type(sub2grid_info) ,intent(in ) :: grd - real(r_kind) ,intent(inout) :: u(grd%lat2,grd%lon2,grd%num_fields) - real(r_kind) ,intent(inout) :: v(grd%lat2,grd%lon2,grd%num_fields) - real(r_kind) ,intent(in ) :: clons(grd%nlon),slons(grd%nlon) - - ! Declare local variables - integer(i_kind) inner_vars,lat2,lon2,nlat,nlon,nvert,kbegin_loc,kend_loc,kend_alloc - integer(i_kind) ii,i,j,k - real(r_kind),allocatable,dimension(:) :: uloc,vloc - real(r_kind),allocatable,dimension(:,:,:,:) :: uwork,vwork - real(r_kind),allocatable,dimension(:,:) :: tempu,tempv - - lat2=grd%lat2 - lon2=grd%lon2 - nlat=grd%nlat - nlon=grd%nlon - nvert=grd%num_fields - inner_vars=grd%inner_vars - kbegin_loc=grd%kbegin_loc - kend_loc=grd%kend_loc - kend_alloc=grd%kend_alloc - allocate(uloc(lat2*lon2*nvert)) - allocate(vloc(lat2*lon2*nvert)) - allocate(uwork(inner_vars,nlat,nlon,kbegin_loc:kend_alloc)) - allocate(vwork(inner_vars,nlat,nlon,kbegin_loc:kend_alloc)) - allocate(tempu(nlat,nlon),tempv(nlat,nlon)) - uwork=zero ; vwork=zero ; uloc=zero ; vloc=zero - ii=0 - do k=1,nvert - do j=1,lon2 - do i=1,lat2 - ii=ii+1 - uloc(ii)=u(i,j,k) - vloc(ii)=v(i,j,k) - enddo - enddo - enddo - call general_sub2grid(grd,uloc,uwork) - call general_sub2grid(grd,vloc,vwork) - - do k=kbegin_loc,kend_loc - do j=1,nlon - do i=1,nlat - tempu(i,j)=uwork(1,i,j,k) - tempv(i,j)=vwork(1,i,j,k) - enddo - enddo - call fillpoles_v_(tempu,tempv,nlon,nlat,clons,slons) - do j=1,nlon - do i=1,nlat - uwork(1,i,j,k)=tempu(i,j) - vwork(1,i,j,k)=tempv(i,j) - enddo - enddo - enddo - call general_grid2sub(grd,uwork,uloc) - call general_grid2sub(grd,vwork,vloc) - ii=0 - do k=1,nvert - do j=1,lon2 - do i=1,lat2 - ii=ii+1 - u(i,j,k)=uloc(ii) - v(i,j,k)=vloc(ii) - enddo - enddo - enddo - - deallocate(uloc,uwork,tempu) - deallocate(vloc,vwork,tempv) - -end subroutine update_vector_poles_ - -subroutine ens_io_partition_(n_ens,io_pe,n_io_pe_s,n_io_pe_e,n_io_pe_em,i_ens) +subroutine ens_io_partition_(n_ens,ntindex,io_pe,n_io_pe_s,n_io_pe_e,n_io_pe_em,i_ens) ! do computation on all processors, then assign final local processor ! values. @@ -711,7 +563,7 @@ subroutine ens_io_partition_(n_ens,io_pe,n_io_pe_s,n_io_pe_e,n_io_pe_em,i_ens) implicit none ! Declare passed variables - integer(i_kind),intent(in ) :: n_ens + integer(i_kind),intent(in ) :: n_ens,ntindex integer(i_kind),intent( out) :: io_pe,n_io_pe_s,n_io_pe_e,n_io_pe_em,i_ens ! Declare local variables @@ -739,9 +591,11 @@ subroutine ens_io_partition_(n_ens,io_pe,n_io_pe_s,n_io_pe_e,n_io_pe_em,i_ens) endif ipe=ipe+jskip enddo - do n=1,n_ens - if(mype==0) write(6,'(2(a,1x,i5,1x))') 'reading ensemble member', n, 'on pe', io_pe0(n) - enddo + if(mype==0)then + do n=1,n_ens + write(6,'(3(a,1x,i5,1x))') 'reading ensemble member', n,' time level',ntindex,'on pe', io_pe0(n) + enddo + end if do n=1,n_ens if(mype==io_pe0(n)) then @@ -758,10 +612,10 @@ end subroutine ens_io_partition_ subroutine parallel_read_nemsio_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsig, & ias,jas,mas, & iasm,iaemz,jasm,jaemz,kasm,kaemz,masm,maemz, & - filename,init_head,clons,slons,filenamesfc) + filename,init_head,filenamesfc) use kinds, only: i_kind,r_kind,r_single - use constants, only: r60,r3600,zero,one,half,pi,deg2rad + use constants, only: r60,r3600,zero,one,half,deg2rad use nemsio_module, only: nemsio_init,nemsio_open,nemsio_close use ncepnems_io, only: error_msg,imp_physics use nemsio_module, only: nemsio_gfile,nemsio_getfilehead,nemsio_readrecv @@ -783,7 +637,6 @@ subroutine parallel_read_nemsio_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsi character(len=*), intent(in ) :: filename character(len=*), optional, intent(in) :: filenamesfc logical, intent(in ) :: init_head - real(r_kind), intent(inout) :: clons(nlon),slons(nlon) ! Declare local variables integer(i_kind) i,ii,j,jj,k,lonb,latb,levs,latb2,lonb2 @@ -809,8 +662,9 @@ subroutine parallel_read_nemsio_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsi real(r_kind) :: fhour type(nemsio_gfile) :: gfile type(nemsio_gfile) :: gfilesfc - real(r_kind),allocatable,dimension(:) :: rlats,rlons - real(r_single),allocatable,dimension(:) :: r4lats,r4lons + real(r_kind),allocatable,dimension(:) :: rlons + real(r_kind) :: clons(nlon),slons(nlon) + real(r_single),allocatable,dimension(:) :: r4lons if ( init_head)call nemsio_init(iret=iret) if (iret /= 0) call error_msg(trim(myname_),trim(filename),null,'init',istop,iret,.true.) @@ -849,24 +703,19 @@ subroutine parallel_read_nemsio_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsi endif endif -! obtain r4lats,r4lons,rlats,rlons,clons,slons exactly as computed in general_read_gfsatm_nems: +! obtain r4lons,rlons,clons,slons exactly as computed in general_read_gfsatm_nems: - allocate(rlats(latb+2),rlons(lonb),r4lats(lonb*latb),r4lons(lonb*latb)) - call nemsio_getfilehead(gfile,lat=r4lats,iret=iret) + allocate(rlons(lonb),r4lons(lonb*latb)) call nemsio_getfilehead(gfile,lon=r4lons,iret=iret) - do j=1,latb - rlats(latb+2-j)=deg2rad*r4lats(lonb/2+(j-1)*lonb) - enddo do j=1,lonb rlons(j)=deg2rad*r4lons(j) enddo - deallocate(r4lats,r4lons) - rlats(1)=-half*pi - rlats(latb+2)=half*pi + deallocate(r4lons) do j=1,lonb clons(j)=cos(rlons(j)) slons(j)=sin(rlons(j)) enddo + deallocate(rlons) fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600 odate(1) = idate(4) !hour @@ -913,38 +762,48 @@ subroutine parallel_read_nemsio_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsi end if end if call move1_(work,temp3(:,:,k,k3),nlon,nlat) + call fillpoles_ss_(temp3(:,:,k,k3),nlon,nlat) elseif(trim(cvars3d(k3))=='ql') then call nemsio_readrecv(gfile,'clwmr','mid layer',k,work,iret=iret) if (iret /= 0) call error_msg(trim(myname_),trim(filename),'clwmr','read',istop+8,iret) call move1_(work,temp3(:,:,k,k3),nlon,nlat) + call fillpoles_ss_(temp3(:,:,k,k3),nlon,nlat) elseif(trim(cvars3d(k3))=='qi') then call nemsio_readrecv(gfile,'icmr','mid layer',k,work,iret=iret) if (iret /= 0) call error_msg(trim(myname_),trim(filename),'icmr','read',istop+9,iret) call move1_(work,temp3(:,:,k,k3),nlon,nlat) + call fillpoles_ss_(temp3(:,:,k,k3),nlon,nlat) elseif(trim(cvars3d(k3))=='qr') then call nemsio_readrecv(gfile,'rwmr','mid layer',k,work,iret=iret) if (iret /= 0) call error_msg(trim(myname_),trim(filename),'rwmr','read',istop+10,iret) call move1_(work,temp3(:,:,k,k3),nlon,nlat) + call fillpoles_ss_(temp3(:,:,k,k3),nlon,nlat) elseif(trim(cvars3d(k3))=='qs') then call nemsio_readrecv(gfile,'snmr','mid layer',k,work,iret=iret) + call fillpoles_ss_(temp3(:,:,k,k3),nlon,nlat) if (iret /= 0) call error_msg(trim(myname_),trim(filename),'snmr','read',istop+11,iret) call move1_(work,temp3(:,:,k,k3),nlon,nlat) + call fillpoles_ss_(temp3(:,:,k,k3),nlon,nlat) elseif(trim(cvars3d(k3))=='qg') then call nemsio_readrecv(gfile,'grle','mid layer',k,work,iret=iret) if (iret /= 0) call error_msg(trim(myname_),trim(filename),'grle','read',istop+12,iret) call move1_(work,temp3(:,:,k,k3),nlon,nlat) + call fillpoles_ss_(temp3(:,:,k,k3),nlon,nlat) elseif(trim(cvars3d(k3))=='oz') then call nemsio_readrecv(gfile,'o3mr','mid layer',k,work,iret=iret) if (iret /= 0) call error_msg(trim(myname_),trim(filename),'o3mr','read',istop+5,iret,.true.) call move1_(work,temp3(:,:,k,k3),nlon,nlat) + call fillpoles_ss_(temp3(:,:,k,k3),nlon,nlat) elseif(trim(cvars3d(k3))=='q') then call nemsio_readrecv(gfile,'spfh','mid layer',k,work,iret=iret) if (iret /= 0) call error_msg(trim(myname_),trim(filename),trim(cvars3d(k3)),'read',istop+4,iret,.true.) call move1_(work,temp3(:,:,k,k3),nlon,nlat) + call fillpoles_ss_(temp3(:,:,k,k3),nlon,nlat) elseif(trim(cvars3d(k3))=='t') then call nemsio_readrecv(gfile,'tmp','mid layer',k,work,iret=iret) if (iret /= 0) call error_msg(trim(myname_),trim(filename),'tmp','read',istop+3,iret,.true.) call move1_(work,temp3(:,:,k,k3),nlon,nlat) + call fillpoles_ss_(temp3(:,:,k,k3),nlon,nlat) elseif(trim(cvars3d(k3))=='sf') then call nemsio_readrecv(gfile,'ugrd','mid layer',k,work,iret=iret) if (iret /= 0) call error_msg(trim(myname_),trim(filename),'ugrd','read',istop+1,iret,.true.) @@ -956,6 +815,9 @@ subroutine parallel_read_nemsio_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsi endif enddo enddo + do k=1,nsig + call fillpoles_sv_(temp3(:,:,k,k3u),temp3(:,:,k,k3v),nlon,nlat,clons,slons) + end do ! if (k3u==0.or.k3v==0.or.k3t==0.or.k3q==0.or.k3cw==0.or.k3oz==0) & if (k3u==0.or.k3v==0.or.k3t==0.or.k3q==0.or.k3oz==0) & write(6,'(" WARNING, problem with one of k3-")') @@ -975,6 +837,7 @@ subroutine parallel_read_nemsio_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsi if (iret /= 0) call error_msg(trim(myname_),trim(filename),'hgt','read',istop+8,iret,.true.) !work=r0_001*work ! convert Pa to cb ! postpone this calculation call move1_(work,temp2(:,:,k2),nlon,nlat) + call fillpoles_ss_(temp2(:,:,k2),nlon,nlat) endif enddo deallocate(work) @@ -997,6 +860,7 @@ subroutine parallel_read_nemsio_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsi enddo enddo enddo + deallocate(temp3) do k2=1,nc2d m_cvars2d(k2)=kf+1 kf=kf+1 @@ -1011,7 +875,6 @@ subroutine parallel_read_nemsio_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsi enddo enddo - deallocate(temp3) deallocate(temp2) end subroutine parallel_read_nemsio_state_ @@ -1019,7 +882,7 @@ end subroutine parallel_read_nemsio_state_ subroutine parallel_read_gfsnc_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsig, & ias,jas,mas, & iasm,iaemz,jasm,jaemz,kasm,kaemz,masm,maemz, & - filename,init_head,clons,slons) + filename) !$$$ subprogram documentation block ! . . . . ! subprogram: parallel_read_gfsnc_state_ read GFS netCDF ensemble member @@ -1031,7 +894,7 @@ subroutine parallel_read_gfsnc_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsig !$$$ use kinds, only: i_kind,r_kind,r_single - use constants, only: r60,r3600,zero,one,half,pi,deg2rad + use constants, only: r60,r3600,zero,one,half,deg2rad use control_vectors, only: cvars2d,cvars3d,nc2d,nc3d use general_sub2grid_mod, only: sub2grid_info use module_fv3gfs_ncio, only: Dataset, Variable, Dimension, open_dataset,& @@ -1048,18 +911,17 @@ subroutine parallel_read_gfsnc_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsig integer(i_kind), intent(inout) :: m_cvars2d(nc2d),m_cvars3d(nc3d) real(r_single), intent(inout) :: en_full(iasm:iaemz,jasm:jaemz,kasm:kaemz,masm:maemz) character(len=*), intent(in ) :: filename - logical, intent(in ) :: init_head - real(r_kind), intent(inout) :: clons(nlon),slons(nlon) ! Declare local variables integer(i_kind) i,ii,j,jj,k,lonb,latb,levs,kr integer(i_kind) k2,k3,k3u,k3v,k3t,k3q,k3cw,k3oz,kf character(len=120) :: myname_ = 'parallel_read_gfsnc_state_' - real(r_single),allocatable,dimension(:,:,:) :: temp2, rwork3d1, rwork3d2 - real(r_single),allocatable,dimension(:,:) :: rwork2d + real(r_single),allocatable,dimension(:,:,:) :: rwork3d1, rwork3d2 + real(r_single),allocatable,dimension(:,:) :: temp2,rwork2d real(r_single),allocatable,dimension(:,:,:,:) :: temp3 - real(r_kind),allocatable,dimension(:) :: rlats,rlons - real(r_kind),allocatable,dimension(:) :: rlats_tmp,rlons_tmp + real(r_kind),allocatable,dimension(:) :: rlons_tmp + real(r_kind) :: clons(nlon),slons(nlon) + real(r_kind) :: rlons type(Dataset) :: atmges type(Dimension) :: ncdim @@ -1079,70 +941,66 @@ subroutine parallel_read_gfsnc_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsig call die(myname_, ': ***ERROR*** incorrect resolution',101) endif -! obtain rlats_tmp,rlons_tnp,rlats,rlons,clons,slons exactly as computed in general_read_gfsatm_nems: +! obtain rlons_tnp,rlons,clons,slons exactly as computed in general_read_gfsatm_nems: - allocate(rlats(latb+2),rlons(lonb)) call read_vardata(atmges, 'grid_xt', rlons_tmp) - call read_vardata(atmges, 'grid_yt', rlats_tmp) - do j=1,latb - rlats(latb+2-j)=deg2rad*rlats_tmp(j) - enddo do j=1,lonb - rlons(j)=deg2rad*rlons_tmp(j) - enddo - deallocate(rlats_tmp,rlons_tmp) - rlats(1)=-half*pi - rlats(latb+2)=half*pi - do j=1,lonb - clons(j)=cos(rlons(j)) - slons(j)=sin(rlons(j)) + rlons=deg2rad*rlons_tmp(j) + clons(j)=cos(rlons) + slons(j)=sin(rlons) enddo + deallocate(rlons_tmp) - allocate(rwork3d2(nlon,(nlat-2),nsig)) + allocate(rwork3d1(nlon,(nlat-2),nsig)) allocate(temp3(nlat,nlon,nsig,nc3d)) - allocate(temp2(nlat,nlon,nc2d)) k3u=0 ; k3v=0 ; k3t=0 ; k3q=0 ; k3cw=0 ; k3oz=0 do k3=1,nc3d - if(cvars3d(k3)=='sf') k3u=k3 - if(cvars3d(k3)=='vp') k3v=k3 - if(cvars3d(k3)=='t') k3t=k3 - if(cvars3d(k3)=='q') k3q=k3 - if(cvars3d(k3)=='cw') k3cw=k3 - if(cvars3d(k3)=='oz') k3oz=k3 if (trim(cvars3d(k3))=='cw') then + k3cw=k3 call read_vardata(atmges, 'clwmr', rwork3d1) - rwork3d2 = 0 + allocate(rwork3d2(nlon,(nlat-2),nsig)) + rwork3d2 = 0._r_single call read_vardata(atmges, 'icmr', rwork3d2) rwork3d1 = rwork3d1 + rwork3d2 + deallocate(rwork3d2) do k=1,nsig kr = levs+1-k call move1_(rwork3d1(:,:,kr),temp3(:,:,k,k3),nlon,nlat) + call fillpoles_ss_(temp3(:,:,k,k3),nlon,nlat) end do else if(trim(cvars3d(k3))=='oz') then + k3oz=k3 call read_vardata(atmges, 'o3mr', rwork3d1) do k=1,nsig kr = levs+1-k call move1_(rwork3d1(:,:,kr),temp3(:,:,k,k3),nlon,nlat) + call fillpoles_ss_(temp3(:,:,k,k3),nlon,nlat) end do else if(trim(cvars3d(k3))=='q') then + k3q=k3 call read_vardata(atmges, 'spfh', rwork3d1) do k=1,nsig kr = levs+1-k call move1_(rwork3d1(:,:,kr),temp3(:,:,k,k3),nlon,nlat) + call fillpoles_ss_(temp3(:,:,k,k3),nlon,nlat) end do else if(trim(cvars3d(k3))=='t') then + k3t=k3 call read_vardata(atmges, 'tmp', rwork3d1) do k=1,nsig kr = levs+1-k call move1_(rwork3d1(:,:,kr),temp3(:,:,k,k3),nlon,nlat) + call fillpoles_ss_(temp3(:,:,k,k3),nlon,nlat) end do else if(trim(cvars3d(k3))=='sf') then + k3u=k3 call read_vardata(atmges, 'ugrd', rwork3d1) do k=1,nsig kr = levs+1-k call move1_(rwork3d1(:,:,kr),temp3(:,:,k,k3),nlon,nlat) end do else if(trim(cvars3d(k3))=='vp') then + k3v=k3 call read_vardata(atmges, 'vgrd', rwork3d1) do k=1,nsig kr = levs+1-k @@ -1150,18 +1008,14 @@ subroutine parallel_read_gfsnc_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsig end do end if enddo + deallocate(rwork3d1) + if (k3u==0.or.k3v==0.or.k3t==0.or.k3q==0.or.k3cw==0.or.k3oz==0) & write(6,'(" WARNING, problem with one of k3-")') - temp2=zero - do k2=1,nc2d - if(trim(cvars2d(k2))=='ps') then - call read_vardata(atmges, 'pressfc', rwork2d) - call move1_(rwork2d,temp2(:,:,k2),nlon,nlat) - endif - enddo - deallocate(rwork2d, rwork3d1) - deallocate(rwork3d2) + do k=1,nsig + call fillpoles_sv_(temp3(:,:,k,k3u),temp3(:,:,k,k3v),nlon,nlat,clons,slons) + end do ! move temp2,temp3 to en_full kf=0 @@ -1180,29 +1034,41 @@ subroutine parallel_read_gfsnc_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsig enddo enddo enddo + + deallocate(temp3) + allocate(temp2(nlat,nlon)) + allocate(rwork2d(nlon,(nlat-2))) do k2=1,nc2d - m_cvars2d(k2)=kf+1 + if(trim(cvars2d(k2))=='ps') then + call read_vardata(atmges, 'pressfc', rwork2d) + call move1_(rwork2d,temp2,nlon,nlat) + call fillpoles_ss_(temp2,nlon,nlat) + else + temp2=zero + endif + kf=kf+1 + m_cvars2d(k2)=kf jj=jas-1 do j=1,nlon jj=jj+1 ii=ias-1 do i=1,nlat ii=ii+1 - en_full(ii,jj,kf,mas)=temp2(i,j,k2) + en_full(ii,jj,kf,mas)=temp2(i,j) enddo enddo enddo - deallocate(temp3) + deallocate(rwork2d) deallocate(temp2) end subroutine parallel_read_gfsnc_state_ -subroutine fillpoles_s_(temp,nlon,nlat) +subroutine fillpoles_ss_(temp,nlon,nlat) !$$$ subprogram documentation block ! . . . . -! subprogram: fillpoles_s_ make pole points average of nearest pole row +! subprogram: fillpoles_ss_ make pole points average of nearest pole row ! prgmmr: parrish org: emc/ncep date: 2016-10-14 ! ! abstract: make pole points average of nearest pole row. @@ -1225,13 +1091,13 @@ subroutine fillpoles_s_(temp,nlon,nlat) ! !$$$ - use kinds, only: i_kind,r_kind + use kinds, only: i_kind,r_kind,r_single use constants, only: zero,one implicit none integer(i_kind),intent(in ) :: nlon,nlat - real(r_kind), intent(inout) :: temp(nlat,nlon) + real(r_single), intent(inout) :: temp(nlat,nlon) integer(i_kind) nlatm1,i real(r_kind) sumn,sums,rnlon @@ -1254,12 +1120,12 @@ subroutine fillpoles_s_(temp,nlon,nlat) temp(nlat,i)=sumn end do -end subroutine fillpoles_s_ +end subroutine fillpoles_ss_ -subroutine fillpoles_v_(tempu,tempv,nlon,nlat,clons,slons) +subroutine fillpoles_sv_(tempu,tempv,nlon,nlat,clons,slons) !$$$ subprogram documentation block ! . . . . -! subprogram: fillpoles_v_ create vector values at pole from nearest pole row +! subprogram: fillpoles_sv_ create vector values at pole from nearest pole row ! prgmmr: parrish org: emc/ncep date: 2016-10-14 ! ! abstract: create vector values at pole from nearest pole row. @@ -1283,13 +1149,13 @@ subroutine fillpoles_v_(tempu,tempv,nlon,nlat,clons,slons) ! !$$$ - use kinds, only: i_kind,r_kind + use kinds, only: i_kind,r_kind,r_single use constants, only: zero implicit none integer(i_kind),intent(in ) :: nlon,nlat - real(r_kind), intent(inout) :: tempu(nlat,nlon),tempv(nlat,nlon) + real(r_single), intent(inout) :: tempu(nlat,nlon),tempv(nlat,nlon) real(r_kind), intent(in ) :: clons(nlon),slons(nlon) integer(i_kind) i @@ -1317,7 +1183,7 @@ subroutine fillpoles_v_(tempu,tempv,nlon,nlat,clons,slons) tempv(1,i )= polsu*slons(i)-polsv*clons(i) end do -end subroutine fillpoles_v_ +end subroutine fillpoles_sv_ subroutine move1_(work,temp,nlon,nlat) !$$$ subprogram documentation block @@ -1356,8 +1222,7 @@ subroutine move1_(work,temp,nlon,nlat) integer(i_kind) ii,i,j ii=0 - temp(1,:)=zero - temp(nlat,:)=zero +! Polar points will be filled in later do i=nlat-1,2,-1 do j=1,nlon ii=ii+1 diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90 index 45a2c83bea..d03715059f 100644 --- a/src/gsi/crtm_interface.f90 +++ b/src/gsi/crtm_interface.f90 @@ -854,7 +854,7 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo if (n_actual_aerosols_wk>0 .or. n_clouds_fwd_wk>0 .and. imp_physics==11) then if (mype==0) write(6,*)myname_,':initial and load GFDL saturation water vapor pressure tables' - + allocate(table (length)) allocate(table2(length)) allocate(tablew(length)) diff --git a/src/gsi/general_spectral_transforms.f90 b/src/gsi/general_spectral_transforms.f90 index 86cf88e42c..d4f0959489 100644 --- a/src/gsi/general_spectral_transforms.f90 +++ b/src/gsi/general_spectral_transforms.f90 @@ -399,7 +399,6 @@ subroutine sfilter(grd,sp,filter,grid) end do return - return end subroutine sfilter diff --git a/src/gsi/genex_mod.f90 b/src/gsi/genex_mod.f90 index a9e7ea2142..278dbdf2f5 100644 --- a/src/gsi/genex_mod.f90 +++ b/src/gsi/genex_mod.f90 @@ -775,27 +775,6 @@ subroutine genex_destroy_info(s) if(associated(s%kabe_rc)) deallocate(s%kabe_rc) if(associated(s%mabs_rc)) deallocate(s%mabs_rc) if(associated(s%mabe_rc)) deallocate(s%mabe_rc) -! following is probably redundant - s%lefts => NULL() - s%rights => NULL() - s%numl => NULL() - s%numrc => NULL() - s%iabs_l => NULL() - s%iabe_l => NULL() - s%jabs_l => NULL() - s%jabe_l => NULL() - s%kabs_l => NULL() - s%kabe_l => NULL() - s%mabs_l => NULL() - s%mabe_l => NULL() - s%iabs_rc => NULL() - s%iabe_rc => NULL() - s%jabs_rc => NULL() - s%jabe_rc => NULL() - s%kabs_rc => NULL() - s%kabe_rc => NULL() - s%mabs_rc => NULL() - s%mabe_rc => NULL() s%lallocated=.false. end subroutine genex_destroy_info diff --git a/src/gsi/get_gefs_ensperts_dualres.f90 b/src/gsi/get_gefs_ensperts_dualres.f90 index 39ec7541da..5e657349eb 100644 --- a/src/gsi/get_gefs_ensperts_dualres.f90 +++ b/src/gsi/get_gefs_ensperts_dualres.f90 @@ -51,7 +51,6 @@ subroutine get_gefs_ensperts_dualres use mpeu_util, only: die use gridmod, only: idsl5 use hybrid_ensemble_parameters, only: n_ens,write_ens_sprd,oz_univ_static,ntlevs_ens - use hybrid_ensemble_parameters, only: sst_staticB use hybrid_ensemble_parameters, only: en_perts,ps_bar,nelen use constants,only: zero,zero_single,half,fv,rd_over_cp,one,qcmin use mpimod, only: mpi_comm_world,mype,npe @@ -76,20 +75,17 @@ subroutine get_gefs_ensperts_dualres real(r_kind),pointer,dimension(:,:,:) :: q ! real(r_kind),dimension(grd_ens%nlat,grd_ens%nlon):: sst_full,dum real(r_kind),pointer,dimension(:,:,:):: p3 - real(r_kind),pointer,dimension(:,:):: p2 - real(r_single),pointer,dimension(:,:,:):: w3 - real(r_single),pointer,dimension(:,:):: w2 - real(r_kind),pointer,dimension(:,:,:):: x3 real(r_kind),pointer,dimension(:,:):: x2 type(gsi_bundle),allocatable,dimension(:) :: en_read - type(gsi_bundle),allocatable,dimension(:) :: en_bar + type(gsi_bundle):: en_bar ! type(gsi_grid) :: grid_ens - real(r_kind) bar_norm,sig_norm,kapr,kap1,rh - real(r_kind),allocatable,dimension(:,:):: z,sst2 + real(r_kind) bar_norm,sig_norm,kapr,kap1 +! real(r_kind),allocatable,dimension(:,:):: z,sst2 real(r_kind),allocatable,dimension(:,:,:) :: tsen,prsl,pri,qs ! integer(i_kind),dimension(grd_ens%nlat,grd_ens%nlon):: idum - integer(i_kind) istatus,iret,i,ic2,ic3,j,k,n,mm1,iderivative,im,jm,km,m,ipic + integer(i_kind) istatus,iret,i,ic3,j,k,n,iderivative,im,jm,km,m,ipic +! integer(i_kind) mm1 integer(i_kind) ipc3d(nc3d),ipc2d(nc2d) integer(i_kind) ier ! integer(i_kind) il,jl @@ -119,24 +115,20 @@ subroutine get_gefs_ensperts_dualres call stop2(999) endif - mm1=mype+1 - kap1=rd_over_cp+one - kapr=one/rd_over_cp im=en_perts(1,1)%grid%im jm=en_perts(1,1)%grid%jm km=en_perts(1,1)%grid%km + bar_norm = one/float(n_ens) + sig_norm=sqrt(one/max(one,n_ens-one)) ! Create temporary communication information for read ensemble routines call gsi_enscoupler_create_sub2grid_info(grd_tmp,km,npe,grd_ens) ! Allocate bundle to hold mean of ensemble members - allocate(en_bar(ntlevs_ens)) - do m=1,ntlevs_ens - call gsi_bundlecreate(en_bar(m),en_perts(1,1)%grid,'ensemble',istatus,names2d=cvars2d,names3d=cvars3d) - if ( istatus /= 0 ) & - call die('get_gefs_ensperts_dualres',': trouble creating en_bar bundle, istatus =',istatus) - end do + call gsi_bundlecreate(en_bar,en_perts(1,1)%grid,'ensemble',istatus,names2d=cvars2d,names3d=cvars3d) + if ( istatus /= 0 ) & + call die('get_gefs_ensperts_dualres',': trouble creating en_bar bundle, istatus =',istatus) ! Allocate bundle used for reading members allocate(en_read(n_ens)) @@ -146,15 +138,21 @@ subroutine get_gefs_ensperts_dualres call die('get_gefs_ensperts_dualres',': trouble creating en_read bundle, istatus =',istatus) end do - allocate(z(im,jm)) - allocate(sst2(im,jm)) +! allocate(z(im,jm)) +! allocate(sst2(im,jm)) - sst2=zero ! for now, sst not used in ensemble perturbations, so if sst array is called for +! sst2=zero ! for now, sst not used in ensemble perturbations, so if sst array is called for ! then sst part of en_perts will be zero when sst2=zero - ntlevs_ens_loop: do m=1,ntlevs_ens +!$omp parallel do schedule(dynamic,1) private(m,n) + do m=1,ntlevs_ens + do n=1,n_ens + en_perts(n,m)%valuesr4=zero_single + end do + end do - en_bar(m)%values=zero + + ntlevs_ens_loop: do m=1,ntlevs_ens call gsi_enscoupler_get_user_Nens(grd_tmp,n_ens,m,en_read,iret) @@ -168,11 +166,11 @@ subroutine get_gefs_ensperts_dualres cycle endif - n_ens_loop: do n=1,n_ens - - en_perts(n,m)%valuesr4=zero + if (.not.q_hyb_ens) then !use RH + kap1=rd_over_cp+one + kapr=one/rd_over_cp + do n=1,n_ens - if (.not.q_hyb_ens) then !use RH call gsi_bundlegetpointer(en_read(n),'ps',ps,ier);istatus=ier call gsi_bundlegetpointer(en_read(n),'t' ,tv,ier);istatus=istatus+ier call gsi_bundlegetpointer(en_read(n),'q' ,q ,ier);istatus=istatus+ier @@ -209,10 +207,24 @@ subroutine get_gefs_ensperts_dualres ice=.true. iderivative=0 call genqsat(qs,tsen,prsl,im,jm,km,ice,iderivative) - deallocate(tsen,prsl) - end if + do k=1,km + do j=1,jm + do i=1,im + q(i,j,k)=q(i,j,k)/qs(i,j,k) + end do + end do + end do + deallocate(tsen,prsl,qs) + enddo + end if + + + en_bar%values=zero + + n_ens_loop: do n=1,n_ens -!_$omp parallel do schedule(dynamic,1) private(i,k,j,ic3,rh) + +!$omp parallel do schedule(dynamic,1) private(i,k,j,ic3,hydrometeor,istatus,p3) do ic3=1,nc3d hydrometeor = trim(cvars3d(ic3))=='cw' .or. trim(cvars3d(ic3))=='ql' .or. & @@ -222,145 +234,56 @@ subroutine get_gefs_ensperts_dualres call gsi_bundlegetpointer(en_read(n),trim(cvars3d(ic3)),p3,istatus) if(istatus/=0) then - write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' from read in member ',m - call stop2(999) - end if - call gsi_bundlegetpointer(en_perts(n,m),trim(cvars3d(ic3)),w3,istatus) - if(istatus/=0) then - write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for ensemble member ',n - call stop2(999) - end if - call gsi_bundlegetpointer(en_bar(m),trim(cvars3d(ic3)),x3,istatus) - if(istatus/=0) then - write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for en_bar' + write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' from read in member ',n,m call stop2(999) end if - if ( trim(cvars3d(ic3)) == 'q' ) then - if (.not.q_hyb_ens) then !use RH - do k=1,km - do j=1,jm - do i=1,im - rh=p3(i,j,k)/qs(i,j,k) - w3(i,j,k) = rh - x3(i,j,k)=x3(i,j,k)+rh - end do - end do - end do - cycle - end if - end if if ( hydrometeor ) then -!$omp parallel do schedule(dynamic,1) private(i,j,k) do k=1,km do j=1,jm do i=1,im - w3(i,j,k) = max(p3(i,j,k),qcmin) - x3(i,j,k)=x3(i,j,k)+max(p3(i,j,k),qcmin) + p3(i,j,k) = max(p3(i,j,k),qcmin) end do end do end do - cycle - end if - if ( trim(cvars3d(ic3)) == 'oz' .and. oz_univ_static ) then - w3 = zero_single - cycle + else if ( trim(cvars3d(ic3)) == 'oz' .and. oz_univ_static ) then + p3 = zero end if -!$omp parallel do schedule(dynamic,1) private(i,j,k) - do k=1,km - do j=1,jm - do i=1,im - w3(i,j,k) = p3(i,j,k) - x3(i,j,k)=x3(i,j,k)+p3(i,j,k) - end do - end do - end do - end do !c3d - if (.not.q_hyb_ens) deallocate(qs) - -!_$omp parallel do schedule(dynamic,1) private(i,j,ic2,ipic) - do ic2=1,nc2d + do i=1,nelen + en_perts(n,m)%valuesr4(i)=en_read(n)%values(i) + en_bar%values(i)=en_bar%values(i)+en_read(n)%values(i) + end do - call gsi_bundlegetpointer(en_read(n),trim(cvars2d(ic2)),p2,istatus) - if(istatus/=0) then - write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' from read in member ',m - call stop2(999) - end if - call gsi_bundlegetpointer(en_perts(n,m),trim(cvars2d(ic2)),w2,istatus) - if(istatus/=0) then - write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for ens-pert member ',n - call stop2(999) - end if - call gsi_bundlegetpointer(en_bar(m),trim(cvars2d(ic2)),x2,istatus) - if(istatus/=0) then - write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for en_bar' - call stop2(999) - end if -!$omp parallel do schedule(dynamic,1) private(i,j) - do j=1,jm - do i=1,im - w2(i,j)=p2(i,j) - x2(i,j)=x2(i,j)+p2(i,j) - end do - end do - - if (sst_staticB.and.trim(cvars2d(ic2))=='sst') then - w2 = zero - x2 = zero ! NOTE: if anyone implements alternative use of SST (as from sst2) care need ! be given to those applications getting SST directly from the members of ! the ensemble for which this code is already handling - i.e., I don't ! know who would want to commented out code below but be mindful ! of how it interacts with option sst_staticB, please - Todling. -!_$omp parallel do schedule(dynamic,1) private(i,j) -! do j=1,jm -! do i=1,im -! w2(i,j) = sst2(i,j) -! x2(i,j)=x2(i,j)+sst2(i,j) -! end do -! end do - cycle - end if - - end do - end do n_ens_loop ! end do over ensemble - end do ntlevs_ens_loop !end do over bins - - do n=n_ens,1,-1 - call gsi_bundledestroy(en_read(n),istatus) - if ( istatus /= 0 ) & - call die('get_gefs_ensperts_dualres',': trouble destroying en_read bundle, istatus = ', istatus) - end do - deallocate(en_read) - call gsi_enscoupler_destroy_sub2grid_info(grd_tmp) + end do n_ens_loop ! end do over ensemble -! Copy pbar to module array. ps_bar may be needed for vertical localization -! in terms of scale heights/normalized p/p -! Convert to mean - bar_norm = one/float(n_ens) - sig_norm=sqrt(one/max(one,n_ens-one)) -!$omp parallel do schedule(dynamic,1) private(i,j,k,n,m,ic2,ic3,ipic,x2) - do m=1,ntlevs_ens do i=1,nelen - en_bar(m)%values(i)=en_bar(m)%values(i)*bar_norm + en_bar%values(i)=en_bar%values(i)*bar_norm end do ! Before converting to perturbations, get ensemble spread - !-- if (m == 1 .and. write_ens_sprd ) call ens_spread_dualres(en_bar(1),1) + !-- if (m == 1 .and. write_ens_sprd ) call ens_spread_dualres(en_bar,1) !!! it is not clear of the next statement is thread/$omp safe. - if (write_ens_sprd ) call ens_spread_dualres(en_bar(m),m) + if (write_ens_sprd ) call ens_spread_dualres(en_bar,m) - call gsi_bundlegetpointer(en_bar(m),'ps',x2,istatus) + call gsi_bundlegetpointer(en_bar,'ps',x2,istatus) if(istatus/=0) & call die('get_gefs_ensperts_dualres:',' error retrieving pointer to (ps) for en_bar, istatus = ', istatus) +! Copy pbar to module array. ps_bar may be needed for vertical localization +! in terms of scale heights/normalized p/p +! Convert to mean do j=1,jm do i=1,im ps_bar(i,j,m)=x2(i,j) @@ -369,14 +292,15 @@ subroutine get_gefs_ensperts_dualres ! Convert ensemble members to perturbations +!$omp parallel do schedule(dynamic,1) private(n,i,ic3,ipic,k,j) do n=1,n_ens do i=1,nelen - en_perts(n,m)%valuesr4(i)=en_perts(n,m)%valuesr4(i)-en_bar(m)%values(i) + en_perts(n,m)%valuesr4(i)=en_perts(n,m)%valuesr4(i)-en_bar%values(i) end do if(.not. q_hyb_ens) then do ic3=1,nc3d - ipic=ipc3d(ic3) if(trim(cvars3d(ic3)) == 'q' .or. trim(cvars3d(ic3)) == 'Q')then + ipic=ipc3d(ic3) do k=1,km do j=1,jm do i=1,im @@ -392,8 +316,17 @@ subroutine get_gefs_ensperts_dualres en_perts(n,m)%valuesr4(i)=en_perts(n,m)%valuesr4(i)*sig_norm end do end do + end do ntlevs_ens_loop !end do over bins + + do n=n_ens,1,-1 + call gsi_bundledestroy(en_read(n),istatus) + if ( istatus /= 0 ) & + call die('get_gefs_ensperts_dualres',': trouble destroying en_read bundle, istatus = ', istatus) end do + deallocate(en_read) + call gsi_enscoupler_destroy_sub2grid_info(grd_tmp) +! mm1=mype+1 ! since initial version is ignoring sst perturbations, skip following code for now. revisit ! later--creating general_read_gfssfc, analogous to general_read_gfsatm above. !! GET SST PERTURBATIONS HERE @@ -436,17 +369,8 @@ subroutine get_gefs_ensperts_dualres ! end do ! end do - do m=ntlevs_ens,1,-1 - call gsi_bundledestroy(en_bar(m),istatus) - if(istatus/=0) then - write(6,*)' in get_gefs_ensperts_dualres: trouble destroying en_bar bundle' - call stop2(999) - endif - end do - - deallocate(sst2) - deallocate(z) - deallocate(en_bar) +! deallocate(sst2) +! deallocate(z) return end subroutine get_gefs_ensperts_dualres diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 997f76bb73..e1cc923f18 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -95,7 +95,7 @@ module gsimod factql,factqi,factqr,factqs,factqg, & 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,lgschmidt,diag_precon,step_start,pseudo_q2,& + bcoption,diurnalbc,print_diag_pcg,tsensible,diag_precon,step_start,pseudo_q2,& clip_supersaturation,cnvw_option use state_vectors, only: init_anasv,final_anasv use control_vectors, only: init_anacv,final_anacv,nrf,nvars,nrf_3d,cvars3d,cvars2d,& @@ -546,8 +546,6 @@ module gsimod ! preserve_restart_date - if true, then do not update regional restart file date. ! tsensible - option to use sensible temperature as the analysis variable. works ! only for twodvar_regional=.true. -! lgschmidt - option for re-biorthogonalization of the {gradx} and {grady} sets -! from pcgsoi when twodvar_regional=.true. ! hilbert_curve - option for hilbert-curve based cross-validation. works only ! with twodvar_regional=.true. ! neutral_stability_windfact_2dvar - option to use simple, similarity @@ -679,7 +677,7 @@ module gsimod nwrvecs,iorthomax,ladtest,ladtest_obs, lgrtest,lobskeep,lsensrecompute,jsiga,ltcost, & lobsensfc,lobsensjb,lobsensincr,lobsensadj,lobsensmin,iobsconv, & idmodel,iwrtinc,lwrite4danl,nhr_anal,jiterstart,jiterend,lobserver,lanczosave,llancdone, & - lferrscale,print_diag_pcg,tsensible,lgschmidt,lread_obs_save,lread_obs_skip, & + lferrscale,print_diag_pcg,tsensible,lread_obs_save,lread_obs_skip, & use_gfs_ozone,check_gfs_ozone_date,regional_ozone,lwrite_predterms,& lwrite_peakwt,use_gfs_nemsio,use_gfs_ncio,sfcnst_comb,liauon,use_prepb_satwnd,l4densvar,ens_nstarthr,& use_gfs_stratosphere,pblend0,pblend1,step_start,diag_precon,lrun_subdirs,& @@ -1454,12 +1452,6 @@ subroutine gsimain_initialize call die(myname_,'Options l4dvar and reduce_diag not allowed together',99) end if -! Diagonal preconditioning is necessary for new bias correction - if(newpc4pred .and. .not. diag_precon)then - diag_precon=.true. - step_start=8.e-4_r_kind - end if - if( (.not.l4dvar) .and. (.not.l4densvar) ) ljcdfi=.false. if (l4dvar.and.lsensrecompute) then diff --git a/src/gsi/hybrid_ensemble_isotropic.F90 b/src/gsi/hybrid_ensemble_isotropic.F90 index 9231392f6b..e0a6e3fd1f 100644 --- a/src/gsi/hybrid_ensemble_isotropic.F90 +++ b/src/gsi/hybrid_ensemble_isotropic.F90 @@ -2593,8 +2593,8 @@ subroutine sqrt_beta_s_mult_cvec(grady) call stop2(999) endif -!$omp parallel do schedule(dynamic,1) private(ic3,ic2,k,j,i,ii) ! multiply by sqrt_beta_s +!$omp parallel do schedule(dynamic,1) private(ic3,ic2,k,j,i,ii) do j=1,lon2 do ii=1,nsubwin do ic3=1,nc3d @@ -2691,8 +2691,8 @@ subroutine sqrt_beta_s_mult_bundle(grady) call stop2(999) endif -!$omp parallel do schedule(dynamic,1) private(ic3,ic2,k,j,i) ! multiply by sqrt_beta_s +!$omp parallel do schedule(dynamic,1) private(ic3,ic2,k,j,i) do j=1,lon2 do ic3=1,nc3d ! check for ozone and skip if oz_univ_static = true @@ -2770,8 +2770,8 @@ subroutine sqrt_beta_e_mult_cvec(grady) ! Initialize timer call timer_ini('sqrt_beta_e_mult') -!$omp parallel do schedule(dynamic,1) private(nn,k,j,i,ii) ! multiply by sqrt_beta_e +!$omp parallel do schedule(dynamic,1) private(nn,k,j,i,ii) do j=1,grd_ens%lon2 do ii=1,nsubwin do nn=1,n_ens @@ -2834,8 +2834,8 @@ subroutine sqrt_beta_e_mult_bundle(aens) ! Initialize timer call timer_ini('sqrt_beta_e_mult') -!$omp parallel do schedule(dynamic,1) private(nn,k,j,i) ! multiply by sqrt_beta_e +!$omp parallel do schedule(dynamic,1) private(nn,k,j,i) do j=1,grd_ens%lon2 do nn=1,n_ens do k=1,nsig @@ -3430,7 +3430,7 @@ subroutine get_new_alpha_beta(aspect,ng,fmat_out,fmat0_out) end subroutine get_new_alpha_beta -subroutine bkerror_a_en(gradx,grady) +subroutine bkerror_a_en(grady) !$$$ subprogram documentation block ! . . . . ! subprogram: bkerror_a_en copy of bkerror for hybrid ensemble @@ -3444,10 +3444,10 @@ subroutine bkerror_a_en(gradx,grady) ! 2010-05-20 todling update to use bundle ! ! input argument list: -! gradx - input field +! grady - input field ! ! output -! grady - background structure * gradx +! grady - background structure * grady ! ! attributes: ! language: f90 @@ -3463,11 +3463,10 @@ subroutine bkerror_a_en(gradx,grady) implicit none ! Declare passed variables - type(control_vector),intent(inout) :: gradx type(control_vector),intent(inout) :: grady ! Declare local variables - integer(i_kind) ii,nn,ip,istatus + integer(i_kind) ii,ip,istatus if (lsqrtb) then write(6,*)'bkerror_a_en: not for use with lsqrtb' @@ -3483,12 +3482,6 @@ subroutine bkerror_a_en(gradx,grady) write(6,*)'bkerror_a_en: trouble getting pointer to ensemble CV' call stop2(317) endif -!$omp parallel do schedule(dynamic,1) private(nn,ii) - do nn=1,n_ens - do ii=1,nsubwin - grady%aens(ii,nn)%r3(ip)%q=gradx%aens(ii,nn)%r3(ip)%q - enddo - enddo ! multiply by sqrt_beta_e_mult call sqrt_beta_e_mult(grady) diff --git a/src/gsi/intall.f90 b/src/gsi/intall.f90 index 43c2b9f350..0f8faa89f8 100644 --- a/src/gsi/intall.f90 +++ b/src/gsi/intall.f90 @@ -319,7 +319,6 @@ subroutine intall(sval,sbias,rval,rbias) end do end if if(npclen > 0)then - do i=1,npclen rbias%predp(i)=qpred(nsclen+i) end do diff --git a/src/gsi/intrad.f90 b/src/gsi/intrad.f90 index 6d3f0f3307..77e92d4a54 100644 --- a/src/gsi/intrad.f90 +++ b/src/gsi/intrad.f90 @@ -307,7 +307,7 @@ subroutine intrad_(radhead,rval,sval,rpred,spred) real(r_quad),dimension(npred*jpch_rad),intent(inout) :: rpred ! Declare local variables - integer(i_kind) j1,j2,j3,j4,i1,i2,i3,i4,n,k,ic,ix,nn,mm,ncr1,ncr2 + integer(i_kind) i1,i2,i3,i4,n,k,ic,ix,nn,mm,ncr1,ncr2 integer(i_kind) ier,istatus integer(i_kind),dimension(nsig) :: i1n,i2n,i3n,i4n real(r_kind),allocatable,dimension(:):: val @@ -388,20 +388,16 @@ subroutine intrad_(radhead,rval,sval,rpred,spred) !radptr => radhead radptr => radNode_typecast(radhead) do while (associated(radptr)) - j1=radptr%ij(1) - j2=radptr%ij(2) - j3=radptr%ij(3) - j4=radptr%ij(4) + i1n(1) = radptr%ij(1) + i2n(1) = radptr%ij(2) + i3n(1) = radptr%ij(3) + i4n(1) = radptr%ij(4) w1=radptr%wij(1) w2=radptr%wij(2) w3=radptr%wij(3) w4=radptr%wij(4) ! Begin Forward model ! calculate temperature, q, ozone, sst vector at observation location - i1n(1) = j1 - i2n(1) = j2 - i3n(1) = j3 - i4n(1) = j4 do k=2,nsig i1n(k) = i1n(k-1)+latlon11 i2n(k) = i2n(k-1)+latlon11 @@ -409,6 +405,7 @@ subroutine intrad_(radhead,rval,sval,rpred,spred) i4n(k) = i4n(k-1)+latlon11 enddo +!$omp parallel do schedule(dynamic,1) private(k,i1,i2,i3,i4,mm) do k=1,nsig i1 = i1n(k) i2 = i2n(k) @@ -423,7 +420,7 @@ subroutine intrad_(radhead,rval,sval,rpred,spred) w3* sq(i3)+w4* sq(i4) endif if(luseoz)then - tdir(ioz+k)=w1* soz(i1)+w2* soz(i2)+ & + tdir(ioz+k)=w1* soz(i1)+w2* soz(i2)+ & w3* soz(i3)+w4* soz(i4) end if if(lusecw)then @@ -454,19 +451,21 @@ subroutine intrad_(radhead,rval,sval,rpred,spred) tdir(iqs+k)=w1* sqs(i1)+w2* sqs(i2)+ & w3* sqs(i3)+w4* sqs(i4) end if + if(k == 1)then + if(luseu)then + tdir(ius+1)= w1* su(i1) +w2* su(i2)+ & + w3* su(i3) +w4* su(i4) + endif + if(lusev)then + tdir(ivs+1)= w1* sv(i1) +w2* sv(i2)+ & + w3* sv(i3) +w4* sv(i4) + endif + if(lusesst)then + tdir(isst+1)=w1*sst(i1) +w2*sst(i2)+ & + w3*sst(i3) +w4*sst(i4) + end if + end if end do - if(luseu)then - tdir(ius+1)= w1* su(j1) +w2* su(j2)+ & - w3* su(j3) +w4* su(j4) - endif - if(lusev)then - tdir(ivs+1)= w1* sv(j1) +w2* sv(j2)+ & - w3* sv(j3) +w4* sv(j4) - endif - if(lusesst)then - tdir(isst+1)=w1*sst(j1) +w2*sst(j2)+ & - w3*sst(j3) +w4*sst(j4) - end if @@ -476,6 +475,7 @@ subroutine intrad_(radhead,rval,sval,rpred,spred) if (.not. ladtest_obs) then allocate(biasvect(radptr%nchan)) +!$omp parallel do schedule(dynamic,1) private(nn,n,ic1,ix1,val_quad) do nn=1,radptr%nchan ic1=radptr%icx(nn) ix1=(ic1-1)*npred @@ -486,8 +486,8 @@ subroutine intrad_(radhead,rval,sval,rpred,spred) biasvect(nn) = val_quad end do end if - ncr1=0 +!$omp parallel do schedule(dynamic,1) private(nn,ic,ix,k) do nn=1,radptr%nchan ic=radptr%icx(nn) ix=(ic-1)*npred @@ -499,8 +499,11 @@ subroutine intrad_(radhead,rval,sval,rpred,spred) do k=1,nsigradjac val(nn)=val(nn)+tdir(k)*radptr%dtb_dvar(k,nn) end do + end do + ncr1=0 ! Include contributions from remaining bias correction terms + do nn=1,radptr%nchan if( .not. ladtest_obs) then if(radptr%use_corr_obs)then val_quad = zero_quad @@ -551,7 +554,6 @@ subroutine intrad_(radhead,rval,sval,rpred,spred) if( .not. ladtest_obs) then if (radptr%use_corr_obs) then - ncr2 = 0 ncr1 = 0 do mm=1,radptr%nchan ncr1 = ncr1 + mm @@ -566,17 +568,15 @@ subroutine intrad_(radhead,rval,sval,rpred,spred) if(radptr%luse)then if(radptr%use_corr_obs)then - do mm=1,radptr%nchan - ic1=radptr%icx(mm) - ix1=(ic1-1)*npred + do nn=1,radptr%nchan + ix=(radptr%icx(nn)-1)*npred do n=1,npred - rpred(ix1+n)=rpred(ix1+n)+biasvect(mm)*radptr%pred(n,mm) + rpred(ix+n)=rpred(ix+n)+biasvect(nn)*radptr%pred(n,nn) enddo enddo else do nn=1,radptr%nchan - ic=radptr%icx(nn) - ix=(ic-1)*npred + ix=(radptr%icx(nn)-1)*npred do n=1,npred rpred(ix+n)=rpred(ix+n)+radptr%pred(n,nn)*val(nn) end do @@ -591,6 +591,7 @@ subroutine intrad_(radhead,rval,sval,rpred,spred) ! Begin adjoint if (l_do_adjoint) then +!$omp parallel do schedule(dynamic,1) private(k,nn) do k=1,nsigradjac tval(k)=zero @@ -603,31 +604,34 @@ subroutine intrad_(radhead,rval,sval,rpred,spred) ! Distribute adjoint contributions over surrounding grid points - if(luseu) then - ru(j1)=ru(j1)+w1*tval(ius+1) - ru(j2)=ru(j2)+w2*tval(ius+1) - ru(j3)=ru(j3)+w3*tval(ius+1) - ru(j4)=ru(j4)+w4*tval(ius+1) - endif - if(lusev) then - rv(j1)=rv(j1)+w1*tval(ivs+1) - rv(j2)=rv(j2)+w2*tval(ivs+1) - rv(j3)=rv(j3)+w3*tval(ivs+1) - rv(j4)=rv(j4)+w4*tval(ivs+1) - endif - - if (lusesst) then - rst(j1)=rst(j1)+w1*tval(isst+1) - rst(j2)=rst(j2)+w2*tval(isst+1) - rst(j3)=rst(j3)+w3*tval(isst+1) - rst(j4)=rst(j4)+w4*tval(isst+1) - end if +!$omp parallel do schedule(dynamic,1) private(k,i1,i2,i3,i4,mm) do k=1,nsig i1 = i1n(k) i2 = i2n(k) i3 = i3n(k) i4 = i4n(k) + if(k == 1)then + if(luseu) then + ru(i1)=ru(i1)+w1*tval(ius+1) + ru(i2)=ru(i2)+w2*tval(ius+1) + ru(i3)=ru(i3)+w3*tval(ius+1) + ru(i4)=ru(i4)+w4*tval(ius+1) + endif + if(lusev) then + rv(i1)=rv(i1)+w1*tval(ivs+1) + rv(i2)=rv(i2)+w2*tval(ivs+1) + rv(i3)=rv(i3)+w3*tval(ivs+1) + rv(i4)=rv(i4)+w4*tval(ivs+1) + endif + + if (lusesst) then + rst(i1)=rst(i1)+w1*tval(isst+1) + rst(i2)=rst(i2)+w2*tval(isst+1) + rst(i3)=rst(i3)+w3*tval(isst+1) + rst(i4)=rst(i4)+w4*tval(isst+1) + end if + end if if(luset)then mm=itv+k diff --git a/src/gsi/intw.f90 b/src/gsi/intw.f90 index 06d13cec2e..5f437ac8ef 100644 --- a/src/gsi/intw.f90 +++ b/src/gsi/intw.f90 @@ -114,7 +114,7 @@ subroutine intw_(whead,rval,sval) integer(i_kind) ib,ik ! real(r_kind) penalty real(r_kind) valu,valv,w1,w2,w3,w4,w5,w6,w7,w8 - real(r_kind) gu,gv,wu,wv,ww + real(r_kind) gu,gv,wu,wv,ww,con real(r_kind) cg_w,p0,gradu,gradv,wnotgross,wgross,term,w_pg real(r_kind),pointer,dimension(:) :: su,sv real(r_kind),pointer,dimension(:) :: ru,rv @@ -198,27 +198,29 @@ subroutine intw_(whead,rval,sval) gradu = valu*wptr%raterr2*wptr%err2 gradv = valv*wptr%raterr2*wptr%err2 else if (njqc .and. wptr%jb > tiny_r_kind .and. wptr%jb <10.0_r_kind) then - valu=sqrt(two*wptr%jb)*tanh(sqrt(wptr%err2)*valu/sqrt(two*wptr%jb)) - valv=sqrt(two*wptr%jb)*tanh(sqrt(wptr%err2)*valv/sqrt(two*wptr%jb)) - gradu = valu*wptr%raterr2*sqrt(wptr%err2) - gradv = valv*wptr%raterr2*sqrt(wptr%err2) + con=sqrt(wptr%err2) + valu=sqrt(two*wptr%jb)*tanh(con*valu/sqrt(two*wptr%jb)) + valv=sqrt(two*wptr%jb)*tanh(con*valv/sqrt(two*wptr%jb)) + gradu = valu*wptr%raterr2*con + gradv = valv*wptr%raterr2*con else if (nvqc .and. wptr%ib >0) then ib=wptr%ib ik=wptr%ik - ww=valu*sqrt(wptr%err2) + con=sqrt(wptr%err2) + ww=valu*con if(hub_norm) then call vqch(ib,ik,ww,gu,wu) else call vqcs(ib,ik,ww,gu,wu) endif - gradu =wu*ww*sqrt(wptr%err2)*wptr%raterr2 - ww=valv*sqrt(wptr%err2) + gradu =wu*ww*con*wptr%raterr2 + ww=valv*con if(hub_norm) then call vqch(ib,ik,ww,gv,wv) else call vqcs(ib,ik,ww,gv,wv) endif - gradv =wv*ww*sqrt(wptr%err2)*wptr%raterr2 + gradv =wv*ww*con*wptr%raterr2 else gradu = valu*wptr%raterr2*wptr%err2 gradv = valv*wptr%raterr2*wptr%err2 diff --git a/src/gsi/jfunc.f90 b/src/gsi/jfunc.f90 index 94a2272851..1f1adc5702 100644 --- a/src/gsi/jfunc.f90 +++ b/src/gsi/jfunc.f90 @@ -28,8 +28,6 @@ module jfunc ! - defer GMAO diurnal bias correction changes. ! 2008-12-01 todling - bring in Tremolet's changes ! 2009-06-01 pondeca,sato - add tsensible initialization. used in 2dvar mode only -! 2009-06-01 pondeca - add lgschmidt initalization. this variable controls the B-norm -! re-orthogonalization of the gradx vectors in 2dvar mode ! 2010-02-20 parrish - add change to get correct nval_len when using hybrid ensemble with dual resolution. ! 2010-02-20 zhu - add nrf_levb and nrf_leve ! 2010-03-23 derber - remove rhgues (not used) @@ -56,6 +54,7 @@ module jfunc ! anav_info - control variables information ! sub read_guess_solution - read guess solution ! sub write_guess_solution - write guess solution +! sub strip1 - strip off halo from subdomain arrays for 1 field ! sub strip2 - strip off halo from subdomain arrays ! sub set_pointer - set location indices for components of vectors ! @@ -82,14 +81,12 @@ module jfunc ! def nvals_levs - number of 2d (x/y) state-vector variables ! def nvals_len - number of 2d state-vector variables * subdomain size (with buffer) ! def nval_levs - number of 2d (x/y) control-vector variables +! def nval_levs_ens - number of 2d (x/y) control-vector variables including +! ensembles ! def nval_len - number of 2d control-vector variables * subdomain size (with buffer) ! def print_diag_pcg - option for turning on GMAO diagnostics in pcgsoi ! def tsensible - option to use sensible temperature as the control variable. applicable ! to the 2dvar mode only -! def lgschmidt - option to re-biorthogonalyze the gradx and grady vectors during the -! inner iteration using the modified gram-schmidt method. useful for -! estimating the analysis error via the projection method. -! ! def ntracer - total number of tracer variables ! def nrft - total number of time tendencies for upper level control variables ! def nrft_ - order of time tendencies for 3d control variables @@ -124,6 +121,7 @@ module jfunc public :: destroy_jfunc public :: read_guess_solution public :: write_guess_solution + public :: strip1 public :: strip2 public :: set_pointer public :: set_sqrt_2dsize @@ -131,22 +129,22 @@ module jfunc public :: nrclen,npclen,nsclen,ntclen,qoption,nval_lenz,tendsflag,tsensible,cwoption,varcw public :: switch_on_derivatives,jiterend,jiterstart,jiter,iter,niter,miter public :: diurnalbc,bcoption,biascor,nval2d,xhatsave,first - public :: factqmax,factqmin,clip_supersaturation,last,yhatsave,nvals_len,nval_levs,iout_iter,nclen + public :: factqmax,factqmin,clip_supersaturation,last,yhatsave,nvals_len,nval_levs,nval_levs_ens,iout_iter,nclen public :: factql,factqi,factqr,factqs,factqg - public :: niter_no_qc,print_diag_pcg,lgschmidt,penorig,gnormorig,iguess + public :: niter_no_qc,print_diag_pcg,penorig,gnormorig,iguess public :: factg,factv,factp,factl,R_option,factw10m,facthowv,factcldch,diag_precon,step_start public :: pseudo_q2 public :: varq public :: cnvw_option - logical first,last,switch_on_derivatives,tendsflag,print_diag_pcg,tsensible,lgschmidt,diag_precon + logical first,last,switch_on_derivatives,tendsflag,print_diag_pcg,tsensible,diag_precon logical clip_supersaturation,R_option logical pseudo_q2 logical cnvw_option integer(i_kind) iout_iter,miter,iguess,nclen,qoption,cwoption integer(i_kind) jiter,jiterstart,jiterend,iter integer(i_kind) nvals_len,nvals_levs - integer(i_kind) nval_len,nval_lenz,nval_levs + integer(i_kind) nval_len,nval_lenz,nval_levs,nval_levs_ens integer(i_kind) nclen1,nclen2,nrclen,nsclen,npclen,ntclen integer(i_kind) nval2d,nclenz @@ -198,7 +196,6 @@ subroutine init_jfunc tendsflag=.false. print_diag_pcg=.false. tsensible=.false. - lgschmidt=.false. diag_precon=.false. step_start=1.e-4_r_kind R_option=.false. @@ -358,7 +355,7 @@ subroutine destroy_jfunc return end subroutine destroy_jfunc - subroutine read_guess_solution(dirx,diry,mype) + subroutine read_guess_solution(diry,mype,success) !$$$ subprogram documentation block ! . . . . ! subprogram: read_guess_solution @@ -376,10 +373,11 @@ subroutine read_guess_solution(dirx,diry,mype) ! ! input argument list: ! mype - mpi task id -! dirx,diry - dynamic components, e.g. %values(:), must be already allocated +! diry - dynamic components, e.g. %values(:), must be already allocated ! ! output argument list: -! dirx,diry +! diry +! success - logical flag determining if input file was successfully read ! ! attributes: ! language: f90 @@ -395,16 +393,17 @@ subroutine read_guess_solution(dirx,diry,mype) implicit none integer(i_kind) ,intent(in ) :: mype - type(control_vector),intent(inout) :: dirx,diry + type(control_vector),intent(inout) :: diry + logical ,intent(inout) :: success integer(i_kind) i,k,mm1,myper,kk,i1,i2 - integer(i_kind) nlatg,nlong,nsigg,nrcleng + integer(i_kind) nlatg,nlong,nsigg + integer(i_kind) nxval,nxrclen,nxsclen,nxpclen integer(i_kind),dimension(5):: iadateg - real(r_single),dimension(max(iglobal,itotsub)):: fieldx,fieldy - real(r_single),dimension(nlat,nlon):: xhatsave_g,yhatsave_g - real(r_single),dimension(nclen):: xhatsave_r4,yhatsave_r4 + real(r_single),dimension(max(iglobal,itotsub)):: fieldy + real(r_single),dimension(nlat,nlon):: yhatsave_g + real(r_single),dimension(nclen):: yhatsave_r4 - jiterstart = 1 mm1=mype+1 myper=0 @@ -415,60 +414,62 @@ subroutine read_guess_solution(dirx,diry,mype) nlatg=0 nlong=0 nsigg=0 - nrcleng=0 - read(12,end=1234)iadateg,nlatg,nlong,nsigg,nrcleng + read(12,end=1234)iadateg,nlatg,nlong,nsigg,nxval,nxrclen,nxsclen,nxpclen if(iadate(1) == iadateg(1) .and. iadate(2) == iadate(2) .and. & iadate(3) == iadateg(3) .and. iadate(4) == iadateg(4) .and. & iadate(5) == iadateg(5) .and. nlat == nlatg .and. & + nxval == nval_levs_ens .and. & nlon == nlong .and. nsig == nsigg) then if(mype == 0) write(6,*)'READ_GUESS_SOLUTION: read guess solution for ',& - iadateg,nlatg,nlong,nsigg,nrcleng - jiterstart=0 + iadateg,nlatg,nlong,nsigg,nxval,nxrclen,nxsclen,nxpclen ! Let all tasks read gesfile_in to pick up bias correction (second read) + if(nxsclen+nxpclen > 0)then + if ( nxsclen==nsclen .and. nxpclen == npclen) then + read(12,end=1236) (yhatsave_r4(i),i=nclen1+1,nclen-ntclen) + do i=nclen+1,nclen-ntclen + diry%values(i)=real(yhatsave_r4(i),r_kind) + end do + else + read(12) + if(mype == 0) then + write(6,*) 'READ_GUESS_SOLUTION: INCOMPATABLE RADIANCE COMPONENT in GESFILE, gesfile_in' + write(6,*) 'READ_GUESS_SOLUTION: nrclen: input,current=',nxsclen,nsclen,nxpclen,npclen + write(6,*) 'READ_GUESS_SOLUTION: ignoring previous radiance guess' + endif + endif + else + if(mype == 0) write(6,*) ' No bias correction in restart file ' + endif ! Loop to read input guess fields. After reading in each field & level, ! scatter the grid to the appropriate location in the xhat and yhatsave ! arrays. - do k=1,nval_levs - read(12,end=1236) xhatsave_g,yhatsave_g - do kk=1,itotsub - i1=ltosi_s(kk); i2=ltosj_s(kk) - fieldx(kk)=xhatsave_g(i1,i2) - fieldy(kk)=yhatsave_g(i1,i2) - end do + do k=1,nval_levs_ens + if(mype == myper)then + read(12,end=1236) yhatsave_g + do kk=1,itotsub + i1=ltosi_s(kk); i2=ltosj_s(kk) + fieldy(kk)=yhatsave_g(i1,i2) + end do + end if i=(k-1)*latlon11 + 1 - call mpi_scatterv(fieldx,ijn_s,displs_s,mpi_real4,& - xhatsave_r4(i),ijn_s(mm1),mpi_real4,myper,mpi_comm_world,ierror) call mpi_scatterv(fieldy,ijn_s,displs_s,mpi_real4,& yhatsave_r4(i),ijn_s(mm1),mpi_real4,myper,mpi_comm_world,ierror) - end do !end do over nval_levs + end do !end do over val_levs_ens` -! Read radiance and precipitation bias correction terms - if ( nrcleng==nrclen ) then - read(12,end=1236) (xhatsave_r4(i),i=nclen1+1,nclen),(yhatsave_r4(i),i=nclen1+1,nclen) - else - if(mype == 0) then - write(6,*) 'READ_GUESS_SOLUTION: INCOMPATABLE RADIANCE COMPONENT in GESFILE, gesfile_in' - write(6,*) 'READ_GUESS_SOLUTION: nrclen: input,current=',nrcleng,nrclen - write(6,*) 'READ_GUESS_SOLUTION: ignoring previous radiance guess' - endif - do i=nclen1+1,nclen - xhatsave_r4(i)=zero - yhatsave_r4(i)=zero - end do - endif - do i=1,nclen - dirx%values(i)=real(xhatsave_r4(i),r_kind) + do i=1,nclen1 diry%values(i)=real(yhatsave_r4(i),r_kind) end do + success = .true. else + read(12) if(mype == 0) then write(6,*) 'READ_GUESS_SOLUTION: INCOMPATABLE GUESS FILE, gesfile_in' write(6,*) 'READ_GUESS_SOLUTION: iguess,iadate,iadateg=',iguess,iadate,iadateg write(6,*) 'READ_GUESS_SOLUTION: nlat,nlatg,nlon,nlong,nsig,nsigg=',& - nlat,nlatg,nlon,nlong,nsig,nsigg + nlat,nlatg,nlon,nlong,nsig,nsigg,nval_levs_ens,nxval end if endif close(12) @@ -535,63 +536,118 @@ subroutine write_guess_solution(mype) integer(i_kind),intent(in ) :: mype integer(i_kind) i,j,k,mm1,mypew,kk,i1,i2,ie,is - real(r_single),dimension(lat1,lon1,2):: field - real(r_single),dimension(max(iglobal,itotsub)):: fieldx,fieldy - real(r_single),dimension(nlat,nlon):: xhatsave_g,yhatsave_g - real(r_single),dimension(nrclen):: xhatsave4,yhatsave4 + real(r_single),dimension(lat1,lon1):: field + real(r_single),dimension(max(iglobal,itotsub)):: fieldy + real(r_single),dimension(nlat,nlon):: yhatsave_g + real(r_single),dimension(nrclen):: yhatsave4 + integer :: nxr,nxs,nxp + logical :: writebias mm1=mype+1 mypew=0 + writebias=.false. + if(writebias)then + nxr=nrclen + nxs=nsclen + nxp=npclen + else + nxr=0 + nxs=0 + nxp=0 + end if ! Write header record to output file if (mype==mypew) then open(51,file='gesfile_out',form='unformatted') - write(51) iadate,nlat,nlon,nsig,nrclen + write(51) iadate,nlat,nlon,nsig,nval_levs_ens,nxr,nxs,nxp endif +! Write radiance and precipitation bias correction terms to output file + if (mype==mypew .and. nsclen+npclen > 0 .and. writebias) then + do i=1,nsclen+npclen + yhatsave4(i)=yhatsave%values(nclen1+i) + end do + write(51) (yhatsave4(i),i=1,nsclen+npclen) + end if + ! Loop over levels. Gather guess solution and write to output - do k=1,nval_levs + do k=1,nval_levs_ens ie=(k-1)*latlon11 + 1 is=ie+latlon11 - call strip2(xhatsave%values(ie:is),yhatsave%values(ie:is),field) - call mpi_gatherv(field(1,1,1),ijn(mm1),mpi_real4,& - fieldx,ijn,displs_g,mpi_real4,mypew,& - mpi_comm_world,ierror) - call mpi_gatherv(field(1,1,2),ijn(mm1),mpi_real4,& + call strip1(yhatsave%values(ie:is),field) + call mpi_gatherv(field,ijn(mm1),mpi_real4,& fieldy,ijn,displs_g,mpi_real4,mypew,& mpi_comm_world,ierror) + if(mype == mypew)then ! Transfer to global arrays - do j=1,nlon - do i=1,nlat - xhatsave_g(i,j)=zero - yhatsave_g(i,j)=zero + do j=1,nlon + do i=1,nlat + yhatsave_g(i,j)=zero + end do + end do + do kk=1,iglobal + i1=ltosi(kk); i2=ltosj(kk) + yhatsave_g(i1,i2)=fieldy(kk) end do - end do - do kk=1,iglobal - i1=ltosi(kk); i2=ltosj(kk) - xhatsave_g(i1,i2)=fieldx(kk) - yhatsave_g(i1,i2)=fieldy(kk) - end do ! Write level record - if (mype==mypew) write(51) xhatsave_g,yhatsave_g - end do !end do over nval_levs + write(51) yhatsave_g + end if + end do !end do over nval_levs_ens -! Write radiance and precipitation bias correction terms to output file - if (mype==mypew) then - do i=1,nrclen - xhatsave4(i)=xhatsave%values(nclen1+i) - yhatsave4(i)=yhatsave%values(nclen1+i) - end do - write(51) (xhatsave4(i),i=1,nrclen),(yhatsave4(i),i=1,nrclen) + if(mype==mypew)then close(51) write(6,*)'WRITE_GUESS_SOLUTION: write guess solution for ',& - iadate,nlat,nlon,nsig,nrclen + iadate,nlat,nlon,nsig,nxr,nxs,nxp endif return end subroutine write_guess_solution + subroutine strip1(field_in1,field_out) +!$$$ subprogram documentation block +! . . . . +! subprogram: strip1 +! prgmmr: treadon org: np23 date: 2003-11-24 +! +! abstract: strip off halo from two subdomain arrays & combine into +! single output array +! +! program history log: +! 2003-11-24 treadon +! 2004-05-18 kleist, documentation +! 2008-05-12 safford - rm unused uses +! +! input argument list: +! field_in1 - subdomain field with halo +! +! output argument list: +! field_out - combined subdomain fields with halo stripped +! +! attributes: +! language: f90 +! machine: ibm rs/6000 sp +! +!$$$ + use kinds, only: r_single + use gridmod, only: lat1,lon1,lat2,lon2 + implicit none + + real(r_single),dimension(lat1,lon1),intent( out) :: field_out + real(r_kind) ,dimension(lat2,lon2),intent(in ) :: field_in1 + + integer(i_kind) i,j,jp1 + + + do j=1,lon1 + jp1 = j+1 + do i=1,lat1 + field_out(i,j)=field_in1(i+1,jp1) + end do + end do + + return + end subroutine strip1 subroutine strip2(field_in1,field_in2,field_out) !$$$ subprogram documentation block @@ -691,6 +747,7 @@ subroutine set_pointer nval_len=nval_levs*latlon11 if(l_hyb_ens) then nval_len=nval_len+n_ens*nsig*grd_ens%latlon11 + nval_levs_ens=nval_levs+n_ens*nsig end if nsclen=npred*jpch_rad npclen=npredp*npcptype diff --git a/src/gsi/pcgsoi.f90 b/src/gsi/pcgsoi.f90 index a628f17eb5..bce3aa0835 100644 --- a/src/gsi/pcgsoi.f90 +++ b/src/gsi/pcgsoi.f90 @@ -127,13 +127,12 @@ subroutine pcgsoi() use obsmod, only: destroyobs,oberror_tune,luse_obsdiag use jfunc, only: iter,jiter,jiterstart,niter,miter,iout_iter,& nclen,penorig,gnormorig,xhatsave,yhatsave,& - iguess,read_guess_solution,diag_precon,step_start, & - niter_no_qc,print_diag_pcg,lgschmidt + iguess,read_guess_solution, & + niter_no_qc,print_diag_pcg use gsi_4dvar, only: nobs_bins, nsubwin, l4dvar, iwrtinc, ladtest, & iorthomax - use gridmod, only: twodvar_regional + use gridmod, only: twodvar_regional,periodic use constants, only: zero,one,tiny_r_kind - use anberror, only: anisotropic use mpimod, only: mype use mpl_allreducemod, only: mpl_allreduce use intallmod, only: intall @@ -145,13 +144,11 @@ subroutine pcgsoi() use state_vectors, only : allocate_state,deallocate_state,& prt_state_norms,inquire_state use bias_predictors, only: allocate_preds,deallocate_preds,predictors,assignment(=) + use anberror, only: anisotropic use bias_predictors, only: update_bias_preds use xhat_vordivmod, only : xhat_vordiv_init, xhat_vordiv_calc, xhat_vordiv_clean use timermod, only: timer_ini,timer_fnl - use projmethod_support, only: init_mgram_schmidt, & - mgram_schmidt,destroy_mgram_schmidt - use hybrid_ensemble_parameters,only : l_hyb_ens,aniso_a_en,ntlevs_ens - use hybrid_ensemble_isotropic, only: bkerror_a_en + use hybrid_ensemble_parameters,only : l_hyb_ens,ntlevs_ens use gsi_bundlemod, only : gsi_bundle use gsi_bundlemod, only : self_add,assignment(=) use gsi_bundlemod, only : gsi_bundleprint @@ -159,7 +156,6 @@ subroutine pcgsoi() use rapidrefresh_cldsurf_mod, only: i_gsdcldanal_type use gsi_io, only: verbose use berror, only: vprecond - use stpjomod, only: stpjo_setup @@ -183,7 +179,7 @@ subroutine pcgsoi() real(r_kind),dimension(3):: gnorm real(r_kind) :: zgini,zfini,fjcost(4),fjcostnew(4),zgend,zfend real(r_kind) :: fjcost_e - type(control_vector) :: xhat,gradx,grady,dirx,diry,ydiff,xdiff + type(control_vector) :: gradx,grady,dirx,diry,ydiff,xdiff type(gsi_bundle) :: sval(nobs_bins), rval(nobs_bins) type(gsi_bundle) :: eval(ntlevs_ens) type(gsi_bundle) :: mval(nsubwin) @@ -193,7 +189,7 @@ subroutine pcgsoi() type(control_vector), allocatable, dimension(:) :: cglworkhat integer(i_kind) :: iortho logical :: print_verbose - logical:: lanlerr + logical :: lanlerr,read_success ! Step size diagnostic strings data step /'good', 'SMALL'/ @@ -210,16 +206,17 @@ subroutine pcgsoi() ! Set constants. Initialize variables. restart=.false. - if (jiter==0 .and. (iguess==1 .or. iguess==2)) restart=.true. + if (jiter==jiterstart .and. (iguess==1 .or. iguess==2)) restart=.true. pennorm=10.e50_r_double iout_6=.true. if (iout_iter==6) iout_6=.false. - stp=step_start - if(diag_precon)stp=one + stp=one small_step=1.e-2_r_kind*stp end_iter=.false. llouter=.false. + gnorm=zero gsave=zero + read_success=.false. ! Convergence criterion needs to be relaxed a bit for anisotropic mode, @@ -236,7 +233,6 @@ subroutine pcgsoi() call init_ if(print_diag_pcg)call prt_guess('guess') - if ( lanlerr .and. lgschmidt ) call init_mgram_schmidt nlnqc_iter=.false. call stpjo_setup(nobs_bins) @@ -252,6 +248,10 @@ subroutine pcgsoi() cglworkhat(ii)=zero END DO end if + do ii=1,nobs_bins + sval(ii)=zero + end do + sbias=zero ! Perform inner iteration inner_iteration: do iter=0,niter(jiter) @@ -266,39 +266,16 @@ subroutine pcgsoi() varqc_iter=one endif end if - +! 1. Calculate gradient do ii=1,nobs_bins rval(ii)=zero end do + rbias=zero gradx=zero - llprt=(mype==0).and.(iter<=1) - -! Convert from control space directly to physical -! space for comparison with obs. - call control2state(xhat,mval,sbias) - if (l4dvar) then - if (l_hyb_ens) then - call ensctl2state(xhat,mval(1),eval) - mval(1)=eval(1) - end if - -! Perform test of AGCM TLM and ADM - call gsi_4dcoupler_grtests(mval,sval,nsubwin,nobs_bins) -! Run TL model to fill sval - call model_tl(mval,sval,llprt) - else - if (l_hyb_ens) then - call ensctl2state(xhat,mval(1),eval) - do ii=1,nobs_bins - sval(ii)=eval(ii) - end do - else - do ii=1,nobs_bins - sval(ii)=mval(1) - end do - end if - end if + llprt=(mype==0).and.(iter<=1) +! Control to state +! call c2s(xhat,sval,sbias,llprt,.true.) if (iter<=1 .and. print_diag_pcg) then do ii=1,nobs_bins @@ -315,35 +292,8 @@ subroutine pcgsoi() enddo endif -! Adjoint of convert control var to physical space - if (l4dvar) then -! Run adjoint model - call model_ad(mval,rval,llprt) - - if (l_hyb_ens) then - eval(1)=mval(1) - call ensctl2state_ad(eval,mval(1),gradx) - end if - else - -! Convert to control space directly from physical space. - if (l_hyb_ens) then - do ii=1,nobs_bins - eval(ii)=rval(ii) - end do - call ensctl2state_ad(eval,mval(1),gradx) - else - mval(1)=rval(1) - if (nobs_bins > 1 ) then - do ii=2,nobs_bins - call self_add(mval(1),rval(ii)) - enddo - end if - end if - - end if - call control2state_ad(mval,rbias,gradx) -! End adjoint of convert control var to physical space +! Adjoint of control to state + call c2s_ad(gradx,rval,rbias,llprt) ! Print initial Jo table if (iter==0 .and. print_diag_pcg .and. luse_obsdiag) then @@ -356,6 +306,7 @@ subroutine pcgsoi() do i=1,nclen gradx%values(i)=gradx%values(i)+yhatsave%values(i) end do +! End of gradient calculation ! Re-orthonormalization if requested if(iorthomax>0) then @@ -370,26 +321,8 @@ subroutine pcgsoi() end if end if -! Multiply by background error - if(anisotropic) then - call anbkerror(gradx,grady) - if(lanlerr .and. lgschmidt) call mgram_schmidt(gradx,grady) - else - call bkerror(gradx,grady) - end if - -! If hybrid ensemble run, then multiply ensemble control variable a_en -! by its localization correlation - if(l_hyb_ens) then - if(aniso_a_en) then - ! call anbkerror_a_en(gradx,grady) ! not available yet - write(6,*)' ANBKERROR_A_EN not written yet, program stops' - stop - else - call bkerror_a_en(gradx,grady) - end if - - end if +! 2. Multiply by background error + call multb(lanlerr,gradx,grady) if(iorthomax>0) then ! save gradients @@ -402,14 +335,11 @@ subroutine pcgsoi() end if end if - if (iter==0 .and. print_diag_pcg) then - call prt_control_norms(grady,'grady') - endif -! Calculate new norm of gradients - if (iter>0) gsave=gnorm(3) +! 3. Calculate new norm of gradients and factors going into b calculation dprod(1) = qdot_prod_sub(gradx,grady) - if(diag_precon)then + if(iter > 0)then + gsave=gnorm(3) if (lanlerr) then ! xdiff used as a temporary array do i=1,nclen @@ -438,99 +368,72 @@ subroutine pcgsoi() gnorm(3)=dprod(4) end if else - if (lanlerr) then - call mpl_allreduce(1,qpvals=dprod) - gnorm(2)=dprod(1) - gnorm(3)=dprod(1) - else - do i=1,nclen - xdiff%values(i)=gradx%values(i)-xdiff%values(i) - ydiff%values(i)=grady%values(i)-ydiff%values(i) - end do - dprod(2) = qdot_prod_sub(xdiff,grady) - dprod(3) = qdot_prod_sub(ydiff,gradx) - call mpl_allreduce(3,qpvals=dprod) -! Two dot products in gnorm(2) should be same, but are slightly -! different due to round off, so use average. - gnorm(2)=0.5_r_quad*(dprod(2)+dprod(3)) - gnorm(3)=dprod(1) - end if +! xdiff used as a temporary array + do i=1,nclen + xdiff%values(i)=vprecond(i)*gradx%values(i) + end do + dprod(2) = qdot_prod_sub(xdiff,grady) + call mpl_allreduce(2,qpvals=dprod) + if(print_diag_pcg) call prt_control_norms(grady,'grady') + gnorm(3)=dprod(2) + end if gnorm(1)=dprod(1) + if(mype == 0)write(iout_iter,*)'Minimization iteration',iter + +! 4. Calculate b and new search direction b=zero - if (gsave>1.e-16_r_kind .and. iter>0) b=gnorm(2)/gsave - if(mype == 0)write(iout_iter,*)'Minimization iteration',iter - if (b7.0_r_kind) then - if (mype==0) then - if (iout_6) write(6,105) gnorm(2),gsave,b - write(iout_iter,105) gnorm(2),gsave,b + if (.not. restart .or. iter > 0) then + if (gsave>1.e-16_r_kind .and. iter>0) b=gnorm(2)/gsave + if (b7.0_r_kind) then + if (mype==0) then + if (iout_6) write(6,105) gnorm(2),gsave,b + write(iout_iter,105) gnorm(2),gsave,b + endif + b=zero endif - b=zero - endif - if (mype==0 .and. print_verbose) write(6,888)'pcgsoi: gnorm(1:3),b=',gnorm,b + if (mype==0 .and. print_verbose) write(6,888)'pcgsoi: gnorm(1:3),b=',gnorm,b + if(read_success .and. iter == 1)b=zero -! Calculate new search direction - if (.not. restart) then - if(.not. lanlerr)then + if (.not. lanlerr) then do i=1,nclen xdiff%values(i)=gradx%values(i) ydiff%values(i)=grady%values(i) end do end if - if(diag_precon)then - do i=1,nclen - dirx%values(i)=-vprecond(i)*grady%values(i)+b*dirx%values(i) - diry%values(i)=-vprecond(i)*gradx%values(i)+b*diry%values(i) - end do - else - do i=1,nclen - dirx%values(i)=-grady%values(i)+b*dirx%values(i) - diry%values(i)=-gradx%values(i)+b*diry%values(i) - end do - end if +! Calculate new search direction + do i=1,nclen + dirx%values(i)=-vprecond(i)*grady%values(i)+b*dirx%values(i) + diry%values(i)=-vprecond(i)*gradx%values(i)+b*diry%values(i) + end do else ! If previous solution available, transfer into local arrays. - if( .not. lanlerr)then - xdiff=zero - ydiff=zero - end if - call read_guess_solution(dirx,diry,mype) - stp=one +! Fill with grady first so that if we read in part of diry there is something + do i=1,nclen + diry%values(i)=-vprecond(i)*gradx%values(i) + end do + call read_guess_solution(diry,mype,read_success) +! Multiply by background error + call multb(lanlerr,diry,dirx) endif +! 5. Calculate stepsize and update solution ! Convert search direction from control space to physical space - call control2state(dirx,mval,rbias) - if (l4dvar) then - if (l_hyb_ens) then - call ensctl2state(dirx,mval(1),eval) - mval(1)=eval(1) - end if - - call model_tl(mval,rval,llprt) - else - - if (l_hyb_ens) then - call ensctl2state(dirx,mval(1),eval) - do ii=1,nobs_bins - rval(ii)=eval(ii) - end do - else - do ii=1,nobs_bins - rval(ii)=mval(1) - end do - end if - - end if + do ii=1,nobs_bins + rval(ii)=zero + end do + rbias=zero + call c2s(dirx,rval,rbias,.false.,.true.) ! Calculate stepsize - call stpcalc(stp,sval,sbias,xhat,dirx,rval,rbias, & + call stpcalc(stp,sval,sbias,dirx,rval,rbias, & diry,penalty,penaltynew,fjcost,fjcostnew,end_iter) if (lanlerr) call writeout_gradients(gradx,grady,niter(jiter),stp,b,mype) -! Diagnostic calculations +! 6. Diagnostic calculations if (iter==0) then if(jiter==jiterstart .or. oberror_tune) then gnormorig=gnorm(1) @@ -572,7 +475,7 @@ subroutine pcgsoi() 9992 format(A,2(1X,I3),4(1X,ES25.18),1x,A6) 9993 format(A,2(1X,I3),2(1X,ES25.18),A1) -! Check for convergence or failure of algorithm +! 7. Check for convergence or failure of algorithm if(gnormx < converge .or. penalty < converge .or. & penx >= pennorm .or. end_iter)then @@ -627,22 +530,11 @@ subroutine pcgsoi() enddo deallocate(cglworkhat) end if - if (lanlerr .and. lgschmidt) call destroy_mgram_schmidt ! Calculate adjusted observation error factor if( oberror_tune .and. (.not.l4dvar) ) then if (mype == 0) write(6,*) 'PCGSOI: call penal for obs perturbation' - call control2state(xhat,mval,sbias) - if (l_hyb_ens) then - call ensctl2state(xhat,mval(1),eval) - do ii=1,nobs_bins - sval(ii)=eval(ii) - end do - else - do ii=1,nobs_bins - sval(ii)=mval(1) - end do - end if +! call c2s(xhat,sval,sbias,.false.,.false.) call penal(sval(1)) xhatsave=zero @@ -656,25 +548,7 @@ subroutine pcgsoi() if (l_tlnmc .and. baldiag_inc) call strong_baldiag_inc(sval,size(sval)) llprt=(mype==0) - call control2state(xhat,mval,sbias) - if (l4dvar) then - if (l_hyb_ens) then - call ensctl2state(xhat,mval(1),eval) - mval(1)=eval(1) - end if - call model_tl(mval,sval,llprt) - else - if (l_hyb_ens) then - call ensctl2state(xhat,mval(1),eval) - do ii=1,nobs_bins - sval(ii)=eval(ii) - end do - else - do ii=1,nobs_bins - sval(ii)=mval(1) - end do - end if - end if +! call c2s(xhat,sval,sbias,llprt,.false.) if(print_diag_pcg)then @@ -686,28 +560,7 @@ subroutine pcgsoi() end do call intall(sval,sbias,rval,rbias) gradx=zero - if (l4dvar) then - call model_ad(mval,rval,llprt) - if (l_hyb_ens) then - eval(1)=mval(1) - call ensctl2state_ad(eval,mval(1),gradx) - end if - else - if (l_hyb_ens) then - do ii=1,nobs_bins - eval(ii)=rval(ii) - end do - call ensctl2state_ad(eval,mval(1),gradx) - else - mval(1)=rval(1) - if (nobs_bins > 1 ) then - do ii=2,nobs_bins - call self_add(mval(1),rval(ii)) - enddo - end if - end if - end if - call control2state_ad(mval,rbias,gradx) + call c2s_ad(gradx,rval,rbias,llprt) ! Add contribution from background term do i=1,nclen @@ -715,24 +568,7 @@ subroutine pcgsoi() end do ! Multiply by background error - if(anisotropic) then - call anbkerror(gradx,grady) - else - call bkerror(gradx,grady) - end if - -! If hybrid ensemble run, then multiply ensemble control variable a_en -! by its localization correlation - if(l_hyb_ens) then - if(aniso_a_en) then - ! call anbkerror_a_en(gradx,grady) ! not available yet - write(6,*)' ANBKERROR_A_EN not written yet, program stops' - stop - else - call bkerror_a_en(gradx,grady) - end if - - end if + call multb(lanlerr,gradx,grady) ! Print final Jo table zgend=dot_product(gradx,grady,r_quad) @@ -866,11 +702,9 @@ subroutine init_ ! !$$$ end documentation block - use jfunc, only: diag_precon implicit none ! Allocate local variables - call allocate_cv(xhat) call allocate_cv(gradx) call allocate_cv(grady) call allocate_cv(dirx) @@ -897,7 +731,6 @@ subroutine init_ diry=zero ydiff=zero xdiff=zero - xhat=zero end subroutine init_ @@ -923,7 +756,6 @@ subroutine clean_ ! !$$$ end documentation block - use jfunc, only: diag_precon use m_obsdiags, only: obsdiags_reset use obsmod, only: destroyobs,lobsdiagsave implicit none @@ -933,7 +765,6 @@ subroutine clean_ if (.not.l4dvar) call obsdiags_reset(obsdiags_keep=lobsdiagsave) ! replacing destroyobs() ! Release state-vector memory - call deallocate_cv(xhat) call deallocate_cv(gradx) call deallocate_cv(grady) call deallocate_cv(dirx) @@ -958,6 +789,241 @@ subroutine clean_ ! call inquire_state end subroutine clean_ + +subroutine periodic_(gradx) +!$$$ subprogram documentation block +! . . . . +! subprogram: periodic_ ensure grad x is periodic +! prgmmr: Todling +! +! abstract: ensure gradx is periodic +! +! program history log: +! 2021-02-02 Derber +! +! input argument list: +! gradx - gradient of x +! +! output argument list: +! gradx - gradient of x +! +! attributes: +! language: f90 +! machine: +! +!$$$ end documentation block + + use gridmod, only: nlat,nlon + use general_commvars_mod, only: s2g_cv + use general_sub2grid_mod, only: general_sub2grid,general_grid2sub + implicit none + + type(control_vector),intent(inout) :: gradx + real(r_kind),dimension(nlat*nlon*s2g_cv%nlevs_alloc)::workcv + +! If dealing with periodic (sub)domain, gather full domain grids, +! account for periodicity, and redistribute to subdomains. This +! only needs to be done when running with a single mpi task and +! then only for array gradx. + do ii=1,nsubwin + call general_sub2grid(s2g_cv,gradx%step(ii)%values,workcv) + call general_grid2sub(s2g_cv,workcv,gradx%step(ii)%values) + end do + + +end subroutine periodic_ + +subroutine multb(lanlerr,vec1,vec2) +!$$$ subprogram documentation block +! . . . . +! subprogram: multb multiply vec1 by background error to equal vec2 +! prgmmr: derber +! +! abstract: multply vec1 by background error +! +! program history log: +! 2021-01-25 derber +! +! input argument list: +! vec1 - input vector +! +! output argument list: +! vec2 - output vector +! +! attributes: +! language: f90 +! machine: +! +!$$$ end documentation block + + use hybrid_ensemble_parameters,only : l_hyb_ens,aniso_a_en + use hybrid_ensemble_isotropic, only: bkerror_a_en + use control_vectors, only: control_vector + implicit none + + type(control_vector),intent(inout) :: vec1 + type(control_vector),intent(inout) :: vec2 + logical ,intent(in ) :: lanlerr + + if(periodic)call periodic_(vec1) +! start by setting vec2=vec1 and then operate on vec2 (unless gram_schmidt) + vec2=vec1 +! Multiply by background error + if(anisotropic) then + call anbkerror(vec2) + else + call bkerror(vec2) + end if + +! If hybrid ensemble run, then multiply ensemble control variable a_en +! by its localization correlation + if(l_hyb_ens) then + if(aniso_a_en) then + ! call anbkerror_a_en(grady) ! not available yet + write(6,*)' ANBKERROR_A_EN not written yet, program stops' + stop + else + call bkerror_a_en(vec2) + end if + + end if + return +end subroutine multb +subroutine c2s(hat,val,bias,llprt,ltest) +!$$$ subprogram documentation block +! . . . . +! subprogram: multb control2state for all options +! prgmmr: derber +! +! abstract: generalized control2state +! +! program history log: +! 2021-01-25 derber +! +! input argument list: +! vec1 - input vector +! +! output argument list: +! vec2 - output vector +! +! attributes: +! language: f90 +! machine: +! +!$$$ end documentation block + + use hybrid_ensemble_parameters,only : l_hyb_ens + use hybrid_ensemble_isotropic, only: bkerror_a_en + use control_vectors, only: control_vector + use bias_predictors, only: predictors + use gsi_bundlemod, only : gsi_bundle,assignment(=) + use gsi_4dvar, only: nobs_bins, nsubwin, l4dvar + use gsi_4dcouplermod, only : gsi_4dcoupler_grtests + implicit none + + type(control_vector) ,intent(inout) :: hat + type(gsi_bundle) ,dimension(nobs_bins),intent(inout) :: val + type(predictors) ,intent(inout) :: bias + logical ,intent(in ) :: llprt,ltest + + +! Convert from control space directly to physical +! space for comparison with obs. + call control2state(hat,mval,bias) + if (l4dvar) then + if (l_hyb_ens) then + call ensctl2state(hat,mval(1),eval) + mval(1)=eval(1) + end if + +! Perform test of AGCM TLM and ADM + if(ltest)call gsi_4dcoupler_grtests(mval,val,nsubwin,nobs_bins) + +! Run TL model to fill val + call model_tl(mval,val,llprt) + else + if (l_hyb_ens) then + call ensctl2state(hat,mval(1),eval) + do ii=1,nobs_bins + val(ii)=eval(ii) + end do + else + do ii=1,nobs_bins + val(ii)=mval(1) + end do + end if + end if + return +end subroutine c2s +subroutine c2s_ad(hat,val,bias,llprt) +!$$$ subprogram documentation block +! . . . . +! subprogram: multb control2state for all options +! prgmmr: derber +! +! abstract: generalized control2state +! +! program history log: +! 2021-01-25 derber +! +! input argument list: +! vec1 - input vector +! +! output argument list: +! vec2 - output vector +! +! attributes: +! language: f90 +! machine: +! +!$$$ end documentation block + + use hybrid_ensemble_parameters,only : l_hyb_ens + use hybrid_ensemble_isotropic, only: bkerror_a_en + use control_vectors, only: control_vector + use bias_predictors, only: predictors + use gsi_bundlemod, only : gsi_bundle,assignment(=) + use gsi_bundlemod, only : self_add + use gsi_4dvar, only: nobs_bins, nsubwin, l4dvar + implicit none + + type(control_vector) ,intent(inout) :: hat + type(gsi_bundle) ,dimension(nobs_bins),intent(inout) :: val + type(predictors) ,intent(inout) :: bias + logical ,intent(in ) :: llprt + + +! Adjoint of convert control var to physical space + if (l4dvar) then +! Run adjoint model + call model_ad(mval,val,llprt) + + if (l_hyb_ens) then + eval(1)=mval(1) + call ensctl2state_ad(eval,mval(1),hat) + end if + else + +! Convert to control space directly from physical space. + if (l_hyb_ens) then + do ii=1,nobs_bins + eval(ii)=val(ii) + end do + call ensctl2state_ad(eval,mval(1),hat) + else + mval(1)=val(1) + if (nobs_bins > 1 ) then + do ii=2,nobs_bins + call self_add(mval(1),val(ii)) + enddo + end if + end if + + end if + call control2state_ad(mval,bias,hat) +! End adjoint of convert control var to physical space + return +end subroutine c2s_ad end subroutine pcgsoi diff --git a/src/gsi/projmethod_support.f90 b/src/gsi/projmethod_support.f90 index 156161dbc8..51b5e78633 100644 --- a/src/gsi/projmethod_support.f90 +++ b/src/gsi/projmethod_support.f90 @@ -1,394 +1,3 @@ -module projmethod_support -!$$$ module documentation block -! . . . . -! module: projmethod_support -! -! abstract: -! -! program history log: -! 2007-10-01 pondeca -! 2010-08-19 lueken - add only to module use -! -! subroutines included: -! sub init_mgram_schmidt -! sub destroy_mgram_schmidt -! sub mgram_schmidt -! -! functions included: -! dplev_mask -! fast_dplev -! dplev5 -! -! variable definition: -! -! attributes: -! language: f90 -! machine: -! -!$$$ end documentation block - -! Uses: - use kinds, only: r_kind - use control_vectors, only: control_vector - - implicit none - -PRIVATE -PUBLIC init_mgram_schmidt, mgram_schmidt, & - destroy_mgram_schmidt - -! Declare local variables - real(r_kind),allocatable,dimension(:,:)::gx,gy - -contains - -subroutine init_mgram_schmidt -!$$$ subprogram documentation block -! . . . . -! subprogram: init_mgram_schmidt -! prgmmr: -! -! abstract: -! -! program history log: -! 2009-09-22 lueken - added subprogram doc block -! -! input argument list: -! -! output argument list: -! -! attributes: -! language: f90 -! machine: -! -!$$$ end documentation block - - use jfunc, only: nclen,jiter,niter - - implicit none - -! Declare passed variables - -! Declare local variables - - allocate(gx(nclen,0:niter(jiter))) - allocate(gy(nclen,0:niter(jiter))) - -end subroutine init_mgram_schmidt - -subroutine destroy_mgram_schmidt -!$$$ subprogram documentation block -! . . . . -! subprogram: destroy_mgram_schmidt -! prgmmr: -! -! abstract: -! -! program history log: -! 2009-09-22 lueken - added subprogram doc block -! -! input argument list: -! -! output argument list: -! -! attributes: -! language: f90 -! machine: -! -!$$$ end documentation block - - implicit none - -! Declare passed variables - - deallocate(gx) - deallocate(gy) - -end subroutine destroy_mgram_schmidt - -subroutine mgram_schmidt(gradx,grady) -!$$$ subprogram documentation block -! . . . . -! subprogram: mgram_schmidt -! prgmmr: pondeca org: np22 date: 2007-08-01 -! -! abstract: apply modified gram-schmidt orthoginalization procedure -! to set of bi-orthogonal vectors -! -! program history log: -! 2007-08-01 pondeca -! -! input argument list: -! gradx - gradient of the cost function w.r.t control variable -! grady - B*(gradx) where B is background error covariance matrix -! -! -! output argument list: -! gradx - modified gradient of the cost function -! grady - modified B*(gradx) -! -! -! remarks: -! -! attributes: -! language: f90 -! machine: ibm RS/6000 SP -! -!$$$ - use kinds, only: i_kind - use jfunc, only: iter,nclen - use constants, only: tiny_r_kind - use mpimod, only: mype - - implicit none - -! Declare passed variables - type(control_vector),intent(inout) :: gradx,grady - -! Declare local variables - integer(i_kind) i,k - real(r_kind) prd0 -!********************************************************************** - - gx(1:nclen,iter)=gradx%values(1:nclen) - gy(1:nclen,iter)=grady%values(1:nclen) - -!==> orthogonalization + renormalization - - do k=1,2 - do i=0,iter-1 - prd0=dplev_mask(gy(1,iter),gx(1,i),mype) - gx(1:nclen,iter)=gx(1:nclen,iter)-gx(1:nclen,i)*prd0 - gy(1:nclen,iter)=gy(1:nclen,iter)-gy(1:nclen,i)*prd0 - enddo - prd0=dplev_mask(gx(1,iter),gy(1,iter),mype) - if (prd0 <= tiny_r_kind) then - if (mype==0) then - print*,'in mgram_schmidt: unable to bi-orthogonalize due to round-off error for iter,k=',iter,k - print*,'in mgram_schmidt: likely to happen when using fast version of inner product' - print*,'in mgram_schmidt: iter,k,prd0=',iter,k,prd0 - endif - return - endif - gx(1:nclen,iter)=gx(1:nclen,iter)/sqrt(prd0) - gy(1:nclen,iter)=gy(1:nclen,iter)/sqrt(prd0) - enddo - -!==> update gradx and grady and put correct B-norm back - - prd0=dplev_mask(gradx%values,grady%values,mype) - gradx%values(1:nclen)=gx(1:nclen,iter)*sqrt(prd0) - grady%values(1:nclen)=gy(1:nclen,iter)*sqrt(prd0) - -contains - -real(r_kind) function dplev_mask(dx,dy,mype) -!$$$ subprogram documentation block -! . . . . -! subprogram: dplev_mask -! prgmmr: -! -! abstract: -! -! program history log: -! 2009-09-22 lueken - added subprogram doc block -! -! input argument list: -! mype -! dx,dy -! -! output argument list: -! -! attributes: -! language: f90 -! machine: -! -!$$$ end documentation block - - use jfunc, only: nval_levs - use gridmod, only: lat2,lon2,twodvar_regional - implicit none - -! Declar passed variables - real(r_kind),dimension(lat2,lon2,nval_levs),intent(in ) :: dx,dy - integer(i_kind) ,intent(in ) :: mype - -! Declare local variables - logical mask(nval_levs) - logical fast - - fast=.false. - mask=.true. -! set fast to .true. for twodvar_regional, -! substantially faster, but no roundoff error reduction and -! results differ for different number of processors. -! if(twodvar_regional) then -!! fast=.true. -! mask(5)=.false. -! mask(6)=.false. -! mask(8)=.false. -! end if - - if(fast) then - dplev_mask=fast_dplev(dx,dy,mask) - else - dplev_mask=dplev5(dx,dy,mype,mask) - end if - -end function dplev_mask - -real(r_kind) function fast_dplev(dx,dy,mask) -!$$$ subprogram documentation block -! . . . . -! subprogram: fast_dplev -! prgmmr: -! -! abstract: -! -! program history log: -! 2009-09-22 lueken - added subprogram doc block -! -! input argument list: -! dx,dy -! mask -! -! output argument list: -! -! attributes: -! language: f90 -! machine: -! -!$$$ end documentation block - - use jfunc, only: nval_levs - use gridmod, only: lat2,lon2 - use mpimod, only: npe,mpi_comm_world,ierror,mpi_rtype - use constants, only: zero - implicit none - -! Declar passed variables - real(r_kind),dimension(lat2,lon2,nval_levs),intent(in ) :: dx,dy - logical ,intent(in ) :: mask(nval_levs) - -! Declare local variables - real(r_kind),dimension(npe):: sumall - - integer(i_kind) i,j,k - real(r_kind) sum - - sum=zero - do k=1,nval_levs - if(.not.mask(k)) cycle - do j=2,lon2-1 - do i=2,lat2-1 - sum=sum+dx(i,j,k)*dy(i,j,k) - end do - end do - end do - - call mpi_allgather(sum,1,mpi_rtype,sumall,1,mpi_rtype,mpi_comm_world,ierror) - fast_dplev=zero - do i=1,npe - fast_dplev=fast_dplev+sumall(i) - end do - -end function fast_dplev - -real(r_kind) function dplev5(dx,dy,mype,mask) -!$$$ subprogram documentation block -! . . . . -! subprogram: dplev calculates dot product for data on subdomain -! prgmmr: derber org: np23 date: 2004-05-13 -! -! abstract: calculates dot product for data on subdomain. Note loops over -! streamfunction, velocity potential, temperature, etc. Also, only -! over interior of subdomain. -! -! program history log: -! 2004-05-13 derber, document -! 2004-07-28 treadon - add only on use declarations; add intent in/out -! 2010-04-01 treadon - move strip to grimod -! 2013-10-25 todling - reposition ltosi and others to commvars -! -! input argument list: -! dx - input vector 1 -! dy - input vector 2 -! mype -! mask -! -! output argument list -! -! attributes: -! language: f90 -! machine: ibm RS/6000 SP -! -!$$$ - - use jfunc, only: nval_levs - use gridmod, only: nlat,nlon,lat2,lon2,lat1,lon1,& - iglobal,ijn,displs_g,strip - use general_commvars_mod, only: ltosi,ltosj - use mpimod, only: mpi_comm_world,ierror,mpi_rtype - use constants, only: zero - implicit none - -! Declare passed variables - real(r_kind),dimension(lat2,lon2,nval_levs),intent(in ) :: dx,dy - integer(i_kind) ,intent(in ) :: mype - logical ,intent(in ) :: mask(nval_levs) - -! Declare local variables - real(r_kind),dimension(lat1*lon1):: zsm - real(r_kind),dimension(iglobal):: work1 - real(r_kind),dimension(lat2,lon2):: sum - real(r_kind),dimension(nlat,nlon):: sumall - - integer(i_kind) i,j,k,mm1 - real(r_kind) e,y,temp - - mm1=mype+1 - sum=zero - do k=1,nval_levs - if(.not.mask(k)) cycle - do j=1,lon2 - do i=1,lat2 - sum(i,j)=sum(i,j)+dx(i,j,k)*dy(i,j,k) - end do - end do - end do - do j=1,lon1*lat1 - zsm(j)=zero - end do - - call strip(sum,zsm) - - call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,& - work1,ijn,displs_g,mpi_rtype,& - mpi_comm_world,ierror) - - do k=1,iglobal - i=ltosi(k) ; j=ltosj(k) - sumall(i,j)=work1(k) - end do - dplev5=zero - e=zero - do j=1,nlon - do i=1,nlat -! Compensated summation version of sum - temp=dplev5 - y=sumall(i,j)+e - dplev5=temp+y - e=(temp-dplev5)+y -! dplev=dplev+sumall(i,j) - end do - end do - -end function dplev5 - -end subroutine mgram_schmidt - -end module projmethod_support - subroutine writeout_gradients(dx,dy,nv,alpha,gamma,mype) !$$$ subprogram documentation block ! . . . . diff --git a/src/gsi/read_atms.f90 b/src/gsi/read_atms.f90 index df2612a420..a0dce48f3a 100644 --- a/src/gsi/read_atms.f90 +++ b/src/gsi/read_atms.f90 @@ -234,7 +234,7 @@ subroutine read_atms(mype,val_tovs,ithin,isfcalc,& ! Set various variables depending on type of data to be read if (obstype /= 'atms') then - write(*,*) 'READ_ATMS called for obstype '//obstype//': RETURNING' + write(6,*) 'READ_ATMS called for obstype '//obstype//': RETURNING' return end if @@ -266,7 +266,7 @@ subroutine read_atms(mype,val_tovs,ithin,isfcalc,& elseif (jsatid == 'n21') then kidsat = 226 else - write(*,*) 'READ_ATMS: Unrecognized value for jsatid '//jsatid//': RETURNING' + write(6,*) 'READ_ATMS: Unrecognized value for jsatid '//jsatid//': RETURNING' return end if @@ -498,7 +498,7 @@ subroutine read_atms(mype,val_tovs,ithin,isfcalc,& num_obs = iob-1 if (num_obs <= 0) then - write(*,*) 'READ_ATMS: No ATMS Data were read in' + write(6,*) 'READ_ATMS: No ATMS Data were read in' return end if @@ -507,14 +507,14 @@ subroutine read_atms(mype,val_tovs,ithin,isfcalc,& ALLOCATE(Relative_Time_In_Seconds(Num_Obs)) ALLOCATE(IScan(Num_Obs)) Relative_Time_In_Seconds = 3600.0_r_kind*T4DV_Save(1:Num_Obs) - write(*,*) 'Calling ATMS_Spatial_Average' + write(6,*) 'Calling ATMS_Spatial_Average' CALL ATMS_Spatial_Average(Num_Obs, NChanl, IFOV_Save(1:Num_Obs), & Relative_Time_In_Seconds, BT_Save(1:nchanl,1:Num_Obs), IScan, IRet) - write(*,*) 'ATMS_Spatial_Average Called with IRet=',IRet + write(6,*) 'ATMS_Spatial_Average Called with IRet=',IRet DEALLOCATE(Relative_Time_In_Seconds) IF (IRet /= 0) THEN - write(*,*) 'Error Calling ATMS_Spatial_Average from READ_ATMS' + write(6,*) 'Error Calling ATMS_Spatial_Average from READ_ATMS' RETURN END IF diff --git a/src/gsi/setupt.f90 b/src/gsi/setupt.f90 index 276659995d..06d50497c3 100644 --- a/src/gsi/setupt.f90 +++ b/src/gsi/setupt.f90 @@ -1007,6 +1007,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav endif my_head%luse = luse(i) my_head%tv_ob = iqtflg + my_head%idx = 0 if (aircraft_t_bc_pof .or. aircraft_t_bc) then effective=upd_pred_t*pof_idx diff --git a/src/gsi/smoothrf.f90 b/src/gsi/smoothrf.f90 index 455ed29636..ec5931db4a 100644 --- a/src/gsi/smoothrf.f90 +++ b/src/gsi/smoothrf.f90 @@ -1171,7 +1171,7 @@ subroutine sqrt_smoothrf(z,work,nlevs) else do j=1,nhscrf - totwgt(j)=sqrt(hswgt(j)*hzscl(j)*hzscl(j)) + totwgt(j)=sqrt(hswgt(j))*hzscl(j) end do ! zero output array @@ -1330,7 +1330,7 @@ subroutine sqrt_smoothrf_ad(z,work,nlevs) if(nrf_var(nvar_id(k))=='sf'.or.nrf_var(nvar_id(k))=='vp')then totwgt(3)=sqrt(half)*totwgt(3) end if - + call sqrt_rfxyyx_ad(zloc,work(1,1,k),ny,nx,ii(1,1,1,k),& jj(1,1,1,k),slw(1,k),totwgt) @@ -1348,7 +1348,7 @@ subroutine sqrt_smoothrf_ad(z,work,nlevs) else do j=1,nhscrf - totwgt(j)=sqrt(hswgt(j)*hzscl(j)*hzscl(j)) + totwgt(j)=sqrt(hswgt(j))*hzscl(j) end do diff --git a/src/gsi/stpcalc.f90 b/src/gsi/stpcalc.f90 index a3ea5e67bd..8b67e35742 100644 --- a/src/gsi/stpcalc.f90 +++ b/src/gsi/stpcalc.f90 @@ -40,7 +40,7 @@ module stpcalcmod contains -subroutine stpcalc(stpinout,sval,sbias,xhat,dirx,dval,dbias, & +subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & diry,penalty,penaltynew,pjcost,pjcostnew,end_iter) !$$$ subprogram documentation block @@ -188,7 +188,6 @@ subroutine stpcalc(stpinout,sval,sbias,xhat,dirx,dval,dbias, & ! input argument list: ! stpinout - guess stepsize ! sval - current solution -! xhat - current solution ! dirx - search direction for x ! diry - search direction for y (B-1 dirx) ! end_iter - end iteration flag @@ -196,7 +195,6 @@ subroutine stpcalc(stpinout,sval,sbias,xhat,dirx,dval,dbias, & ! sbias,dbias ! ! output argument list: -! xhat ! stpinout - final estimate of stepsize ! penalty - penalty current solution ! penaltynew - estimate of penalty for new solution @@ -222,7 +220,7 @@ subroutine stpcalc(stpinout,sval,sbias,xhat,dirx,dval,dbias, & use constants, only: zero,one_quad,zero_quad use gsi_4dvar, only: nobs_bins,ltlint,ibin_anl use jfunc, only: iout_iter,nclen,xhatsave,yhatsave,& - iter + iter,nrclen use jcmod, only: ljcpdry,ljc4tlevs,ljcdfi,ljclimqc use gsi_obOperTypeManager, only: nobs_type => obOper_count use stpjcmod, only: stplimq,stplimg,stplimv,stplimp,stplimw10m,& @@ -247,11 +245,11 @@ subroutine stpcalc(stpinout,sval,sbias,xhat,dirx,dval,dbias, & real(r_kind) ,intent( out) :: penalty,penaltynew real(r_kind) ,intent( out) :: pjcost(4),pjcostnew(4) - type(control_vector),intent(inout) :: xhat type(control_vector),intent(in ) :: dirx,diry - type(gsi_bundle) ,intent(in ) :: sval(nobs_bins) + type(gsi_bundle) ,intent(inout) :: sval(nobs_bins) type(gsi_bundle) ,intent(in ) :: dval(nobs_bins) - type(predictors) ,intent(in ) :: sbias,dbias + type(predictors) ,intent(inout) :: sbias + type(predictors) ,intent(in ) :: dbias ! Declare local parameters @@ -923,9 +921,17 @@ subroutine stpcalc(stpinout,sval,sbias,xhat,dirx,dval,dbias, & endif ! Update solution + do i=1,nrclen + sbias%values(i)=sbias%values(i)+stpinout*dbias%values(i) + end do +!$omp parallel do schedule(dynamic,1) private(i,ii) + do ii=1,nobs_bins + do i=1,sval(ii)%ndim + sval(ii)%values(i)=sval(ii)%values(i)+stpinout*dval(ii)%values(i) + end do + end do !DIR$ IVDEP do i=1,nclen - xhat%values(i)=xhat%values(i)+stpinout*dirx%values(i) xhatsave%values(i)=xhatsave%values(i)+stpinout*dirx%values(i) yhatsave%values(i)=yhatsave%values(i)+stpinout*diry%values(i) end do diff --git a/src/gsi/stpt.f90 b/src/gsi/stpt.f90 index f793c6a214..27f5385ac1 100644 --- a/src/gsi/stpt.f90 +++ b/src/gsi/stpt.f90 @@ -193,10 +193,10 @@ subroutine stpt(thead,dval,xval,out,sges,nstep,rpred,spred) val2=w1*stv(j1)+w2*stv(j2)+w3*stv(j3)+w4*stv(j4)+ & w5*stv(j5)+w6*stv(j6)+w7*stv(j7)+w8*stv(j8) else - val= w1* rt(j1)+w2* rt(j2)+w3* rt(j3)+w4* rt(j4)+ & - w5* rt(j5)+w6* rt(j6)+w7* rt(j7)+w8* rt(j8) - val2=w1* st(j1)+w2* st(j2)+w3* st(j3)+w4* st(j4)+ & - w5* st(j5)+w6* st(j6)+w7* st(j7)+w8* st(j8) + val= w1* rt(j1)+w2* rt(j2)+w3* rt(j3)+w4* rt(j4)+ & + w5* rt(j5)+w6* rt(j6)+w7* rt(j7)+w8* rt(j8) + val2=w1* st(j1)+w2* st(j2)+w3* st(j3)+w4* st(j4)+ & + w5* st(j5)+w6* st(j6)+w7* st(j7)+w8* st(j8) end if ! contribution from bias correction @@ -233,9 +233,9 @@ subroutine stpt(thead,dval,xval,out,sges,nstep,rpred,spred) ts_prime=tt(kk) tg_prime=valsst2+sges(kk)*valsst qs_prime=valq2+sges(kk)*valq - us_prime=valu2+sges(kk)*val - vs_prime=valv2+sges(kk)*val - psfc_prime=val2+sges(1)*val + us_prime=valu2+sges(kk)*valu + vs_prime=valv2+sges(kk)*valv + psfc_prime=valp2+sges(1)*valp tt(kk)=psfc_prime*tptr%tlm_tsfc(1) + tg_prime*tptr%tlm_tsfc(2) + & ts_prime *tptr%tlm_tsfc(3) + qs_prime*tptr%tlm_tsfc(4) + &