From 860a459a2f4cec71733a68c9f7e7b7b1054d52e5 Mon Sep 17 00:00:00 2001 From: Terra Ladwig Date: Tue, 22 Jun 2021 23:00:58 +0000 Subject: [PATCH 1/3] Adding print statements to track/debug setupcldtot execution. --- src/gsi/gsi_cldtotOper.F90 | 2 ++ src/gsi/setupcldtot.F90 | 3 +++ src/gsi/setuprhsall.f90 | 1 + 3 files changed, 6 insertions(+) diff --git a/src/gsi/gsi_cldtotOper.F90 b/src/gsi/gsi_cldtotOper.F90 index 5a63e24765..b0285d1098 100644 --- a/src/gsi/gsi_cldtotOper.F90 +++ b/src/gsi/gsi_cldtotOper.F90 @@ -96,6 +96,8 @@ subroutine setup_(self, lunin, mype, is, nobs, init_pass,last_pass) if(nobs == 0) return ! try data header + write (*,*) 'cldtot_Oper', mype, i_cloud_q_innovation + write (*,*) obstype,isis,nreal,nchanl read(lunin,iostat=ier) obstype,isis,nreal,nchanl if(ier/=0) then call perr(myname_,'read(obstype,...), iostat =',ier) diff --git a/src/gsi/setupcldtot.F90 b/src/gsi/setupcldtot.F90 index fc3978e3d2..c328056ac2 100755 --- a/src/gsi/setupcldtot.F90 +++ b/src/gsi/setupcldtot.F90 @@ -245,6 +245,9 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di ! If require guess vars available, extract from bundle ... call init_vars_ + write (*,*) 'setupcldtot', mype + + allocate(h_bk(nsig)) allocate(t_bk(nsig)) allocate(q_bk(nsig)) diff --git a/src/gsi/setuprhsall.f90 b/src/gsi/setuprhsall.f90 index 2a06e95510..d33b6d9d43 100644 --- a/src/gsi/setuprhsall.f90 +++ b/src/gsi/setuprhsall.f90 @@ -484,6 +484,7 @@ subroutine setuprhsall(ndata,mype,init_pass,last_pass) is_obOper => obOper_create(dtype(is)) if(associated(is_obOper)) then + write(*,*) 'setuprhsall', trim(dtype(is)), is call is_obOper%setup(lunin,mype, is, nsat1(is), init_pass,last_pass) call obOper_destroy(is_obOper) From f95d790d73de60fb35e445feece3fee42d8aaf0b Mon Sep 17 00:00:00 2001 From: Terra Ladwig Date: Thu, 24 Jun 2021 13:57:16 +0000 Subject: [PATCH 2/3] Update to cloud operator, namelist control, and remove debug prints. --- src/gsi/gsdcloudlib_pseudoq_mod.f90 | 33 +++++++++++++++++++++++----- src/gsi/gsi_cldtotOper.F90 | 4 +--- src/gsi/gsimod.F90 | 4 +++- src/gsi/rapidrefresh_cldsurf_mod.f90 | 4 +++- src/gsi/setupcldtot.F90 | 23 ++++++++++--------- src/gsi/setuprhsall.f90 | 5 ++--- 6 files changed, 50 insertions(+), 23 deletions(-) diff --git a/src/gsi/gsdcloudlib_pseudoq_mod.f90 b/src/gsi/gsdcloudlib_pseudoq_mod.f90 index 5afba70452..707be38f99 100644 --- a/src/gsi/gsdcloudlib_pseudoq_mod.f90 +++ b/src/gsi/gsdcloudlib_pseudoq_mod.f90 @@ -32,7 +32,8 @@ module gsdcloudlib_pseudoq_mod contains -SUBROUTINE cloudCover_Surface_col(mype,nsig,& +SUBROUTINE cloudCover_Surface_col(mype,nsig, & + i_cloud_q_innovation,& cld_bld_hgt,h_bk,zh, & NVARCLD_P,ocld,Oelvtn,& wthr_type,pcp_type_obs, & @@ -57,6 +58,7 @@ SUBROUTINE cloudCover_Surface_col(mype,nsig,& ! input argument list: ! mype - processor ID ! nsig - no. of levels +! i_cloud_q_innovation - flag to control building/clearing/both ! cld_bld_hgt - Height below which cloud building is done ! ! h_bk - 3D background height (m) @@ -96,6 +98,7 @@ SUBROUTINE cloudCover_Surface_col(mype,nsig,& integer(i_kind),intent(in) :: mype integer(i_kind),intent(in) :: nsig + integer(i_kind),intent(in) :: i_cloud_q_innovation real(r_kind), intent(in) :: cld_bld_hgt ! ! surface observation @@ -124,7 +127,7 @@ SUBROUTINE cloudCover_Surface_col(mype,nsig,& INTEGER(i_kind) :: k INTEGER(i_kind) :: ic integer(i_kind) :: firstcloud,cl_base_broken_k,obused - integer(i_kind) :: kcld + integer(i_kind) :: kcld,kclr real(r_single) :: underlim REAL(r_kind) :: zdiff REAL(r_kind) :: zlev_clr,cloud_dz,cl_base_ista,betav @@ -133,6 +136,7 @@ SUBROUTINE cloudCover_Surface_col(mype,nsig,& !==================================================================== ! Begin ! + !write(6,*) 'cloudCover_Surface', mype, i_cloud_q_innovation ! set constant names consistent with original RUC code ! vis2qc=-9999.0_r_single @@ -140,6 +144,7 @@ SUBROUTINE cloudCover_Surface_col(mype,nsig,& firstcloud = 0 obused =0 kcld=-9 + kclr=99 ! !***************************************************************** ! analysis of surface/METAR cloud observations @@ -161,6 +166,15 @@ SUBROUTINE cloudCover_Surface_col(mype,nsig,& endif enddo + ! cloud clearing obs + if(i_cloud_q_innovation==20 .or. i_cloud_q_innovation==22) then + do k=3,nsig,5 + if (h_bk(k) < zlev_clr) then + cld_cover_obs(k)=0.0_r_single + endif + enddo + endif + ! -- Now consider non-clear obs ! -------------------------- else @@ -203,10 +217,12 @@ SUBROUTINE cloudCover_Surface_col(mype,nsig,& if(k==8) underlim=95.0_r_kind ! 3000 feet if(k>=9 .and. k= 1.0_r_kind .and. (firstcloud==0 .or. abs(zdiff)= 3) cld_cover_obs(kclr)=0.0_r_single + endif - endif ! end if ocld valid endif ! obused enddo ! end IC loop endif ! end if cloudy ob diff --git a/src/gsi/gsi_cldtotOper.F90 b/src/gsi/gsi_cldtotOper.F90 index b0285d1098..11516daac3 100644 --- a/src/gsi/gsi_cldtotOper.F90 +++ b/src/gsi/gsi_cldtotOper.F90 @@ -96,8 +96,6 @@ subroutine setup_(self, lunin, mype, is, nobs, init_pass,last_pass) if(nobs == 0) return ! try data header - write (*,*) 'cldtot_Oper', mype, i_cloud_q_innovation - write (*,*) obstype,isis,nreal,nchanl read(lunin,iostat=ier) obstype,isis,nreal,nchanl if(ier/=0) then call perr(myname_,'read(obstype,...), iostat =',ier) @@ -109,7 +107,7 @@ subroutine setup_(self, lunin, mype, is, nobs, init_pass,last_pass) diagsave = write_diag(jiter) .and. diag_conv select case(i_cloud_q_innovation) - case(2) + case(20, 21, 22) call setup(self%obsLL(:), self%odiagLL(:), & lunin,mype,bwork,awork(:,iwork),nele,nobs,is,diagsave) diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index af5c77d6c6..73183006f0 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -1201,7 +1201,9 @@ module gsimod ! i_cloud_q_innovation - integer to choose if and how cloud obs are used ! 0= no innovations ! 1= cloud total innovations -! 2= water vapor innovations +! 20= cloud build/clear derived water vapor innovations +! 21= cloud build derived water vapor innovations +! 22= cloud clear derived water vapor innovations ! 3= cloud total & water vapor innovations ! i_ens_mean - integer for setupcldtot behavior ! 0=single model run diff --git a/src/gsi/rapidrefresh_cldsurf_mod.f90 b/src/gsi/rapidrefresh_cldsurf_mod.f90 index cd52133134..1ee35fffba 100644 --- a/src/gsi/rapidrefresh_cldsurf_mod.f90 +++ b/src/gsi/rapidrefresh_cldsurf_mod.f90 @@ -156,7 +156,9 @@ module rapidrefresh_cldsurf_mod ! i_cloud_q_innovation - integer to choose if and how cloud obs are used ! 0= no innovations ! 1= cloud total innovations -! 2= water vapor innovations +! 20= cloud build/clear derived water vapor innovations +! 21= cloud build derived water vapor innovations +! 22= cloud clear derived water vapor innovations ! 3= cloud total & water vapor innovations ! i_ens_mean - integer for setupcldtot behavior ! 0=single model run diff --git a/src/gsi/setupcldtot.F90 b/src/gsi/setupcldtot.F90 index c328056ac2..980dca5314 100755 --- a/src/gsi/setupcldtot.F90 +++ b/src/gsi/setupcldtot.F90 @@ -12,7 +12,7 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di !! . . . . ! subprogram: setupcldtot compute rhs of oi for pseudo moisture observations from ! METAR and Satellite cloud observations -! prgmmr: Ladwag org: GSD date: 2019-06-01 +! prgmmr: Ladwig org: GSD date: 2019-06-01 ! ! abstract: For moisture observations, this routine ! a) reads obs assigned to given mpi task (geographic region), @@ -245,9 +245,6 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di ! If require guess vars available, extract from bundle ... call init_vars_ - write (*,*) 'setupcldtot', mype - - allocate(h_bk(nsig)) allocate(t_bk(nsig)) allocate(q_bk(nsig)) @@ -276,7 +273,7 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di allocate(cdiagbuf(nobs*nsig),rdiagbuf(nreal,nobs*nsig)) rdiagbuf=zero endif - if (i_cloud_q_innovation == 2 .or. i_cloud_q_innovation == 3) then + if (i_cloud_q_innovation .ge. 20 .or. i_cloud_q_innovation == 3) then iip=0 allocate(cdiagbufp(nobs*nsig),rdiagbufp(nreal,nobs*nsig)) cdiagbufp="EMPTY" @@ -464,7 +461,7 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di cycle endif - call cloudCover_surface_col(mype,nsig,cld_bld_hgt,h_bk,z_bk, & + call cloudCover_surface_col(mype,nsig,i_cloud_q_innovation,cld_bld_hgt,h_bk,z_bk, & nvarcld_p,ocld,oelvtn,wthr_type,pcp_type_obs,vis2qc,cld_cover_obs) @@ -519,8 +516,8 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di muse(i)=.true. !******************************************************************************* - if (i_cloud_q_innovation /= 2) then - write(*,*) "Warning - setupcldtot: this code version is only designed for i_cloud_q_innovation == 2" + if (i_cloud_q_innovation .lt. 20 .or. i_cloud_q_innovation .gt. 22 ) then + write(*,*) "Warning - setupcldtot: this code version is only designed for i_cloud_q_innovation == 20,21,22" return else @@ -569,6 +566,8 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di ddiff=qv_ob-q_bk(k) q_build0_count=q_build0_count+1 endif + ! build error = 80% + error=one/(cloudqvis*8.E-01_r_kind) elseif (qob > -0.000001_r_single) then @@ -581,13 +580,16 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di ddiff=qv_ob-q_bk(k) q_clear0_count=q_clear0_count+1 endif + ! clear error = 30% + error=one/(cloudqvis*3.E-01_r_kind) else cycle endif q_obcount=q_obcount+1 - error=one/(cloudqvis*3.E-01_r_kind) + ! all obs errors = 30% + !error=one/(cloudqvis*3.E-01_r_kind) ratio_errors=1.0_r_kind val = error*ddiff @@ -715,7 +717,8 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di !! Write information to diagnostic file if(conv_diagsave)then - if (i_cloud_q_innovation == 2 .and. iip>0) then + if (i_cloud_q_innovation .ge. 20 .and. iip>0) then +! call dtime_show(myname,'diagsave:q',i_q_ob_type) if(netcdf_diag) call nc_diag_write if(binary_diag)then write(7)' q',nchar,nreal,iip,mype,ioff0 diff --git a/src/gsi/setuprhsall.f90 b/src/gsi/setuprhsall.f90 index d33b6d9d43..205e1ad6c2 100644 --- a/src/gsi/setuprhsall.f90 +++ b/src/gsi/setuprhsall.f90 @@ -484,7 +484,6 @@ subroutine setuprhsall(ndata,mype,init_pass,last_pass) is_obOper => obOper_create(dtype(is)) if(associated(is_obOper)) then - write(*,*) 'setuprhsall', trim(dtype(is)), is call is_obOper%setup(lunin,mype, is, nsat1(is), init_pass,last_pass) call obOper_destroy(is_obOper) @@ -540,11 +539,11 @@ subroutine setuprhsall(ndata,mype,init_pass,last_pass) ! luse_obsdiag, sorting could become a problem. Among them, cases of ! l_PBL_pseudo_SurfobsT, l_PBL_pseudo_SurfobsQ, and l_PBL_pseudo_SurfobsUV ! have been fixed since, but it might be better to keep it simple for - ! those applications. The case of i_cloud_q_innovation==2 is new. It is + ! those applications. The case of i_cloud_q_innovation is new. It is ! not sure why it won't work even in case of .not.luse_obsdiag. if(.not.(l_PBL_pseudo_SurfobsT .or. l_PBL_pseudo_SurfobsQ .or. & - l_PBL_pseudo_SurfobsUV .or. (i_cloud_q_innovation==2)) ) then + l_PBL_pseudo_SurfobsUV .or. (i_cloud_q_innovation.ge.1)) ) then call obsdiags_sort() endif From 6af4ef1f1f16f6593571cd7ae5bf0350592e3165 Mon Sep 17 00:00:00 2001 From: Terra Ladwig Date: Thu, 24 Jun 2021 18:47:46 +0000 Subject: [PATCH 3/3] Fix to use conditional symbols, e.g. > or <, to comply with GSI code standards. --- src/gsi/setupcldtot.F90 | 6 +++--- src/gsi/setuprhsall.f90 | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/gsi/setupcldtot.F90 b/src/gsi/setupcldtot.F90 index 980dca5314..22fa31cce1 100755 --- a/src/gsi/setupcldtot.F90 +++ b/src/gsi/setupcldtot.F90 @@ -273,7 +273,7 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di allocate(cdiagbuf(nobs*nsig),rdiagbuf(nreal,nobs*nsig)) rdiagbuf=zero endif - if (i_cloud_q_innovation .ge. 20 .or. i_cloud_q_innovation == 3) then + if (i_cloud_q_innovation >= 20 .or. i_cloud_q_innovation == 3) then iip=0 allocate(cdiagbufp(nobs*nsig),rdiagbufp(nreal,nobs*nsig)) cdiagbufp="EMPTY" @@ -516,7 +516,7 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di muse(i)=.true. !******************************************************************************* - if (i_cloud_q_innovation .lt. 20 .or. i_cloud_q_innovation .gt. 22 ) then + if (i_cloud_q_innovation < 20 .or. i_cloud_q_innovation > 22 ) then write(*,*) "Warning - setupcldtot: this code version is only designed for i_cloud_q_innovation == 20,21,22" return else @@ -717,7 +717,7 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di !! Write information to diagnostic file if(conv_diagsave)then - if (i_cloud_q_innovation .ge. 20 .and. iip>0) then + if (i_cloud_q_innovation >= 20 .and. iip>0) then ! call dtime_show(myname,'diagsave:q',i_q_ob_type) if(netcdf_diag) call nc_diag_write if(binary_diag)then diff --git a/src/gsi/setuprhsall.f90 b/src/gsi/setuprhsall.f90 index 205e1ad6c2..52a7ca5b0b 100644 --- a/src/gsi/setuprhsall.f90 +++ b/src/gsi/setuprhsall.f90 @@ -543,7 +543,7 @@ subroutine setuprhsall(ndata,mype,init_pass,last_pass) ! not sure why it won't work even in case of .not.luse_obsdiag. if(.not.(l_PBL_pseudo_SurfobsT .or. l_PBL_pseudo_SurfobsQ .or. & - l_PBL_pseudo_SurfobsUV .or. (i_cloud_q_innovation.ge.1)) ) then + l_PBL_pseudo_SurfobsUV .or. (i_cloud_q_innovation>0)) ) then call obsdiags_sort() endif