From f4f0e39ce888cc1dadae5e797e71c4442fe069a7 Mon Sep 17 00:00:00 2001 From: gmao-yzhu Date: Wed, 10 Mar 2021 13:55:10 -0500 Subject: [PATCH] Closes #70. Commit changes for PBL height --- .../GEOSgsi_Coupler/cplr_ensemble.F90 | 20 +- .../GEOSgsi_Coupler/g5_pertmod.F90 | 2 + .../GEOSgsi_Coupler/geos_StateIO.F90 | 60 +++- .../GEOSgsi_Coupler/geos_pertState.F90 | 10 +- .../GEOSgsi_Coupler/geos_pertStateIO.F90 | 73 +++- .../GSI_GridComp/GSI_GridCompMod.F90 | 23 +- .../GSI_GridComp/compute_qvar3d.f90 | 24 +- .../GSI_GridComp/control2model.f90 | 3 + .../GSI_GridComp/control2model_ad.f90 | 3 + .../GSI_GridComp/ensctl2model.f90 | 4 +- .../GSI_GridComp/ensctl2model_ad.f90 | 5 +- .../GSI_GridComp/ensctl2state.f90 | 4 +- .../GSI_GridComp/ensctl2state_ad.f90 | 4 +- .../hybrid_ensemble_isotropic.F90 | 41 +++ GEOSaana_GridComp/GSI_GridComp/inc2guess.f90 | 6 +- GEOSaana_GridComp/GSI_GridComp/intpblh.f90 | 11 +- .../GSI_GridComp/m_berror_stats.f90 | 16 +- GEOSaana_GridComp/GSI_GridComp/m_pblhNode.F90 | 28 +- GEOSaana_GridComp/GSI_GridComp/prewgt.f90 | 12 +- GEOSaana_GridComp/GSI_GridComp/prt_guess.f90 | 15 +- GEOSaana_GridComp/GSI_GridComp/read_obs.F90 | 57 ++- GEOSaana_GridComp/GSI_GridComp/read_pblh.f90 | 337 +++++++++++++++++- GEOSaana_GridComp/GSI_GridComp/setuppblh.f90 | 24 +- GEOSaana_GridComp/GSI_GridComp/stppblh.f90 | 4 +- 24 files changed, 712 insertions(+), 74 deletions(-) diff --git a/GEOSaana_GridComp/GEOSgsi_Coupler/cplr_ensemble.F90 b/GEOSaana_GridComp/GEOSgsi_Coupler/cplr_ensemble.F90 index 832565c7..f6f3bf78 100644 --- a/GEOSaana_GridComp/GEOSgsi_Coupler/cplr_ensemble.F90 +++ b/GEOSaana_GridComp/GEOSgsi_Coupler/cplr_ensemble.F90 @@ -397,6 +397,7 @@ subroutine put_geos_ens(this,grd,member,ntindex,pert,iret) ! 2011-10-01 todling - created for testing purposes ! 2012-05-15 el akkraoui/todling - overload for r4/r8; pass time index ! 2015-10-19 todling - add qi/ql/qr/qs and should work for all cases of CV/SV +! 2020-03-09 yzhu - add pblh ! ! input argument list: ! grd - structure variable containing information about grid @@ -448,9 +449,9 @@ subroutine put_geos_ens(this,grd,member,ntindex,pert,iret) type(gsi_bundle) flds ! Declare fields in file (should come from resource file) - integer(i_kind),parameter:: no2d=2 + integer(i_kind),parameter:: no2d=3 integer(i_kind),dimension(no2d)::iptr2d - character(len=4), parameter :: ovars2d(no2d) = (/ 'ps ', 'z '/) + character(len=4), parameter :: ovars2d(no2d) = (/ 'ps ', 'z ', 'pblh'/) integer(i_kind),parameter:: no3d=10 integer(i_kind),dimension(no3d)::iptr3d character(len=5), parameter :: ovars3d(no3d) = (/ 'u ', 'v ',& @@ -511,6 +512,21 @@ subroutine put_geos_ens(this,grd,member,ntindex,pert,iret) ! if(mype==0) write(6,*) myname,': z ens member not in incoming pert' ! call stop2(999) ! endif + +! PBLH + call gsi_bundlegetpointer (pert,'pblh',ipnt,istatus) + if(istatus==0) then + call gsi_bundlegetpointer (flds,'pblh',optr2d,istatus) + if(iamsingle) then + optr2d=pert%r2(ipnt)%qr4 + else + optr2d=pert%r2(ipnt)%qr8 + endif + else + if(mype==0) write(6,*) myname,': pblh ens member not in incoming pert' + call stop2(999) + endif + ! SF/VP->U/V istatus=0 call gsi_bundlegetpointer (pert,'sf',isf,ier);istatus=ier+istatus diff --git a/GEOSaana_GridComp/GEOSgsi_Coupler/g5_pertmod.F90 b/GEOSaana_GridComp/GEOSgsi_Coupler/g5_pertmod.F90 index ff21edd7..2ba48952 100644 --- a/GEOSaana_GridComp/GEOSgsi_Coupler/g5_pertmod.F90 +++ b/GEOSaana_GridComp/GEOSgsi_Coupler/g5_pertmod.F90 @@ -65,6 +65,7 @@ module g5_pertmod use stepon, only : ak ! GEOS-5 ADM/TLM pressure levels use stepon, only : bk ! GEOS-5 ADM/TLM pressure levels use stepon, only : ts +! use stepon, only : pblht use stepon, only : oro use stepon, only : job use stepon, only : nstep @@ -1494,6 +1495,7 @@ subroutine gsi2pgcm0_ ( nymd, nhms, xx, which, stat, & write(fname,'(4a,i3.3)') trim(job), '.', trim(fnxgsi), '_', mycount endif call putpert ( job, nymd, nhms, xpert, fvpsasdt, nstep, & +! ak, bk, Ts, pblht, oro, ps, fname, vectype=vectype_, forceflip=.true. ) ak, bk, Ts, oro, ps, fname, vectype=vectype_, forceflip=.true. ) if(ierr/=0)then stat = 90 diff --git a/GEOSaana_GridComp/GEOSgsi_Coupler/geos_StateIO.F90 b/GEOSaana_GridComp/GEOSgsi_Coupler/geos_StateIO.F90 index 252f7273..bbaa705b 100644 --- a/GEOSaana_GridComp/GEOSgsi_Coupler/geos_StateIO.F90 +++ b/GEOSaana_GridComp/GEOSgsi_Coupler/geos_StateIO.F90 @@ -456,19 +456,21 @@ subroutine gsi2pgcm1_ ( xx, xpert, sg, GSI_xGrid, which, stat ) ! 15May2010 Todling Update to use GSI_Bundle ! 21Feb2011 Todling Adapt to work with MAPL-SimpleBundle ! 28Oct2013 Todling Rename p3d to prse +! 06Mar2020 Zhu Add pblh ! !EOP !----------------------------------------------------------------------- character(len=*), parameter :: myname_ = myname//'*gsi2pgcm_' - real(r_kind), allocatable, dimension(:,:,:) :: sub_tv,sub_u,sub_v,sub_q,sub_delp + real(r_kind), allocatable, dimension(:,:,:) :: sub_tv,sub_u,sub_v,sub_q,sub_delp,sub_pblh3 real(r_kind), allocatable, dimension(:,:,:) :: sub_oz,sub_cw real(r_kind), allocatable, dimension(:,:,:) :: sub_qi,sub_ql,sub_qr,sub_qs - real(r_kind), allocatable, dimension(:,:) :: sub_ps + real(r_kind), allocatable, dimension(:,:) :: sub_ps,sub_pblh + real(4), allocatable, dimension(:,:,:) :: wktmp integer(i_kind) i,j,k,kk,ijk,ij - integer(i_kind) i_u,i_v,i_t,i_q,i_oz,i_cw,i_prse,i_p,i_dp + integer(i_kind) i_u,i_v,i_t,i_q,i_oz,i_cw,i_prse,i_p,i_dp,i_pblh integer(i_kind) i_qi,i_ql,i_qr,i_qs integer(i_kind) ierr,istatus,rc,status character(len=255) :: whatin @@ -611,6 +613,26 @@ subroutine gsi2pgcm1_ ( xx, xpert, sg, GSI_xGrid, which, stat ) sub_delp(:,:,k) = sub_ps ! just a trick to use gsi2pert_ w/o interface change enddo + call gsi_bundlegetpointer(xx,'pblh', i_pblh , istatus) + if (i_pblh>0) then + allocate(sub_pblh(sg%lat2,sg%lon2), stat=ierr ) + if ( ierr/=0 ) then + stat = 91 + if(mype==ROOT) print*, trim(myname_), ': Alloc(sub_pblh)' + return + end if + do j=1,sg%lon2 + do i=1,sg%lat2 + sub_pblh(i,j) = xx%r2(i_pblh)%q(i,j) + enddo + enddo + allocate(sub_pblh3(sg%lat2,sg%lon2,nsig), stat=ierr ) + sub_pblh3=zero + do k=1,nsig + sub_pblh3(:,:,k) = sub_pblh ! just a trick to use gsi2pert_ w/o interface change + enddo + endif + ! Gather from GSI subdomains/Scatter to GCM ! ----------------------------------------- i_p = MAPL_SimpleBundleGetIndex ( xpert, 'ps' , 2, rc=status ) @@ -648,6 +670,13 @@ subroutine gsi2pgcm1_ ( xx, xpert, sg, GSI_xGrid, which, stat ) i_qs= MAPL_SimpleBundleGetIndex ( xpert, 'qstot',3, rc=status ) call gsi2pert_ ( sub_qs, xpert%r3(i_qs)%qr4, ierr ) endif + if (i_pblh>0) then + i_pblh= MAPL_SimpleBundleGetIndex ( xpert, 'PBLH',2, rc=status ) + allocate(wktmp(sg%lat2,sg%lon2,nsig), stat=ierr ) + wktmp = zero + call gsi2pert_ ( sub_pblh3, wktmp, ierr ) + xpert%r2(i_pblh)%qr4 = wktmp(:,:,1) + endif ! Build surface pressure perturbation - output purposes ! ----------------------------------------------------- @@ -670,6 +699,9 @@ subroutine gsi2pgcm1_ ( xx, xpert, sg, GSI_xGrid, which, stat ) if (allocated(sub_ql)) deallocate(sub_ql) if (allocated(sub_qr)) deallocate(sub_qr) if (allocated(sub_qs)) deallocate(sub_qs) + if (allocated(sub_pblh)) deallocate(sub_pblh) + if (allocated(sub_pblh3)) deallocate(sub_pblh3) + if (allocated(wktmp)) deallocate(wktmp) CONTAINS @@ -1161,6 +1193,7 @@ subroutine gcm2gsi1_ ( xpert, xx, sg, GSI_xGrid, stat ) ! 20Feb2011 Todling Adapt to work with MAPL-SimpleBundle ! 30Nov2014 Todling Add some variable exceptions (qi/ql/ph/ts) ! 06May2020 Todling Shuffled pointer check and alloc for increased flexibility +! 06Mar2020 Zhu Add pblh ! !EOP !----------------------------------------------------------------------- @@ -1171,11 +1204,11 @@ subroutine gcm2gsi1_ ( xpert, xx, sg, GSI_xGrid, stat ) real(r_kind), allocatable, dimension(:,:,:) :: sub_u,sub_v,sub_q,sub_tv real(r_kind), allocatable, dimension(:,:,:) :: sub_oz,sub_ci,sub_cl real(r_kind), allocatable, dimension(:,:,:) :: sub_cr,sub_cs - real(r_kind), allocatable, dimension(:,:) :: sub_ps,sub_ph,sub_ts + real(r_kind), allocatable, dimension(:,:) :: sub_ps,sub_ph,sub_ts,sub_pblh real(r_kind), allocatable, dimension(:,:,:) :: sub_3d integer(i_kind) i,j,k,ij,ijk,id,ng_d,ng_s,ni,nj - integer(i_kind) i_u,i_v,i_t,i_q,i_oz,i_cw,i_ps,i_ts,i_tsen,i_ph + integer(i_kind) i_u,i_v,i_t,i_q,i_oz,i_cw,i_ps,i_ts,i_tsen,i_ph,i_pblh integer(i_kind) i_ci,i_cl,i_cr,i_cs integer(i_kind) j_cr,j_cs integer(i_kind) ierr,istatus,status @@ -1200,6 +1233,7 @@ subroutine gcm2gsi1_ ( xpert, xx, sg, GSI_xGrid, stat ) i_ps = MAPL_SimpleBundleGetIndex ( xpert, 'ps' , 2, rc=status ) i_ts = MAPL_SimpleBundleGetIndex ( xpert, 'ts' , 2, rc=status ) i_ph = MAPL_SimpleBundleGetIndex ( xpert, 'phis' , 2, rc=status ) + i_pblh= MAPL_SimpleBundleGetIndex (xpert, 'PBLH' , 2, rc=status ) if (i_u>0.and.i_v>0) then allocate(sub_u(sg%lat2,sg%lon2,nsig), stat=ierr ) call pert2gsi_ ( xpert%r3(i_u)%qr4 , sub_u , ierr ) @@ -1238,6 +1272,10 @@ subroutine gcm2gsi1_ ( xpert, xx, sg, GSI_xGrid, stat ) allocate(sub_ts(sg%lat2,sg%lon2)) call pert2gsi2d_ ( xpert%r2(i_ts)%qr4, sub_ts, ierr ) endif + if (i_pblh>0) then + allocate(sub_pblh(sg%lat2,sg%lon2)) + call pert2gsi2d_ ( xpert%r2(i_pblh)%qr4, sub_pblh, ierr ) + endif if ( ierr/=0 ) then stat = 99 if(mype==ROOT) print*, trim(myname_), ': unfinished convertion ...' @@ -1428,6 +1466,18 @@ subroutine gcm2gsi1_ ( xpert, xx, sg, GSI_xGrid, stat ) deallocate(sub_ts) endif + if (allocated(sub_pblh)) then + call gsi_bundlegetpointer(xx,'pblh', i_pblh, ierr) + if (i_pblh>0) then + do j=1,sg%lon2 + do i=1,sg%lat2 + xx%r2(i_pblh)%q(i,j)= sub_pblh(i,j) + enddo + enddo + endif + deallocate(sub_pblh) + endif + ! check on essential vars if ( ierr/=0 ) then stat = 99 diff --git a/GEOSaana_GridComp/GEOSgsi_Coupler/geos_pertState.F90 b/GEOSaana_GridComp/GEOSgsi_Coupler/geos_pertState.F90 index 6f84ab5e..8db163e1 100644 --- a/GEOSaana_GridComp/GEOSgsi_Coupler/geos_pertState.F90 +++ b/GEOSaana_GridComp/GEOSgsi_Coupler/geos_pertState.F90 @@ -103,6 +103,7 @@ module geos_pertState ! and it has be verified against a GSI built-in ! adjoint-test in GSI evaljgrad(). ! 28Oct2013 Todling Rename p3d to prse +! 06Mar2020 Zhu Add pblh ! !EOP !------------------------------------------------------------------------- @@ -152,13 +153,14 @@ module geos_pertState character(len=*),parameter:: oz_GSIpert='oz' ! = O3, in GG=gram/gram character(len=*),parameter:: cw_GSIpert='cw' ! = QL + QI character(len=*),parameter:: sst_GSIpert='sst' ! = ? + character(len=*),parameter:: pblh_GSIpert='pblh' ! = ZPBL character(len=*),dimension(8),parameter:: var3dList_GSIpert = & (/ u_GSIpert, v_GSIpert,prse_GSIpert, tv_GSIpert, & tsen_GSIpert, q_GSIpert, oz_GSIpert, cw_GSIpert /) - character(len=*),dimension(2),parameter:: var2dList_GSIpert = & - (/ps_GSIpert, sst_GSIpert /) + character(len=*),dimension(3),parameter:: var2dList_GSIpert = & + (/ps_GSIpert, sst_GSIpert, pblh_GSIpert /) ! These are expected from a PGCM state character(len=*),parameter:: u_GCMpert='U' ! = u @@ -169,11 +171,15 @@ module geos_pertState character(len=*),parameter:: ql_GCMpert='QL' ! = cw - qi character(len=*),parameter:: qi_GCMpert='QI' ! = qi (=0) character(len=*),parameter:: o3_GCMpert='O3' ! = oz, in PPMV + character(len=*),parameter:: pblh_GCMpert='PBLH' ! = zpbl character(len=*),dimension(8),parameter:: var3dList_GCMpert = & (/ u_GCMpert, v_GCMpert, tv_GCMpert, dp_GCMpert, & qv_GCMpert, ql_GCMpert, qi_GCMpert, o3_GCMpert /) + character(len=*),dimension(1),parameter:: var2dList_GCMpert = & + (/ pblh_GCMpert /) + integer(i_kind),parameter:: ROOT =0 real(r_kind), parameter :: kPa_per_Pa = 0.001_r_kind real(r_kind), parameter :: Pa_per_kPa = 1000._r_kind diff --git a/GEOSaana_GridComp/GEOSgsi_Coupler/geos_pertStateIO.F90 b/GEOSaana_GridComp/GEOSgsi_Coupler/geos_pertStateIO.F90 index 5fc8d364..ed915f4d 100644 --- a/GEOSaana_GridComp/GEOSgsi_Coupler/geos_pertStateIO.F90 +++ b/GEOSaana_GridComp/GEOSgsi_Coupler/geos_pertStateIO.F90 @@ -337,6 +337,7 @@ subroutine set_gsi2pgcm0_ ( nymd, nhms, stat ) ! !REVISION HISTORY: ! ! 03Mar2020 Todling Broken up from original gsi2pgcm0_ +! 07Mar2020 Zhu Add pblh ! !EOP !----------------------------------------------------------------------- @@ -351,7 +352,7 @@ subroutine set_gsi2pgcm0_ ( nymd, nhms, stat ) character(len=255) fname,bkgfname,tmpl character(len=ESMF_MAXSTR) :: etmpl - character(len=*), parameter :: only_vars='ps,ts,delp,u,v,tv,sphu,ozone,qitot,qltot,qrtot,qstot' + character(len=*), parameter :: only_vars='ps,ts,PBLH,delp,u,v,tv,sphu,ozone,qitot,qltot,qrtot,qstot' integer(i_kind) thistime(6) integer(i_kind) rc,status,ierr integer(i_kind) idim,jdim,kdim,i,j,k,i_dp,i_ps @@ -477,7 +478,7 @@ subroutine gsi2pgcm0_ ( nymd, nhms, xx, which, stat, & character(len=255) fname,bkgfname,tmpl character(len=ESMF_MAXSTR) :: etmpl - character(len=*), parameter :: only_vars='ps,ts,delp,u,v,tv,sphu,ozone,qitot,qltot,qrtot,qstot' + character(len=*), parameter :: only_vars='ps,ts,PBLH,delp,u,v,tv,sphu,ozone,qitot,qltot,qrtot,qstot' integer(i_kind) thistime(6) integer(i_kind) rc,status,ierr integer(i_kind) idim,jdim,kdim,i,j,k,i_dp,i_ps @@ -642,7 +643,7 @@ subroutine final_gsi2pgcm0_ ( stat ) character(len=255) fname,bkgfname,tmpl character(len=ESMF_MAXSTR) :: etmpl - character(len=*), parameter :: only_vars='ps,ts,delp,u,v,tv,sphu,ozone,qitot,qltot,qrtot,qstot' + character(len=*), parameter :: only_vars='ps,ts,PBLH,delp,u,v,tv,sphu,ozone,qitot,qltot,qrtot,qstot' integer(i_kind) thistime(6) integer(i_kind) rc,status,ierr integer(i_kind) idim,jdim,kdim,i,j,k,i_dp,i_ps @@ -760,11 +761,11 @@ subroutine gsi2pgcm1_ ( xx, xpert, which, stat ) real(r_kind), allocatable, dimension(:,:,:) :: sub_oz,sub_cw real(r_kind), allocatable, dimension(:,:,:) :: sub_qi,sub_ql real(r_kind), allocatable, dimension(:,:,:) :: sub_qr,sub_qs - real(r_kind), allocatable, dimension(:,:) :: sub_ps,sub_ts + real(r_kind), allocatable, dimension(:,:) :: sub_ps,sub_ts,sub_pblh real(r_single), allocatable, dimension(:,:,:) :: aux3d integer(i_kind) i,j,k,ijk,ij - integer(i_kind) i_u,i_v,i_t,i_q,i_oz,i_cw,i_prse,i_p,i_tsen,i_ts + integer(i_kind) i_u,i_v,i_t,i_q,i_oz,i_cw,i_prse,i_p,i_tsen,i_ts,i_pblh integer(i_kind) i_qi,i_ql,i_qr,i_qs integer(i_kind) ierr,ierr_adm,istatus,rc,status character(len=255) :: whatin @@ -935,6 +936,20 @@ subroutine gsi2pgcm1_ ( xx, xpert, which, stat ) enddo enddo endif + call gsi_bundlegetpointer(xx,'pblh', i_pblh, istatus) + if (i_pblh>0) then + allocate (sub_pblh(lat2,lon2), stat=ierr ) + if ( ierr/=0 ) then + stat = 91 + if(mype==ROOT) print*, trim(myname_), ': Alloc(sub_pblh)' + return + end if + do j=1,lon2 + do i=1,lat2 + sub_pblh(i,j) = xx%r2(i_pblh)%q(i,j) + enddo + enddo + endif ! Calculate perturbation delp ! --------------------------- @@ -998,6 +1013,16 @@ subroutine gsi2pgcm1_ ( xx, xpert, which, stat ) xpert%r2(i_ts)%qr4=aux3d(:,:,nsig) ! level nsig will contain result give vertical swap deallocate(aux3d) + if (i_pblh>0) then + i_pblh = MAPL_SimpleBundleGetIndex ( xpert, 'PBLH' , 2, rc=ierr ) + sub_q=zero + sub_q(:,:,1)=sub_pblh ! copy to level 1 + allocate(aux3d(size(xpert%r3(i_q)%qr4,1),size(xpert%r3(i_q)%qr4,2),size(xpert%r3(i_q)%qr4,3))) + call gsi2pert_ ( sub_q, aux3d, ierr ) + xpert%r2(i_pblh)%qr4=aux3d(:,:,nsig) ! level nsig will contain result give vertical swap + deallocate(aux3d) + end if + deallocate (sub_tv, sub_u, sub_v, sub_q, sub_delp, sub_ps, sub_ts, stat=ierr ) if ( ierr/=0 ) then stat = 99 @@ -1052,7 +1077,14 @@ subroutine gsi2pgcm1_ ( xx, xpert, which, stat ) return end if endif - + if ( allocated(sub_pblh) ) then + deallocate (sub_pblh, stat=ierr ) + if ( ierr/=0 ) then + stat = 99 + if(mype==ROOT) print*, trim(myname_), ': Dealloc(sub_pblh)' + return + end if + endif @@ -1323,7 +1355,7 @@ subroutine pgcm2gsi0_ ( xx, which, stat, & character(len=*), parameter :: Iam = myname_ character(len=*), parameter :: fname_def ='fsens.eta.nc4' ! full forecast sensitivity from ADM integration character(len=*), parameter :: sens_vars='ts,delp,u,v,tv,sphu,ozone' - character(len=*), parameter :: xinc_vars='ts,ps,u,v,tv,sphu,ozone,qitot,qltot,qrtot,qstot' + character(len=*), parameter :: xinc_vars='ts,ps,PBLH,u,v,tv,sphu,ozone,qitot,qltot,qrtot,qstot' character(len=256) :: fname, ffname(1) character(len=256) :: my_vars character(len=30) :: GSIGRIDNAME @@ -1524,6 +1556,7 @@ subroutine pgcm2gsi1_ ( xpert, xx, which, stat, jgradf ) ! 28Oct2013 Todling Rename p3d to prse ! 12Jul2016 Todling Slightly revisit handle for tracers ! 08Oct2016 M.Kim Allow for input of qr/qs +! 06Mar2020 Zhu Add pblh ! !EOP !----------------------------------------------------------------------- @@ -1533,12 +1566,12 @@ subroutine pgcm2gsi1_ ( xpert, xx, which, stat, jgradf ) real(r_kind), allocatable, dimension(:,:,:) :: sub_u,sub_v,sub_delp,sub_q,sub_tv real(r_kind), allocatable, dimension(:,:,:) :: sub_oz,sub_cw real(r_kind), allocatable, dimension(:,:,:) :: sub_ql,sub_qi,sub_qr,sub_qs - real(r_kind), allocatable, dimension(:,:) :: sub_ps,sub_ts + real(r_kind), allocatable, dimension(:,:) :: sub_ps,sub_ts,sub_pblh real(r_kind), allocatable, dimension(:,:,:) :: aux3d character(len=256) :: failvars integer(i_kind) i,j,k,ij,ijk,id,ng_d,ng_s - integer(i_kind) i_u,i_v,i_t,i_q,i_oz,i_cw,i_prse,i_p,i_tsen,i_ts,i_ql,i_qi,i_qr,i_qs + integer(i_kind) i_u,i_v,i_t,i_q,i_oz,i_cw,i_prse,i_p,i_tsen,i_ts,i_ql,i_qi,i_qr,i_qs,i_pblh integer(i_kind) slat2,slon2 integer(i_kind) ierr,istatus,status logical scaleit @@ -1654,6 +1687,11 @@ subroutine pgcm2gsi1_ ( xpert, xx, which, stat, jgradf ) allocate (sub_ts(slat2,slon2), stat=istatus ); ierr=ierr+istatus call pert2gsi2d_ ( xpert%r2(i_ts)%qr4, sub_ts, ierr ) endif + i_pblh= MAPL_SimpleBundleGetIndex ( xpert, 'PBLH' , 2, rc=status ) + if (i_pblh>0) then + allocate (sub_pblh(lat2,lon2), stat=istatus ); ierr=ierr+istatus + call pert2gsi2d_ ( xpert%r2(i_pblh)%qr4, sub_pblh, ierr ) + end if if (which=='inc') then i_p= MAPL_SimpleBundleGetIndex ( xpert, 'ps', 2, rc=status ) @@ -1936,6 +1974,23 @@ subroutine pgcm2gsi1_ ( xpert, xx, which, stat, jgradf ) return end if endif + + if (allocated(sub_pblh) ) then + call gsi_bundlegetpointer(xx,'pblh',i_pblh,istatus) + if (i_pblh>0) then + do j=1,lon2 + do i=1,lat2 + xx%r2(i_pblh)%q(i,j) = sub_pblh(i,j) + enddo + enddo + end if + deallocate (sub_pblh, stat=ierr ) + if ( ierr/=0 ) then + stat = 99 + if(mype==ROOT) print*, trim(myname_), ': Dealloc(sub_pblh)' + return + end if + end if CONTAINS diff --git a/GEOSaana_GridComp/GSI_GridComp/GSI_GridCompMod.F90 b/GEOSaana_GridComp/GSI_GridComp/GSI_GridCompMod.F90 index eb33db7c..060e74bf 100644 --- a/GEOSaana_GridComp/GSI_GridComp/GSI_GridCompMod.F90 +++ b/GEOSaana_GridComp/GSI_GridComp/GSI_GridCompMod.F90 @@ -235,6 +235,7 @@ module GSI_GridCompMod ! - Remove reference to vars not used here ! 18Sep2012 Akella - Add SwapJI2i_ to GSI_GridCompSwapJI_ for use in writing out TSKIN when doing NST analysis ! 19Oct2013 Todling - metguess now holds background +! 05Mar2021 Zhu - Add pblh ! !EOP !------------------------------------------------------------------------- @@ -1186,6 +1187,7 @@ subroutine Run ( gc, import, export, clock, rc ) ! import state upper air pointers real(r_single),dimension(:,: ), pointer :: hsp ! terrain real(r_single),dimension(:,: ), pointer :: psp ! surf. pressure + real(r_single),dimension(:,: ), pointer :: pblhp ! pbl height real(r_single),dimension(:,:,:), pointer :: up ! u wind real(r_single),dimension(:,:,:), pointer :: vp ! v wind real(r_single),dimension(:,:,:), pointer :: tp ! virtual temp. @@ -1252,6 +1254,7 @@ subroutine Run ( gc, import, export, clock, rc ) real(r_single),dimension(:,: ), pointer :: dts ! skin/surface temperature real(r_single),dimension(:,: ), pointer :: dhs ! terrain real(r_single),dimension(:,: ), pointer :: dps ! surf pressure + real(r_single),dimension(:,: ), pointer :: dpblh ! pbl height real(r_single),dimension(:,:,:), pointer :: ddp ! del pressure real(r_single),dimension(:,:,:), pointer :: du ! u wind real(r_single),dimension(:,:,:), pointer :: dv ! v wind @@ -1563,6 +1566,9 @@ subroutine GSI_GridCompGetPointers_() case ( 'PS' ) call ESMFL_StateGetPointerToData(import, psp,trim(cvar), rc=STATUS) VERIFY_(STATUS) + case ( 'PBLH' ) + call ESMFL_StateGetPointerToData(import, pblhp,trim(cvar), rc=STATUS) + VERIFY_(STATUS) case ( 'PHIS' ) call ESMFL_StateGetPointerToData(import, hsp,trim(cvar), rc=STATUS) VERIFY_(STATUS) @@ -1724,6 +1730,9 @@ subroutine GSI_GridCompGetPointers_() case ('ps') call ESMFL_StateGetPointerToData(export, dps, trim(cvar), alloc=.true., rc=STATUS) VERIFY_(STATUS) + case ('PBLH') + call ESMFL_StateGetPointerToData(export, dpblh, trim(cvar), alloc=.true., rc=STATUS) + VERIFY_(STATUS) case ('u') call ESMFL_StateGetPointerToData(export, du, trim(cvar), alloc=.true., rc=STATUS) VERIFY_(STATUS) @@ -2015,6 +2024,8 @@ subroutine GSI_GridCompCopyImportDyn2Internal_(lit) call GSI_GridCompSwapIJ_(tdelp,GSI_MetGuess_bundle(it)%r2(ipnt)%q) case ('tref') call GSI_GridCompSwapIJ_(trefp ,GSI_MetGuess_bundle(it)%r2(ipnt)%q) + case ('pblh') + call GSI_GridCompSwapIJ_(pblhp ,GSI_MetGuess_bundle(it)%r2(ipnt)%q) end select #ifdef SFCverbose call guess_grids_stats(cvar, GSI_MetGuess_Bundle(it)%r2(ipnt)%q, mype) @@ -2848,6 +2859,8 @@ subroutine GSI_GridCompCopyInternal2Export_(lit) select case (cvar) case('ps') call GSI_GridCompSwapJI_(dps,GSI_MetGuess_Bundle(it)%r2(ipnt)%q) + case('pblh') + call GSI_GridCompSwapJI_(dpblh,GSI_MetGuess_Bundle(it)%r2(ipnt)%q) case('z') call GSI_GridCompSwapJI_(dhs,GSI_MetGuess_Bundle(it)%r2(ipnt)%q) end select @@ -3672,10 +3685,11 @@ subroutine GSI_GridCompSetupSpecs (GC, opthw, rc) ! Declare import 2d-fields ! ------------------------ - integer, parameter :: nin2d=14 + integer, parameter :: nin2d=15 character(len=16), parameter :: insname2d(nin2d) = (/ & 'phis ', & 'ps ', & + 'PBLH ', & 'ts ', & ! 'NCEP_VEGFRAC ', & ! 'NCEP_VEGTYPE ', & @@ -3694,6 +3708,7 @@ subroutine GSI_GridCompSetupSpecs (GC, opthw, rc) character(len=32), parameter :: inlname2d(nin2d) = (/ & 'geopotential height ', & 'surface pressure ', & + 'planetary boundary layer height ', & 'skin temperature ', & ! 'NCEP(CRTM-like) veg fraction ', & ! 'NCEP(CRTM-like) veg type ', & @@ -3712,6 +3727,7 @@ subroutine GSI_GridCompSetupSpecs (GC, opthw, rc) character(len=16), parameter :: inunits2d(nin2d) = (/ & 'm**2/s**2 ', & 'Pa ', & + 'm ', & 'K ', & ! '1 ', & ! '1 ', & @@ -3876,10 +3892,11 @@ subroutine GSI_GridCompSetupSpecs (GC, opthw, rc) ! Declare export 2d-fields ! ------------------------ - integer, parameter :: nex2d=8 + integer, parameter :: nex2d=9 character(len=16), parameter :: exsname2d(nex2d) = (/ & 'phis ', & 'ps ', & + 'PBLH ', & 'ts ', & 'frland ', & 'frlandice ', & @@ -3889,6 +3906,7 @@ subroutine GSI_GridCompSetupSpecs (GC, opthw, rc) character(len=32), parameter :: exlname2d(nex2d) = (/ & 'geopotential height ', & 'surface pressure inc ', & + 'pbl height inc ', & 'skin temperature inc ', & 'fraction of land ', & 'fraction of land ice ', & @@ -3898,6 +3916,7 @@ subroutine GSI_GridCompSetupSpecs (GC, opthw, rc) character(len=16), parameter :: exunits2d(nex2d) = (/ & 'm**2/s**2 ', & 'Pa ', & + 'm ', & 'K ', & '1 ', & '1 ', & diff --git a/GEOSaana_GridComp/GSI_GridComp/compute_qvar3d.f90 b/GEOSaana_GridComp/GSI_GridComp/compute_qvar3d.f90 index 54aa7721..3665de2b 100644 --- a/GEOSaana_GridComp/GSI_GridComp/compute_qvar3d.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/compute_qvar3d.f90 @@ -24,6 +24,7 @@ subroutine compute_qvar3d ! 2015-09-10 zhu - use centralized cloud_names_fwd and n_clouds_fwd in the assignments of cloud ! variances (either cw or individual hydrometerors) for all-sky radiance assimilation ! - remove cwoption1 +! 2021-03-03 zhu - add pblh ! ! input argument list: ! @@ -36,10 +37,10 @@ subroutine compute_qvar3d ! !$$$ use kinds, only: r_kind,i_kind,r_single - use berror, only: dssv + use berror, only: dssv,dssvs use jfunc, only: varq,qoption,varcw,cwoption,clip_supersaturation use derivsmod, only: qsatg,qgues - use control_vectors, only: cvars3d + use control_vectors, only: cvars3d,cvars2d use gridmod, only: lat2,lon2,nsig use constants, only: zero,one,fv,r100,qmin use guess_grids, only: fact_tv,ntguessig,nfldsig,ges_tsen,ges_prsl,ges_qsat @@ -54,13 +55,13 @@ subroutine compute_qvar3d ! Declare local variables logical ice - integer(i_kind) :: i,j,k,it,n,np,iderivative,nrf3_q,nrf3_cw + integer(i_kind) :: i,j,k,it,n,np,iderivative,nrf3_q,nrf3_cw,nrf2_pblh integer(i_kind) :: ic,nrf3_var real(r_kind) d,dn1,dn2 real(r_kind),allocatable,dimension(:,:,:):: rhgues integer(i_kind):: istatus,ier,ier6 - real(r_kind):: cwtmp + real(r_kind):: cwtmp,htmp real(r_kind),pointer,dimension(:,:,:):: ges_var=>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_ql=>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_qi=>NULL() @@ -69,11 +70,13 @@ subroutine compute_qvar3d real(r_kind),pointer,dimension(:,:,:):: ges_qg=>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_qh=>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_q =>NULL() + real(r_kind),pointer,dimension(:,:):: ges_pblh =>NULL() integer(i_kind):: maxvarq1 nrf3_q=getindex(cvars3d,'q') nrf3_cw=getindex(cvars3d,'cw') + nrf2_pblh=getindex(cvars2d,'pblh') ! Calculate qsat independently of presence of q in guess iderivative = 0 @@ -153,6 +156,19 @@ subroutine compute_qvar3d deallocate(rhgues) + if (nrf2_pblh>0) then + call gsi_bundlegetpointer (gsi_metguess_bundle(ntguessig),'pblh',ges_pblh,istatus) + if (istatus/=0) return + do j = 1,lon2 + do i = 1,lat2 + htmp=ges_pblh(i,j) + if (ges_pblh(i,j)0) then diff --git a/GEOSaana_GridComp/GSI_GridComp/control2model.f90 b/GEOSaana_GridComp/GSI_GridComp/control2model.f90 index ec628afe..3705ab4e 100644 --- a/GEOSaana_GridComp/GSI_GridComp/control2model.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/control2model.f90 @@ -93,6 +93,7 @@ subroutine control2model(xhat,sval,bval) ! motley vars (if any) real(r_kind),pointer,dimension(:,:) :: sv_ps=>NULL() real(r_kind),pointer,dimension(:,:) :: sv_sst=>NULL() +real(r_kind),pointer,dimension(:,:) :: sv_pblh=>NULL() real(r_kind),pointer,dimension(:,:,:) :: sv_u=>NULL() real(r_kind),pointer,dimension(:,:,:) :: sv_v=>NULL() real(r_kind),pointer,dimension(:,:,:) :: sv_prse=>NULL() @@ -185,6 +186,7 @@ subroutine control2model(xhat,sval,bval) call gsi_bundlegetpointer (sval(jj),'q' ,sv_q , istatus) call gsi_bundlegetpointer (sval(jj),'oz' ,sv_oz , istatus) call gsi_bundlegetpointer (sval(jj),'sst' ,sv_sst, istatus) + call gsi_bundlegetpointer (sval(jj),'pblh' ,sv_pblh, istatus) ! Copy variables from CV to SV call gsi_bundlegetvar ( wbundle, 'sf' , workst, istatus ) @@ -194,6 +196,7 @@ subroutine control2model(xhat,sval,bval) call gsi_bundlegetvar ( wbundle, 'oz' , sv_oz, istatus ) call gsi_bundlegetvar ( wbundle, 'ps' , sv_ps, istatus ) call gsi_bundlegetvar ( wbundle, 'sst', sv_sst, istatus ) + call gsi_bundlegetvar ( wbundle, 'pblh', sv_pblh, istatus ) ! Same one-to-one map for chemistry-vars; take care of them together do ic=1,ngases diff --git a/GEOSaana_GridComp/GSI_GridComp/control2model_ad.f90 b/GEOSaana_GridComp/GSI_GridComp/control2model_ad.f90 index 639bf87f..31476193 100644 --- a/GEOSaana_GridComp/GSI_GridComp/control2model_ad.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/control2model_ad.f90 @@ -87,6 +87,7 @@ subroutine control2model_ad(rval,bval,grad) real(r_kind),pointer,dimension(:,:) :: rv_ps=>NULL() real(r_kind),pointer,dimension(:,:) :: rv_sst=>NULL() +real(r_kind),pointer,dimension(:,:) :: rv_pblh=>NULL() real(r_kind),pointer,dimension(:,:,:) :: rv_u=>NULL() real(r_kind),pointer,dimension(:,:,:) :: rv_v=>NULL() real(r_kind),pointer,dimension(:,:,:) :: rv_prse=>NULL() @@ -165,6 +166,7 @@ subroutine control2model_ad(rval,bval,grad) call gsi_bundlegetpointer (rval(jj),'q' ,rv_q , istatus) call gsi_bundlegetpointer (rval(jj),'oz' ,rv_oz , istatus) call gsi_bundlegetpointer (rval(jj),'sst' ,rv_sst, istatus) + call gsi_bundlegetpointer (rval(jj),'pblh',rv_pblh, istatus) ! Convert RHS calculations for u,v to st/vp for application of ! background error @@ -204,6 +206,7 @@ subroutine control2model_ad(rval,bval,grad) call gsi_bundleputvar ( wbundle, 'ps', rv_ps, istatus ) call gsi_bundleputvar ( wbundle, 'oz', rv_oz, istatus ) call gsi_bundleputvar ( wbundle, 'sst',rv_sst, istatus ) + call gsi_bundleputvar ( wbundle, 'pblh',rv_pblh, istatus ) if (nclouds>0) then if (cw_to_hydro_ad) then diff --git a/GEOSaana_GridComp/GSI_GridComp/ensctl2model.f90 b/GEOSaana_GridComp/GSI_GridComp/ensctl2model.f90 index 525bf33a..05b69466 100644 --- a/GEOSaana_GridComp/GSI_GridComp/ensctl2model.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/ensctl2model.f90 @@ -71,7 +71,7 @@ subroutine ensctl2model(xhat,mval,eval) character(len=4), parameter :: mysvars(nsvars) = (/ & ! vars from ST needed here 'u ', 'v ', 'prse', 'q ', 'tsen' /) logical :: ls_u,ls_v,ls_prse,ls_q,ls_tsen -real(r_kind),pointer,dimension(:,:) :: sv_ps,sv_sst +real(r_kind),pointer,dimension(:,:) :: sv_ps,sv_sst,sv_pblh real(r_kind),pointer,dimension(:,:,:) :: sv_u,sv_v,sv_prse,sv_q,sv_tsen,sv_tv,sv_oz real(r_kind),pointer,dimension(:,:,:) :: sv_rank3 @@ -161,6 +161,7 @@ subroutine ensctl2model(xhat,mval,eval) call gsi_bundlegetpointer (eval(jj),'q' ,sv_q , istatus) call gsi_bundlegetpointer (eval(jj),'oz' ,sv_oz , istatus) call gsi_bundlegetpointer (eval(jj),'sst' ,sv_sst, istatus) + call gsi_bundlegetpointer (eval(jj),'pblh' ,sv_pblh, istatus) call sqrt_beta_s_mult(wbundle_c) if(dual_res) then @@ -197,6 +198,7 @@ subroutine ensctl2model(xhat,mval,eval) call gsi_bundlegetvar ( wbundle_c, 'oz' , sv_oz, istatus ) call gsi_bundlegetvar ( wbundle_c, 'ps' , sv_ps, istatus ) call gsi_bundlegetvar ( wbundle_c, 'sst', sv_sst, istatus ) + call gsi_bundlegetvar ( wbundle_c, 'pblh', sv_pblh, istatus ) ! Since cloud-vars map one-to-one, take care of them together do ic=1,nclouds diff --git a/GEOSaana_GridComp/GSI_GridComp/ensctl2model_ad.f90 b/GEOSaana_GridComp/GSI_GridComp/ensctl2model_ad.f90 index 8c491b6d..dfc70d1c 100644 --- a/GEOSaana_GridComp/GSI_GridComp/ensctl2model_ad.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/ensctl2model_ad.f90 @@ -71,7 +71,7 @@ subroutine ensctl2model_ad(eval,mval,grad) character(len=4), parameter :: mysvars(nsvars) = (/ & ! vars from ST needed here 'u ', 'v ', 'prse', 'q ', 'tsen' /) logical :: ls_u,ls_v,ls_prse,ls_q,ls_tsen -real(r_kind),pointer,dimension(:,:) :: rv_ps,rv_sst +real(r_kind),pointer,dimension(:,:) :: rv_ps,rv_sst,rv_pblh real(r_kind),pointer,dimension(:,:,:) :: rv_u,rv_v,rv_prse,rv_q,rv_tsen,rv_tv,rv_oz real(r_kind),pointer,dimension(:,:,:) :: rv_rank3 @@ -156,6 +156,7 @@ subroutine ensctl2model_ad(eval,mval,grad) call gsi_bundlegetpointer (eval(jj),'q' ,rv_q , istatus) call gsi_bundlegetpointer (eval(jj),'oz' ,rv_oz , istatus) call gsi_bundlegetpointer (eval(jj),'sst' ,rv_sst, istatus) + call gsi_bundlegetpointer (eval(jj),'pblh' ,rv_pblh, istatus) ! Adjoint of consistency for sensible temperature, calculate sensible temperature if(do_tv_to_tsen_ad) call tv_to_tsen_ad(rv_tv,rv_q,rv_tsen) @@ -178,6 +179,8 @@ subroutine ensctl2model_ad(eval,mval,grad) call gsi_bundleputvar ( wbundle_c, 'ps', rv_ps, istatus ) call gsi_bundleputvar ( wbundle_c, 'oz', rv_oz, istatus ) call gsi_bundleputvar ( wbundle_c, 'sst', rv_sst, istatus ) + call gsi_bundleputvar ( wbundle_c, 'pblh', rv_pblh, istatus ) + if (istatus/=0) print*, 'ensctl2model_ad not find pblh' ! Since cloud-vars map one-to-one, take care of them together do ic=1,nclouds diff --git a/GEOSaana_GridComp/GSI_GridComp/ensctl2state.f90 b/GEOSaana_GridComp/GSI_GridComp/ensctl2state.f90 index f3712320..82bb9fad 100644 --- a/GEOSaana_GridComp/GSI_GridComp/ensctl2state.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/ensctl2state.f90 @@ -77,7 +77,7 @@ subroutine ensctl2state(xhat,mval,eval) logical :: ls_u,ls_v,ls_prse,ls_q,ls_tsen,ls_ql,ls_qi logical :: ls_qr,ls_qs,ls_qg,ls_qh logical :: ls_w,ls_dw -real(r_kind),pointer,dimension(:,:) :: sv_ps,sv_sst +real(r_kind),pointer,dimension(:,:) :: sv_ps,sv_sst,sv_pblh real(r_kind),pointer,dimension(:,:,:) :: sv_u,sv_v,sv_prse,sv_q,sv_tsen,sv_tv,sv_oz real(r_kind),pointer,dimension(:,:,:) :: sv_rank3,sv_w,sv_dw @@ -231,6 +231,7 @@ subroutine ensctl2state(xhat,mval,eval) ! Get pointers to required state variables call gsi_bundlegetpointer (eval(jj),'oz' ,sv_oz , istatus) call gsi_bundlegetpointer (eval(jj),'sst' ,sv_sst, istatus) + call gsi_bundlegetpointer (eval(jj),'pblh',sv_pblh,istatus) if(ls_w)then call gsi_bundlegetpointer (eval(jj),'w' ,sv_w, istatus) if(ls_dw.and.nems_nmmb_regional)then @@ -240,6 +241,7 @@ subroutine ensctl2state(xhat,mval,eval) ! Copy variables call gsi_bundlegetvar ( wbundle_c, 'oz' , sv_oz, istatus ) call gsi_bundlegetvar ( wbundle_c, 'sst', sv_sst, istatus ) + call gsi_bundlegetvar ( wbundle_c, 'pblh', sv_pblh, istatus ) if(lc_w)then call gsi_bundlegetvar ( wbundle_c, 'w' , sv_w, istatus ) if(lc_dw.and.nems_nmmb_regional)then diff --git a/GEOSaana_GridComp/GSI_GridComp/ensctl2state_ad.f90 b/GEOSaana_GridComp/GSI_GridComp/ensctl2state_ad.f90 index 886fe3a3..001f73b3 100644 --- a/GEOSaana_GridComp/GSI_GridComp/ensctl2state_ad.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/ensctl2state_ad.f90 @@ -76,7 +76,7 @@ subroutine ensctl2state_ad(eval,mval,grad) logical :: ls_u,ls_v,ls_prse,ls_q,ls_tsen,ls_ql,ls_qi logical :: ls_qr,ls_qs,ls_qg,ls_qh logical :: ls_w,ls_dw -real(r_kind),pointer,dimension(:,:) :: rv_ps,rv_sst +real(r_kind),pointer,dimension(:,:) :: rv_ps,rv_sst,rv_pblh real(r_kind),pointer,dimension(:,:,:) :: rv_u,rv_v,rv_prse,rv_q,rv_tsen,rv_tv,rv_oz real(r_kind),pointer,dimension(:,:,:) :: rv_rank3,rv_w,rv_dw @@ -197,8 +197,10 @@ subroutine ensctl2state_ad(eval,mval,grad) call gsi_bundlegetpointer (eval(jj),'oz' ,rv_oz , istatus) call gsi_bundlegetpointer (eval(jj),'sst' ,rv_sst, istatus) + call gsi_bundlegetpointer (eval(jj),'pblh' ,rv_pblh, istatus) call gsi_bundleputvar ( wbundle_c, 'oz', rv_oz, istatus ) call gsi_bundleputvar ( wbundle_c, 'sst', rv_sst, istatus ) + call gsi_bundleputvar ( wbundle_c, 'pblh', rv_pblh, istatus ) if(wdw_exist)then call gsi_bundlegetpointer (eval(jj),'w' ,rv_w, istatus) call gsi_bundleputvar ( wbundle_c, 'w', rv_w, istatus ) diff --git a/GEOSaana_GridComp/GSI_GridComp/hybrid_ensemble_isotropic.F90 b/GEOSaana_GridComp/GSI_GridComp/hybrid_ensemble_isotropic.F90 index 0ef4d434..97059153 100644 --- a/GEOSaana_GridComp/GSI_GridComp/hybrid_ensemble_isotropic.F90 +++ b/GEOSaana_GridComp/GSI_GridComp/hybrid_ensemble_isotropic.F90 @@ -48,6 +48,7 @@ module hybrid_ensemble_isotropic ! 2015-04-07 carley - bug fix to allow grd_loc%nlat=grd_loc%nlon ! 2016-05-13 parrish - remove beta12mult ! 2018-02-15 wu - add code for fv3_regional option +! 2021-03-01 zhu - add pblh ! ! subroutines included: ! sub init_rf_z - initialize localization recursive filter (z direction) @@ -1852,6 +1853,17 @@ subroutine ensemble_forward_model(cvec,a_en,ibin) enddo enddo ! enddo n_ens + case('PBLH') + + do n=1,n_ens + do j=1,jm + do i=1,im + cvec%r2(ipic)%q(i,j)=cvec%r2(ipic)%q(i,j) & + +a_en(n)%r3(ipx)%q(i,j,1)*en_perts(n,ibin)%r2(ipic)%qr4(i,j) + enddo + enddo + enddo ! enddo n_ens + end select enddo @@ -2011,6 +2023,17 @@ subroutine ensemble_forward_model_dual_res(cvec,a_en,ibin) enddo enddo ! enddo n_ens + case('PBLH') + + do n=1,n_ens + do j=1,jm + do i=1,im + work_ens%r2(ipic)%q(i,j)=work_ens%r2(ipic)%q(i,j) & + +a_en(n)%r3(ipx)%q(i,j,1)*en_perts(n,ibin)%r2(ipic)%qr4(i,j) + enddo + enddo + enddo ! enddo n_ens + end select enddo @@ -2155,6 +2178,15 @@ subroutine ensemble_forward_model_ad(cvec,a_en,ibin) enddo enddo + case('PBLH') + + do j=1,jm + do i=1,im + a_en(n)%r3(ipx)%q(i,j,1)=a_en(n)%r3(ipx)%q(i,j,1) & + +cvec%r2(ipic)%q(i,j)*en_perts(n,ibin)%r2(ipic)%qr4(i,j) + enddo + enddo + end select enddo enddo ! enddo n_ens @@ -2311,6 +2343,15 @@ subroutine ensemble_forward_model_ad_dual_res(cvec,a_en,ibin) enddo enddo + case('PBLH') + + do j=1,jm + do i=1,im + a_en(n)%r3(ipx)%q(i,j,1)=a_en(n)%r3(ipx)%q(i,j,1) & + +work_ens%r2(ipic)%q(i,j)*en_perts(n,ibin)%r2(ipic)%qr4(i,j) + enddo + enddo + end select enddo enddo ! enddo n_ens diff --git a/GEOSaana_GridComp/GSI_GridComp/inc2guess.f90 b/GEOSaana_GridComp/GSI_GridComp/inc2guess.f90 index 14c892b7..fb3c0afe 100644 --- a/GEOSaana_GridComp/GSI_GridComp/inc2guess.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/inc2guess.f90 @@ -123,7 +123,11 @@ subroutine inc2guess(sval) if (id>0) then call gsi_bundlegetpointer (sval(ii), guess(ic),ptr2dinc,istatus) call gsi_bundlegetpointer (gsi_metguess_bundle(it),guess(ic),ptr2dges,istatus) - ptr2dges = ptr2dinc + if (trim(guess(ic))=='pblh') then + call copy_positive_fldr2_(ptr2dges,ptr2dinc) + else + ptr2dges = ptr2dinc + endif endif enddo diff --git a/GEOSaana_GridComp/GSI_GridComp/intpblh.f90 b/GEOSaana_GridComp/GSI_GridComp/intpblh.f90 index b5e9d715..9d9f8b37 100644 --- a/GEOSaana_GridComp/GSI_GridComp/intpblh.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/intpblh.f90 @@ -64,7 +64,7 @@ subroutine intpblh(pblhhead,rval,sval) use kinds, only: r_kind,i_kind use constants, only: half,one,tiny_r_kind,cg_term use obsmod, only: lsaveobsens, l_do_adjoint,luse_obsdiag - use qcmod, only: nlnqc_iter,varqc_iter + use qcmod, only: nlnqc_iter,varqc_iter,vqc use jfunc, only: jiter use gsi_bundlemod, only: gsi_bundle use gsi_bundlemod, only: gsi_bundlegetpointer @@ -87,6 +87,9 @@ subroutine intpblh(pblhhead,rval,sval) real(r_kind),pointer,dimension(:) :: rpblh type(pblhNode), pointer :: pblhptr +! If no ps data return + if(.not. associated(pblhhead))return + ! Retrieve pointers ! Simply return if any pointer not found ier=0 @@ -126,7 +129,7 @@ subroutine intpblh(pblhhead,rval,sval) if( .not. ladtest_obs) val=val-pblhptr%res ! gradient of nonlinear operator - if (nlnqc_iter .and. pblhptr%pg > tiny_r_kind .and. & + if (vqc .and. nlnqc_iter .and. pblhptr%pg > tiny_r_kind .and. & pblhptr%b > tiny_r_kind) then pg_pblh=pblhptr%pg*varqc_iter cg_pblh=cg_term/pblhptr%b @@ -136,9 +139,9 @@ subroutine intpblh(pblhhead,rval,sval) val = val*(one-p0) endif if( ladtest_obs ) then - grad = val + grad = val else - grad = val*pblhptr%raterr2*pblhptr%err2 + grad = val*pblhptr%raterr2*pblhptr%err2 end if endif diff --git a/GEOSaana_GridComp/GSI_GridComp/m_berror_stats.f90 b/GEOSaana_GridComp/GSI_GridComp/m_berror_stats.f90 index 82d4cf0d..b7de832d 100644 --- a/GEOSaana_GridComp/GSI_GridComp/m_berror_stats.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/m_berror_stats.f90 @@ -300,6 +300,7 @@ subroutine read_wgt(corz,corp,hwll,hwllp,vz,corsst,hsst,varq,qoption,varcw,cwopt ! of cw for allsky radiance ! 09Sept15 - Zhu - use centralized cloud_names_fwd and n_clouds_fwd to add ! flexibility for all-sky radiance assimilation +! 03Mar2021 - Zhu - add pblh !EOP ___________________________________________________________________ character(len=*),parameter :: myname_=myname//'::read_wgt' @@ -504,7 +505,20 @@ subroutine read_wgt(corz,corp,hwll,hwllp,vz,corsst,hsst,varq,qoption,varcw,cwopt end if endif - ! need simliar general template for undefined 2d variables ... + ! corp, hwllp for undefined 2d variables + do n=1,size(cvars2d) + if ( .not.found2d(n) ) then + if ( n>0 ) then + if ( cvars2d(n)=='pblh' ) then + do i=1,nlat + corp(i,n)=one + hwllp(i,n)=0.5_r_kind*hwll(i,1,iq) + enddo + end if + endif + if ( mype==0 ) write(6,*) myname_, ': WARNING, using general Berror template for ', cvars2d(n) + endif + enddo deallocate(found3d,found2d) diff --git a/GEOSaana_GridComp/GSI_GridComp/m_pblhNode.F90 b/GEOSaana_GridComp/GSI_GridComp/m_pblhNode.F90 index a80fed90..6f60c0a9 100644 --- a/GEOSaana_GridComp/GSI_GridComp/m_pblhNode.F90 +++ b/GEOSaana_GridComp/GSI_GridComp/m_pblhNode.F90 @@ -36,15 +36,19 @@ module m_pblhNode type,extends(obsNode):: pblhNode !type(pblh_ob_type),pointer :: llpoint => NULL() type(obs_diag), pointer :: diags => NULL() - real(r_kind) :: res ! pblh residual - real(r_kind) :: err2 ! pblh error squared - real(r_kind) :: raterr2 ! square of ratio of final obs error - ! to original obs error - !real(r_kind) :: time ! observation time in sec - real(r_kind) :: b ! variational quality control parameter - real(r_kind) :: pg ! variational quality control parameter - real(r_kind) :: wij(4) ! horizontal interpolation weights - integer(i_kind) :: ij(4) ! horizontal locations + real(r_kind) :: res =0._r_kind ! pblh residual + real(r_kind) :: err2 =0._r_kind ! pblh error squared + real(r_kind) :: raterr2=0._r_kind ! square of ratio of final obs error + ! to original obs error + !real(r_kind) :: time ! observation time in sec + real(r_kind) :: b =0._r_kind ! variational quality control parameter + real(r_kind) :: pg =0._r_kind ! variational quality control parameter + real(r_kind) :: jb =0._r_kind ! variational quality control parameter + real(r_kind) :: wij(4) =0._r_kind ! horizontal interpolation weights + real(r_kind) :: ppertb =0._r_kind ! random number adding to the obs + integer(i_kind) :: ij(4) =0._r_kind ! horizontal locations + integer(i_kind) :: kx =0_i_kind ! ob type + !logical :: luse ! flag indicating if ob is used in pen. !integer(i_kind) :: idv,iob ! device id and obs index for sorting @@ -158,6 +162,9 @@ subroutine obsNode_xread_(aNode,iunit,istat,diagLookup,skip) aNode%raterr2, & aNode%b , & aNode%pg , & + aNode%jb , & + aNode%ppertb , & + aNode%kx , & aNode%wij , & aNode%ij if (istat/=0) then @@ -191,6 +198,9 @@ subroutine obsNode_xwrite_(aNode,junit,jstat) aNode%raterr2, & aNode%b , & aNode%pg , & + aNode%jb , & + aNode%ppertb , & + aNode%kx , & aNode%wij , & aNode%ij if (jstat/=0) then diff --git a/GEOSaana_GridComp/GSI_GridComp/prewgt.f90 b/GEOSaana_GridComp/GSI_GridComp/prewgt.f90 index e7b9ef9a..d41de613 100644 --- a/GEOSaana_GridComp/GSI_GridComp/prewgt.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/prewgt.f90 @@ -73,6 +73,7 @@ subroutine prewgt(mype) ! 2014-08-02 zhu - set up new background error variance and correlation lengths of cw ! for all-sky radiance assimilation ! 2020-07-14 todling- add adjustozhscl as optional +! 2021-03-03 zhu - add pblh ! ! input argument list: ! mype - mpi task id @@ -126,7 +127,7 @@ subroutine prewgt(mype) integer(i_kind) ix,jx,mlat integer(i_kind) nf2p,istatus integer(i_kind),dimension(0:40):: iblend - integer(i_kind) nrf3_sf,nrf3_q,nrf3_vp,nrf3_t,nrf3_oz,nrf2_ps,nrf2_sst,nrf3_cw + integer(i_kind) nrf3_sf,nrf3_q,nrf3_vp,nrf3_t,nrf3_oz,nrf2_ps,nrf2_sst,nrf2_pblh,nrf3_cw integer(i_kind),allocatable,dimension(:) :: nrf3_loc,nrf2_loc real(r_kind) wlipi,wlipih,df @@ -205,6 +206,7 @@ subroutine prewgt(mype) nrf3_cw = getindex(cvars3d,'cw') nrf2_ps = getindex(cvars2d,'ps') nrf2_sst = getindex(cvars2d,'sst') + nrf2_pblh = getindex(cvars2d,'pblh') ! nrf2_stl = getindex(cvarsmd,'stl') ! nrf2_sti = getindex(cvarsmd,'sti') @@ -326,6 +328,12 @@ subroutine prewgt(mype) end do endif +! pbl height + if(nrf2_pblh>0) then + do i=1,nlat + hwllp(i,nrf2_pblh)=hwllinp(i,nrf2_pblh) + end do + endif ! sea surface temperature, convert from km to m ! also calculate a minimum horizontal length scale for @@ -488,7 +496,7 @@ subroutine prewgt(mype) end do ! Special case of dssv for qoption=2 and cw - if (qoption==2) call compute_qvar3d + if (qoption==2 .or. nrf2_pblh>0) call compute_qvar3d !!!$omp parallel do schedule(dynamic,1) private(i,n,j,jx,ix,loc) do n=1,nc2d diff --git a/GEOSaana_GridComp/GSI_GridComp/prt_guess.f90 b/GEOSaana_GridComp/GSI_GridComp/prt_guess.f90 index ad5cb800..6e4fea92 100644 --- a/GEOSaana_GridComp/GSI_GridComp/prt_guess.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/prt_guess.f90 @@ -55,7 +55,7 @@ subroutine prt_guess(sgrep) character(len=*), intent(in ) :: sgrep ! Declare local variables - integer(i_kind), parameter :: nvars=13 + integer(i_kind), parameter :: nvars=14 integer(i_kind) ii,istatus,ier,icf integer(i_kind) ntsig integer(i_kind) ntsfc @@ -72,6 +72,7 @@ subroutine prt_guess(sgrep) real(r_kind),pointer,dimension(:,:,:)::ges_oz_it=>NULL() real(r_kind),pointer,dimension(:,:,:)::ges_cwmr_it=>NULL() real(r_kind),pointer,dimension(:,:,:)::ges_cf_it=>NULL() + real(r_kind),pointer,dimension(:,: )::ges_pblh_it=>NULL() character(len=4) :: cvar(nvars+3) !******************************************************************************* @@ -103,6 +104,8 @@ subroutine prt_guess(sgrep) ier=ier+istatus if (ier/=0) return ! this is a fundamental routine, when some not found just return + call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'pblh',ges_pblh_it,istatus) + call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'cf',ges_cf_it,icf) if (icf/=0) ges_cf_it =>cfgues @@ -139,9 +142,10 @@ subroutine prt_guess(sgrep) cvar(11)='PRSL' cvar(12)='PS ' cvar(13)='SST ' - cvar(14)='radb' - cvar(15)='pcpb' - cvar(16)='aftb' + cvar(14)='PBLH' + cvar(15)='radb' + cvar(16)='pcpb' + cvar(17)='aftb' zloc(1) = sum (ges_u_it (2:lat1+1,2:lon1+1,1:nsig)) zloc(2) = sum (ges_v_it (2:lat1+1,2:lon1+1,1:nsig)) @@ -156,6 +160,7 @@ subroutine prt_guess(sgrep) zloc(11) = sum (ges_prsl (2:lat1+1,2:lon1+1,1:nsig,ntsig)) zloc(12) = sum (ges_ps_it (2:lat1+1,2:lon1+1 )) zloc(13) = sum (sfct (2:lat1+1,2:lon1+1, ntsfc)) + zloc(14) = sum (ges_pblh_it(2:lat1+1,2:lon1+1 )) zloc(nvars+1) = minval(ges_u_it (2:lat1+1,2:lon1+1,1:nsig)) zloc(nvars+2) = minval(ges_v_it (2:lat1+1,2:lon1+1,1:nsig)) zloc(nvars+3) = minval(ges_tv_it (2:lat1+1,2:lon1+1,1:nsig)) @@ -169,6 +174,7 @@ subroutine prt_guess(sgrep) zloc(nvars+11) = minval(ges_prsl (2:lat1+1,2:lon1+1,1:nsig,ntsig)) zloc(nvars+12) = minval(ges_ps_it (2:lat1+1,2:lon1+1 )) zloc(nvars+13) = minval(sfct (2:lat1+1,2:lon1+1, ntsfc)) + zloc(nvars+14) = minval(ges_pblh_it(2:lat1+1,2:lon1+1 )) zloc(2*nvars+1) = maxval(ges_u_it (2:lat1+1,2:lon1+1,1:nsig)) zloc(2*nvars+2) = maxval(ges_v_it (2:lat1+1,2:lon1+1,1:nsig)) zloc(2*nvars+3) = maxval(ges_tv_it (2:lat1+1,2:lon1+1,1:nsig)) @@ -182,6 +188,7 @@ subroutine prt_guess(sgrep) zloc(2*nvars+11) = maxval(ges_prsl (2:lat1+1,2:lon1+1,1:nsig,ntsig)) zloc(2*nvars+12) = maxval(ges_ps_it (2:lat1+1,2:lon1+1 )) zloc(2*nvars+13) = maxval(sfct (2:lat1+1,2:lon1+1, ntsfc)) + zloc(2*nvars+14) = maxval(ges_pblh_it (2:lat1+1,2:lon1+1 )) zloc(3*nvars+1) = real(lat1*lon1*nsig*ntsig,r_kind) zloc(3*nvars+2) = real(lat1*lon1*ntsig,r_kind) zloc(3*nvars+3) = real(lat1*lon1*nsig*ntsig,r_kind) diff --git a/GEOSaana_GridComp/GSI_GridComp/read_obs.F90 b/GEOSaana_GridComp/GSI_GridComp/read_obs.F90 index 7dbac635..a5b06cd9 100644 --- a/GEOSaana_GridComp/GSI_GridComp/read_obs.F90 +++ b/GEOSaana_GridComp/GSI_GridComp/read_obs.F90 @@ -166,11 +166,12 @@ subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread) integer(i_kind) ,intent(in) :: minuse integer(i_kind) ,intent(inout) :: nread - integer(i_kind) :: lnbufr,idate,idate2,iret,kidsat + integer(i_kind) :: lnbufr,idate,idate2,iret,kidsat,istat integer(i_kind) :: ireadsb,ireadmg,kx,nc,said - real(r_double) :: satid,rtype + real(r_double) :: satid,rtype,stationid character(8) subset logical,parameter:: GMAO_READ=.false. + logical pblexist satid=1 ! debug executable wants default value ??? idate=0 @@ -187,10 +188,15 @@ subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread) if(lexist .and. trim(dtype) /= 'tcp' )then lnbufr = 15 - open(lnbufr,file=trim(filename),form='unformatted',status ='unknown') - call openbf(lnbufr,'IN',lnbufr) - call datelen(10) - call readmg(lnbufr,subset,idate,iret) + if (dtype=='pblh') then + open(lnbufr,file=trim(filename),form='formatted',status ='unknown') + read(lnbufr,*, iostat=iret) idate + else + open(lnbufr,file=trim(filename),form='unformatted',status ='unknown') + call openbf(lnbufr,'IN',lnbufr) + call datelen(10) + call readmg(lnbufr,subset,idate,iret) + end if if(iret == 0)then ! Extract date and check for consistency with analysis date @@ -323,10 +329,16 @@ subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread) kidsat = 0 end if - call closbf(lnbufr) - open(lnbufr,file=trim(filename),form='unformatted',status ='unknown') - call openbf(lnbufr,'IN',lnbufr) - call datelen(10) + if (dtype=='pblh') then + close(lnbufr) + open(lnbufr,file=trim(filename),form='formatted',status ='unknown') + read(lnbufr,*) idate2 + else + call closbf(lnbufr) + open(lnbufr,file=trim(filename),form='unformatted',status ='unknown') + call openbf(lnbufr,'IN',lnbufr) + call datelen(10) + end if if(kidsat /= 0)then lexist = .false. @@ -370,6 +382,21 @@ subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread) end do nread = nread + 1 end do file2loop + else if(trim(filename) == 'wppblh')then + lexist = .false. + pblexist=.true. + wppblhloop: do while(pblexist) + read(lnbufr,iostat=istat) stationid + if (istat /=0) pblexist=.false. + kx=888 + do nc=1,nconvtype + if(trim(ioctype(nc)) == trim(dtype) .and. kx == ictype(nc) .and. icuse(nc) > minuse)then + lexist = .true. + exit wppblhloop + end if + end do + nread = nread + 1 + end do wppblhloop else if(trim(filename) == 'gps_ref' .or. trim(filename) == 'gps_bnd')then lexist = .false. gpsloop: do while(ireadmg(lnbufr,subset,idate2) >= 0) @@ -1249,6 +1276,12 @@ subroutine read_obs(ndata,mype) use_hgtl_full=.true. if(belong(i))use_hgtl_full_proc=.true. end if + if(obstype == 'pblh')then + use_hgtl_full=.true. + if(belong(i))use_hgtl_full_proc=.true. +! use_prsl_full=.true. +! if(belong(i))use_prsl_full_proc=.true. + end if if(obstype == 'sst')then if(belong(i))use_sfc=.true. endif @@ -1583,8 +1616,10 @@ subroutine read_obs(ndata,mype) ! Process pblh else if (obstype == 'pblh') then +! call read_pblh(nread,npuse,nouse,infile,obstype,lunout,twind,sis, & +! nobs_sub1(1,i)) call read_pblh(nread,npuse,nouse,infile,obstype,lunout,twind,sis, & - nobs_sub1(1,i)) + hgtl_full,nobs_sub1(1,i)) string='READ_PBLH' end if conv_obstype_select ! Process swcp and lwcp diff --git a/GEOSaana_GridComp/GSI_GridComp/read_pblh.f90 b/GEOSaana_GridComp/GSI_GridComp/read_pblh.f90 index ff7d2f22..47504289 100755 --- a/GEOSaana_GridComp/GSI_GridComp/read_pblh.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/read_pblh.f90 @@ -1,4 +1,4 @@ - subroutine read_pblh(nread,ndata,nodata,infile,obstype,lunout,twindin,& + subroutine read_pblh0(nread,ndata,nodata,infile,obstype,lunout,twindin,& sis,nobs) !$$$ subprogram documentation block ! . . . . @@ -490,5 +490,340 @@ subroutine read_pblh(nread,ndata,nodata,infile,obstype,lunout,twindin,& write(6,*)'READ_PREPFITS: closbf(',lunin,')' endif + close(lunin) + end subroutine read_pblh0 + subroutine read_pblh(nread,ndata,nodata,infile,obstype,lunout,twindin,& + sis,hgtl_full,nobs) +!$$$ subprogram documentation block +! . . . . +! subprogram: read_pblh read obs from msgs in PREPFITS files (rpf == read aircraft) +! +! program history log: +! 2009-06 whiting - coding rpf +! 2009-10-20 zhu - modify rpf for reading in pblh data in GSI +! 2009-10-21 whiting - modify cnem & pblhob for reading Caterina's files +! 2013-01-26 parrish - change from grdcrd to grdcrd1 (to allow successful debug compile on WCOSS) +! 2015-02-23 Rancic/Thomas - add l4densvar to time window logical +! 2015-10-01 guo - consolidate use of ob location (in deg) +! +! input argument list: +! infile - unit from which to read BUFR data +! obstype - observation type to process +! lunout - unit to which to write data for further processing +! hgtl_full- 3d geopotential height on full domain grid +! +! output argument list: +! nread - number of type "obstype" observations read +! nodata - number of individual "obstype" observations read +! ndata - number of type "obstype" observations retained for further processing +! twindin - input group time window (hours) +! sis - satellite/instrument/sensor indicator +! nobs - array of observations on each subdomain for each processor +! +! attributes: +! language: f90 +! machine: ibm RS/6000 SP +! +!$$$ + use kinds, only: r_kind,r_double,i_kind,r_single + use constants, only: zero,one_tenth,one,deg2rad,rad2deg,three,rearth, & + grav_equator,eccentricity,somigliana,grav_ratio,grav, & + semi_major_axis,flattening,two + use gridmod, only: diagnostic_reg,regional,nlon,nlat,nsig, & + tll2xy,txy2ll,rotate_wind_ll2xy,rotate_wind_xy2ll,& + rlats,rlons + use convinfo, only: nconvtype,ctwind, & + icuse,ictype,ioctype + use gsi_4dvar, only: l4dvar,l4densvar,time_4dvar,winlen + use obsmod, only: iadate,offtime_data,bmiss + use deter_sfc_mod, only: deter_sfc2 + use mpimod, only: npe + implicit none + +! Declare passed variables + character(len=*),intent(in):: infile,obstype + character(20),intent(in):: sis + integer(i_kind),intent(in):: lunout + integer(i_kind),intent(inout):: nread,ndata,nodata + integer(i_kind),dimension(npe),intent(inout):: nobs + real(r_kind),intent(in):: twindin + real(r_kind),dimension(nlat,nlon,nsig),intent(in):: hgtl_full + +! Declare local parameters + real(r_kind),parameter:: r360 = 360.0_r_kind + + integer(i_kind) lunin,idate + + real(r_kind),allocatable,dimension(:,:):: cdata_all + + character(10) date + logical first,outside,inflate_error,lexist,more_data + integer(i_kind) ihh,idd,iret,im,iy,istat,i,j,k + integer(i_kind) ikx,nkx,kx,nreal,ilat,ilon,iout + integer(i_kind) kk,klon1,klat1,klonp1,klatp1 + integer(i_kind) ntest,nchanl + integer(i_kind) pblhqm,maxobs,idomsfc +! integer(i_kind),dimension(8):: obs_time,anal_time,lobs_time +! real(r_kind) ltime,cenlon_tmp + real(r_kind) usage + real(r_kind) pblhob,pblhoe,pblhelev + real(r_kind) dlat,dlon,dlat_earth,dlon_earth,stnelev + real(r_kind) dlat_earth_deg,dlon_earth_deg + real(r_kind) cdist,disterr,disterrmax,rlon00,rlat00 + real(r_kind) :: tsavg,ff10,sfcr,zz +! real(r_kind),dimension(5):: tmp_time + real(r_kind) sin2,termg,termr,termrg + real(r_kind),dimension(nsig):: zges,hges + real(r_kind) dx,dy,dx1,dy1,w00,w10,w01,w11 + + integer(i_kind) idate5(5),minobs,minan + real(r_kind) time_correction,timeobs,time,toff,t4dv,zeps + +! real(r_single) stationid,lat_deg,lon_deg,altitude,localtime,utctime,localday,utcday,pblh_wp + real(r_kind) stationid,lat_deg,lon_deg,altitude,localtime,utctime,localday,utcday,pblh_wp + data lunin /50/ + +! Initialize variables + nreal=14 + ntest=0 + kx= 888 + +! Check if pblh file exists + inquire(file=trim(infile),exist=lexist) + if (.not.lexist) return + + open(lunin,file=trim(infile),form='formatted') + read(lunin, *) idate + + maxobs=0 + first=.true. + more_data=.true. + readfile1: do while(more_data) + read(lunin,*,iostat=istat) stationid, lat_deg, lon_deg, altitude, localtime, utctime, localday, utcday, pblh_wp + if (istat/=0) then + more_data=.false. + else + maxobs=maxobs+1 + end if + +! Time offset + if (first) then + call time_4dvar(idate,toff) + first=.false. + end if + end do readfile1 + + if (maxobs == 0) return + + allocate(cdata_all(nreal,maxobs)) + nread=0 + nchanl=0 + ilon=2 + ilat=3 + close(lunin) + open(lunin,file=trim(infile),form='formatted') + read(lunin, *) idate + print*, 'read_pblh idate=',idate,' maxobs=',maxobs + do i=1,maxobs + + nread=nread+1 + if(kx == 120) nkx= 120 + if(kx == 227) nkx= 181 + if(kx == 888) nkx= 888 + +! Is pblh in convinfo file + ikx=0 + do j=1,nconvtype + if(kx == ictype(j)) then + ikx=j + exit + end if + end do + if(ikx == 0) cycle + + read(lunin,*) stationid, lat_deg, lon_deg, altitude, localtime, utctime, localday, utcday, pblh_wp + print*, 'read_pblh:', stationid, lat_deg, lon_deg, altitude, localtime, utctime, localday, utcday, pblh_wp + + if(lon_deg>= r360)lon_deg=lon_deg-r360 + if(lon_deg < zero)lon_deg=lon_deg+r360 + dlon_earth_deg=lon_deg + dlat_earth_deg=lat_deg + dlon_earth=lon_deg*deg2rad + dlat_earth=lat_deg*deg2rad + if(regional)then + call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) ! convert to rotated coordinate + if(diagnostic_reg) then + call txy2ll(dlon,dlat,rlon00,rlat00) + ntest=ntest+1 + cdist=sin(dlat_earth)*sin(rlat00)+cos(dlat_earth)*cos(rlat00)* & + (sin(dlon_earth)*sin(rlon00)+cos(dlon_earth)*cos(rlon00)) + cdist=max(-one,min(cdist,one)) + disterr=acos(cdist)*rad2deg + disterrmax=max(disterrmax,disterr) + end if + if(outside) cycle ! check to see if outside regional domain + else + dlat = dlat_earth + dlon = dlon_earth + call grdcrd1(dlat,rlats,nlat,1) + call grdcrd1(dlon,rlons,nlon,1) + endif + +! Interpolate guess pressure profile to observation location + klon1= int(dlon); klat1= int(dlat) + dx = dlon-klon1; dy = dlat-klat1 + dx1 = one-dx; dy1 = one-dy + w00=dx1*dy1; w10=dx1*dy; w01=dx*dy1; w11=dx*dy + + klat1=min(max(1,klat1),nlat); klon1=min(max(0,klon1),nlon) + if (klon1==0) klon1=nlon + klatp1=min(nlat,klat1+1); klonp1=klon1+1 + if (klonp1==nlon+1) klonp1=1 + do kk=1,nsig + hges(kk)=w00*hgtl_full(klat1 ,klon1 ,kk) + & + w10*hgtl_full(klatp1,klon1 ,kk) + & + w01*hgtl_full(klat1 ,klonp1,kk) + & + w11*hgtl_full(klatp1,klonp1,kk) + end do + sin2 = sin(dlat_earth)*sin(dlat_earth) + termg = grav_equator * & + ((one+somigliana*sin2)/sqrt(one-eccentricity*eccentricity*sin2)) + termr = semi_major_axis /(one + flattening + grav_ratio - & + two*flattening*sin2) + termrg = (termg/grav)*termr + do k=1,nsig + zges(k) = (termr*hges(k)) / (termrg-hges(k)) + end do + + + if(offtime_data) then + +! in time correction for observations to account for analysis +! time being different from obs file time. + write(date,'( i10)') idate + read (date,'(i4,3i2)') iy,im,idd,ihh + idate5(1)=iy + idate5(2)=im + idate5(3)=idd + idate5(4)=ihh + idate5(5)=0 + call w3fs21(idate5,minobs) ! obs ref time in minutes relative to historic date + idate5(1)=iadate(1) + idate5(2)=iadate(2) + idate5(3)=iadate(3) + idate5(4)=iadate(4) + idate5(5)=0 + call w3fs21(idate5,minan) ! analysis ref time in minutes relative to historic date + +! Add obs reference time, then subtract analysis time to get obs time relative to analysis + + time_correction=float(minobs-minan)/60._r_kind + + else + time_correction=zero + end if + + +! Time check + timeobs=real(utctime, r_double) + time=timeobs + time_correction + t4dv=timeobs + toff + zeps=1.0e-8_r_kind + if (t4dv -zeps) t4dv=zero + if (t4dv>winlen.and.t4dvwinlen) cycle + else + if((real(abs(time)) > real(ctwind(ikx)) .or. real(abs(time)) > real(twindin))) cycle + end if + print*, 'read_pblh: idate, iadate, timeobs, toff, t4dv=', idate, iadate, timeobs, toff, t4dv + + stnelev=altitude + pblhelev=stnelev + pblhob=pblh_wp + pblhqm=0 + if (pblhob .lt. 0.0) pblhqm=15 + +! if (nkx==131 .or. nkx==133 .or. nkx==135) then +! anal_time=0 +! obs_time=0 +! tmp_time=zero +! tmp_time(2)=timeobs +! anal_time(1)=iadate(1) +! anal_time(2)=iadate(2) +! anal_time(3)=iadate(3) +! anal_time(5)=iadate(4) +! call w3movdat(tmp_time,anal_time,obs_time) ! observation time + +! lobs_time=0 +! tmp_time=zero +! cenlon_tmp=hdr(2) +! if (hdr(2) > 180.0) cenlon_tmp=hdr(2)-360.0_r_kind +! tmp_time(2)=cenlon_tmp/15.0_r_kind +! call w3movdat(tmp_time,obs_time,lobs_time) ! local observation time +! ltime = lobs_time(5)+lobs_time(6)/60.0_r_kind+lobs_time(7)/3600.0_r_kind +! if ((ltime.gt.21.0) .and. (ltime.lt.5.0)) pblhqm=3 +! end if + +! Set usage variable + usage = 0. + if(icuse(ikx) <= 0) usage=150. + if(pblhqm == 15 .or. pblhqm == 9) usage=150. + +! Set inflate_error logical + inflate_error=.false. + if (pblhqm == 3) inflate_error=.true. + +!--Outputs + + ndata=ndata+1 + nodata=nodata+1 + iout=ndata + + if(ndata > maxobs) then + write(6,*)'READ_PBLH: ***WARNING*** ndata > maxobs for ',obstype + ndata = maxobs + end if + +! Get information from surface file necessary for conventional data here + call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,tsavg,ff10,sfcr,zz) + +! setup for sort of averaged obs + pblhoe=100.0_r_kind ! temporarily + if (nkx==120) pblhoe=50. + if (nkx==181) pblhoe=100. + if (inflate_error) pblhoe=pblhoe*1.5_r_kind + + cdata_all(1,iout)=pblhoe ! pblh error (cb) + cdata_all(2,iout)=dlon ! grid relative longitude + cdata_all(3,iout)=dlat ! grid relative latitude + cdata_all(4,iout)=pblhelev ! pblh obs elevation + cdata_all(5,iout)=pblhob ! pblh obs + cdata_all(6,iout)=stationid ! station id + cdata_all(7,iout)=t4dv ! time + cdata_all(8,iout)=ikx ! type + cdata_all(9,iout)=localtime ! obs local time + cdata_all(10,iout)=pblhqm ! quality mark + cdata_all(11,iout)=usage ! usage parameter + cdata_all(12,iout)=dlon_earth_deg ! earth relative longitude (degrees) + cdata_all(13,iout)=dlat_earth_deg ! earth relative latitude (degrees) + cdata_all(14,iout)=altitude ! station elevation (m) + + end do + close(lunin) +! Normal exit + +! Write observation to scratch file + call count_obs(ndata,nreal,ilat,ilon,cdata_all,nobs) + write(lunout) obstype,sis,nreal,nchanl,ilat,ilon + write(lunout) ((cdata_all(j,i),j=1,nreal),i=1,ndata) + deallocate(cdata_all) + + if (ndata == 0) then + close(lunin) + write(6,*)'READ_pblh: close(',lunin,')' + endif + close(lunin) end subroutine read_pblh diff --git a/GEOSaana_GridComp/GSI_GridComp/setuppblh.f90 b/GEOSaana_GridComp/GSI_GridComp/setuppblh.f90 index 2ecd1ed5..6fa08555 100644 --- a/GEOSaana_GridComp/GSI_GridComp/setuppblh.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/setuppblh.f90 @@ -81,7 +81,7 @@ subroutine setuppblh(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diag use constants, only: zero,tiny_r_kind,one,half,wgtlim, & two,cg_term,pi,huge_single,r1000 use jfunc, only: jiter,last,miter - use qcmod, only: dfact,dfact1,npres_print + use qcmod, only: dfact,dfact1,npres_print,vqc use convinfo, only: nconvtype,cermin,cermax,cgross,cvar_b,cvar_pg,ictype use convinfo, only: icsubtype use m_dtime, only: dtime_setup, dtime_check @@ -109,7 +109,7 @@ subroutine setuppblh(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diag real(r_kind) pblhges,dlat,dlon,ddiff,dtime,error real(r_kind) scale,val2,ratio,ressw2,ress,residual - real(r_kind) obserrlm,obserror,val,valqc + real(r_kind) obserrlm,obserror,val,valqc,var_jb real(r_kind) term,halfpi,rwgt real(r_kind) cg_pblh,wgross,wnotgross,wgt,arg,exp_arg,rat_err2 real(r_kind) ratio_errors,tfact @@ -120,7 +120,7 @@ subroutine setuppblh(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diag real(r_single),allocatable,dimension(:,:)::rdiagbuf - integer(i_kind) ier,ilon,ilat,ipblh,id,itime,ikx,imaxerr,iqc + integer(i_kind) ier,ilon,ilat,ipblh,id,itime,ikx,iqc,ilocaltime integer(i_kind) iuse,ihgt,ilate,ilone,istnelv integer(i_kind) i,nchar,nreal,k,ii,ikxx,nn,ibin,ioff,ioff0,jj integer(i_kind) l,mm1 @@ -129,8 +129,8 @@ subroutine setuppblh(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diag integer(i_kind),dimension(nobs):: ioid ! initial (pre-distribution) obs ID logical proceed - character(8) station_id - character(8),allocatable,dimension(:):: cdiagbuf + character(15) station_id + character(15),allocatable,dimension(:):: cdiagbuf logical:: in_curbin, in_anybin type(pblhNode),pointer:: my_head @@ -142,7 +142,7 @@ subroutine setuppblh(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diag real(r_kind),allocatable,dimension(:,:,:) :: ges_z real(r_kind),allocatable,dimension(:,:,:) :: ges_pblh - equivalence(rstation_id,station_id) +! equivalence(rstation_id,station_id) type(obsLList),pointer,dimension(:):: pblhhead pblhhead => obsLL(:) @@ -170,7 +170,7 @@ subroutine setuppblh(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diag id=6 ! index of station id itime=7 ! index of observation time in data array ikxx=8 ! index of ob type - imaxerr=9 ! index of pblh max error + ilocaltime=9 ! index of pblh observation localtime iqc=10 ! index of qulaity mark iuse=11 ! index of use parameter ilone=12 ! index of longitude (degrees) @@ -183,7 +183,7 @@ subroutine setuppblh(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diag ! Check for missing data !need obs value and error do i=1,nobs - if (abs(data(ipblh,i)-bmiss) .lt. 10.0_r_kind) then + if (abs(data(ipblh,i)-bmiss) .lt. 10.0_r_kind .or. data(ipblh,i) tiny_r_kind .and. error > tiny_r_kind) then + if (vqc .and. (cvar_pg(ikx) > tiny_r_kind) .and. (error > tiny_r_kind)) then arg = exp(exp_arg) wnotgross= one-cvar_pg(ikx) cg_pblh=cvar_b(ikx) @@ -316,7 +316,7 @@ subroutine setuppblh(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diag rwgt = wgt/wgtlim else term = exp_arg - wgt = wgtlim + wgt = one rwgt = wgt/wgtlim endif valqc = -two*rat_err2*term @@ -373,6 +373,7 @@ subroutine setuppblh(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diag my_head%time = dtime my_head%b = cvar_b(ikx) my_head%pg = cvar_pg(ikx) + my_head%jb = var_jb my_head%luse = luse(i) if (luse_obsdiag) then @@ -388,6 +389,7 @@ subroutine setuppblh(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diag if(conv_diagsave .and. luse(i))then ii=ii+1 rstation_id = data(id,i) + write(station_id,*) int(data(id,i)) err_input = data(ier,i) err_adjst = data(ier,i) if (ratio_errors*error>tiny_r_kind) then diff --git a/GEOSaana_GridComp/GSI_GridComp/stppblh.f90 b/GEOSaana_GridComp/GSI_GridComp/stppblh.f90 index 08e178b4..b6898608 100644 --- a/GEOSaana_GridComp/GSI_GridComp/stppblh.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/stppblh.f90 @@ -55,7 +55,7 @@ subroutine stppblh(pblhhead,rval,sval,out,sges,nstep) ! !$$$ use kinds, only: r_kind,i_kind,r_quad - use qcmod, only: nlnqc_iter,varqc_iter + use qcmod, only: nlnqc_iter,varqc_iter,vqc use constants, only: half,one,two,tiny_r_kind,cg_term,zero_quad use gsi_bundlemod, only: gsi_bundle use gsi_bundlemod, only: gsi_bundlegetpointer @@ -117,7 +117,7 @@ subroutine stppblh(pblhhead,rval,sval,out,sges,nstep) end if ! Modify penalty term if nonlinear QC - if (nlnqc_iter .and. pblhptr%pg > tiny_r_kind .and. & + if (vqc .and. nlnqc_iter .and. pblhptr%pg > tiny_r_kind .and. & pblhptr%b > tiny_r_kind) then pg_pblh=pblhptr%pg*varqc_iter cg_pblh=cg_term/pblhptr%b