From a95d8ec27d7b8f249d9269d93bc9935cf2a3caa2 Mon Sep 17 00:00:00 2001 From: David Dowell Date: Mon, 10 Jul 2023 15:44:12 +0000 Subject: [PATCH 1/7] Trying to add new FED observation operator. --- src/gsi/setupfed.f90 | 1152 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1152 insertions(+) create mode 100644 src/gsi/setupfed.f90 diff --git a/src/gsi/setupfed.f90 b/src/gsi/setupfed.f90 new file mode 100644 index 0000000000..6432791443 --- /dev/null +++ b/src/gsi/setupfed.f90 @@ -0,0 +1,1152 @@ +module fed_setup + implicit none + private + public:: setup + interface setup; module procedure setupfed; end interface + +contains +subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsave,init_pass) +!$$$ subprogram documentation block +! . . . . +! subprogram: setupfed compute rhs of flash extent density +! orig. prgmmr: +! Rong Kong CAPS/OU 2018-01-21 (modified based on setupdbz.f90) +! modified: +! Yaping Wang CIMMS/OU 2019-11-11 +! David Dowell (DCD) NOAA GSL 2021-07-01 +! - added a second option (tanh) for observation operator, based on the +! work of Sebok and Back (2021, unpublished) +! - capped maximum model FED +! +! + use mpeu_util, only: die,perr + use kinds, only: r_kind,r_single,r_double,i_kind + use m_obsdiagNode, only: obs_diag + use m_obsdiagNode, only: obs_diags + use m_obsdiagNode, only: obsdiagLList_nextNode + use m_obsdiagNode, only: obsdiagNode_set + use m_obsdiagNode, only: obsdiagNode_get + use m_obsdiagNode, only: obsdiagNode_assert + use obsmod, only: rmiss_single,& + lobsdiagsave,nobskeep,lobsdiag_allocated,time_offset + use obsmod, only: oberror_tune + use m_obsNode, only: obsNode + use m_fedNode, only: fedNode + use m_fedNode, only: fedNode_appendto + use obsmod, only: luse_obsdiag, netcdf_diag, binary_diag, dirname, ianldate + use nc_diag_write_mod, only: nc_diag_init, nc_diag_header, nc_diag_metadata, & + nc_diag_write, nc_diag_data2d + use nc_diag_read_mod, only: nc_diag_read_init, nc_diag_read_get_dim,nc_diag_read_close + use m_obsLList, only: obsLList + use gsi_4dvar, only: nobs_bins,hr_obsbin + use oneobmod, only: oneobtest,maginnov,magoberr + use guess_grids, only: hrdifsig,nfldsig,ges_prsi + use guess_grids, only: ges_lnprsl, ges_prsl, ges_tsen, geop_hgtl + use gridmod, only: lat2, lon2 + use gridmod, only: nsig, get_ij,get_ijk,regional,tll2xy + use constants, only: flattening,semi_major_axis,grav_ratio,zero,grav,wgtlim + use constants, only: half,one,two,grav_equator,eccentricity,somigliana + use constants, only: rad2deg,deg2rad,r60,tiny_r_kind,cg_term,huge_single + use constants, only: r10,r100,r1000 + use constants, only: rd,grav,tpwcon + use qcmod, only: npres_print,ptop,pbot,ptopq,pbotq,dfact,dfact1 + use jfunc, only: jiter,last,jiterstart,miter + use convinfo, only: nconvtype,cermin,cermax,cgross,cvar_b,cvar_pg,ictype + use convinfo, only: icsubtype + use converr, only: ptabl + use m_dtime, only: dtime_setup, dtime_check, dtime_show + use state_vectors, only: nsdim + + use gsi_bundlemod, only: GSI_BundleGetPointer + use gsi_metguess_mod, only: gsi_metguess_get, GSI_MetGuess_Bundle + + use netcdf + + + implicit none +! Declare passed variables + type(obsLList ),target,dimension(:),intent(in):: obsLL + type(obs_diags),target,dimension(:),intent(in):: odiagLL + + logical ,intent(in ) :: fed_diagsave + integer(i_kind) ,intent(in ) :: lunin,mype,nele,nobs + real(r_kind),dimension(100+7*nsig) ,intent(inout) :: awork + real(r_kind),dimension(npres_print,nconvtype,5,3),intent(inout) :: bwork + integer(i_kind) ,intent(in ) :: is ! ndat index + logical ,intent(in ) :: init_pass ! state of "setup" parameters + + +! Declare local parameters + integer(i_kind),parameter:: fed_obs_ob_shape = 2 ! 1 = linear (Allen et al.) + ! 2 = tanh (Sebok and Back) +! coefficients for tanh operator, from Sebok and Back (2021) +! real(r_kind),parameter:: a_coeff = 8.4_r_kind ! a (flashes/min) in tanh operator +! real(r_kind),parameter:: b_coeff = 12.248_r_kind ! b (flashes/min) in tanh operator +! real(r_kind),parameter:: c_coeff = 5.0e-10_r_kind ! c (radians/kg) in tanh operator +! real(r_kind),parameter:: d_coeff = 1.68e9_r_kind ! d (kg) in tanh operator +! real(r_kind),parameter:: fed_highbnd = 18.0_r_kind ! DCD: Sebok and Back (2021, unpublished) + +! coefficients for tanh operator, from work by A. Back with regional FV3 (2023) + real(r_kind),parameter:: a_coeff = -3.645_r_kind ! a (flashes/min) in tanh operator + real(r_kind),parameter:: b_coeff = 15.75_r_kind ! b (flashes/min) in tanh operator + real(r_kind),parameter:: c_coeff = 1.939e-10_r_kind ! c (radians/kg) in tanh operator + real(r_kind),parameter:: d_coeff = -1.215e9_r_kind ! d (kg) in tanh operator + real(r_kind),parameter:: fed_highbnd = 8.0_r_kind ! DCD: Back (2023, unpublished) for FV3 + + real(r_kind),parameter:: fed_height = 6500.0_r_kind ! assumed height (m) of FED observations + real(r_kind),parameter:: r0_001 = 0.001_r_kind + real(r_kind),parameter:: r8 = 8.0_r_kind + real(r_kind),parameter:: ten = 10.0_r_kind + real(r_kind),parameter:: D608=0.608_r_kind + character(len=*),parameter:: myname='setupfed' + + real(r_kind) :: Cs_tmp, Cg_tmp ! temporary coefficients for check-up + +! Declare external calls for code analysis + external:: tintrp2a1 + external:: tintrp2a11 + external:: tintrp2a1116 + external:: tintrp31 + external:: grdcrd1 + external:: stop2 + +! Declare local variables + real(r_kind) rlow,rhgh,rsig + real(r_kind) dz,denom + real(r_kind) jqg_num,jqg + real(r_kind) wgt_dry, wgt_wet + real(r_kind) jqg_num_dry, jqg_num_wet + real(r_kind) dlnp,pobl,zob + real(r_kind) sin2,termg,termr,termrg + real(r_kind) psges,zsges + real(r_kind),dimension(nsig):: zges,hges + real(r_kind) prsltmp(nsig) + real(r_kind) sfcchk + real(r_kind) residual,obserrlm,obserror,ratio,scale,val2 + real(r_kind) ress,ressw + real(r_kind) val,valqc,rwgt + real(r_kind) cg_w,wgross,wnotgross,wgt,arg,exp_arg,term,rat_err2 + real(r_double) rstation_id + real(r_kind) dlat,dlon,dtime,dpres,ddiff,error,slat,dlat8km,dlon8km + real(r_kind) ratio_errors + real(r_kind) qgges,rhoges + real(r_kind) Ze,rdBZ,presw,fednoise,fednoise_runits + real(r_kind) Ze_orig, Zer, Zes, Zeg + real(r_kind) Zeg_dry, Zeg_wet + real(r_kind) errinv_input,errinv_adjst,errinv_final + real(r_kind) err_input,err_adjst,err_final + real(r_kind) qgexp + real(r_kind),dimension(nele,nobs):: data + real(r_kind),dimension(lat2,lon2,nfldsig)::rp + real(r_single),allocatable,dimension(:,:)::rdiagbuf + real(r_kind),allocatable,dimension(:,:,: ) :: ges_z + real(r_kind),allocatable,dimension(:,:,: ) :: ges_ps + real(r_kind),allocatable,dimension(:,:,:,:) :: ges_q + real(r_kind),allocatable,dimension(:,:,:,:) :: ges_qg,ges_qg_mask + + real(r_kind) :: presq + real(r_kind) :: P1D,T1D,Q1D,RHO + real(r_kind) :: qges,tsenges ! used to calculate tv - virtual temperature + real(r_kind) :: lnprslges ! use log(p) for vertical interpolation + real(r_kind) :: qg_min + real(r_kind) :: glmcoeff = 2.088_r_kind*10.0**(-8.0) ! Allen et al. (2016,MWR) + real(r_kind) :: CM = 0.5_r_kind ! tuning factor in eq. 14 of Kong et al. 2020 + + integer(i_kind) i,nchar,nreal,k,j,k1,ii,nii,jj,im,jm,km + integer(i_kind) mm1,k2,isli + integer(i_kind) jsig,ikxx,nn,ibin,ioff,ioff0 + integer(i_kind) ier,ilat,ilon,ihgt,ifedob,ikx,itime,iuse + integer(i_kind) ielev,id,itilt,iazm,ilone,ilate,irange + integer(i_kind) ier2,ifednoise,it,istatus + integer(i_kind) ier_b + integer(i_kind) ijk + + integer(i_kind) i4,j4,k4,n4 + integer(i_kind) nlat_ll,nlon_ll,nsig_ll,nfld_ll + + integer(i_kind) ipres,iqmax,iqc,icat,itemp + integer(i_kind) istnelv,iobshgt,izz,iprvd,isprvd,iptrb + integer(i_kind) idomsfc,iskint,isfcr,iff10 + integer(i_kind) nguess + + character(8) station_id + character(8),allocatable,dimension(:):: cdiagbuf + character(80):: string + character(128):: diag_file + logical :: diagexist + logical,dimension(nobs):: luse,muse + integer(i_kind),dimension(nobs):: ioid ! initial (pre-distribution) obs ID + logical proceed + + equivalence(rstation_id,station_id) + real(r_kind) wrange + integer(i_kind) numequal,numnotequal,kminmin,kmaxmax,istat + + logical:: in_curbin, in_anybin + + type(fedNode),pointer:: my_head + type(obs_diag),pointer:: my_diag + type(obs_diags),pointer:: my_diagLL + + real(r_kind),dimension(nsig+1):: prsitmp + + +!------------------------------------------------! + + integer(i_kind) :: ibgn, iend, jbgn, jend, ips,ipe,jps,jpe,itmp,jtmp,ktmp + character(256) :: binfilename + + integer(i_kind), parameter :: ntimesfed=1 + character(256) ::fedfilename + integer(i_kind),parameter :: nxfed=99, nyfed=99, nzfed=1, nfldfed=3 + real(4) :: a(ntimesfed,nfldfed,nzfed,nyfed,nxfed) + real(4) :: gga(ntimesfed,nfldfed,nzfed,nyfed,nxfed) + integer(i_kind) irec1, irec2, irec3, irec4, itot + integer(i_kind) :: la, iobs, lt, nnnnn + real(r_kind),dimension(nobs) :: FEDMdiag,FEDMdiagTL + real(r_kind),dimension(nobs) :: FEDMdiag2D + integer(i_kind) :: npt + integer(i_kind) :: nobsfed + real(r_kind) :: dlat_earth,dlon_earth + logical :: outside + +! YPW added the next lines + logical :: l_set_oerr_ratio_fed=.False. + logical :: l_gpht2gmht = .True. + integer(i_kind) :: ncid,status,x_dimid,y_dimid,z_dimid,varid,x_varid,y_varid + integer(i_kind),dimension(3):: dimids + integer(i_kind),dimension(2):: dimids_2d + character(256) :: outfile + real(r_kind),dimension(nobs) :: dlatobs,dlonobs + real(r_kind),dimension(4):: wgrd + integer(i_kind),dimension(4):: jgrd + integer(i_kind):: ngx,ngy,igx,jgy + real(r_kind):: dx_m, dy_m + + type(obsLList),pointer,dimension(:):: fedhead + fedhead => obsLL(:) + +!============================================================================================ +! Read and reformat observations in work arrays. + read(lunin)data,luse,ioid + + write(6,*)myname,'(pe=',mype,') nele nobs =',nele,nobs, & + 'luse_obsdiag=',luse_obsdiag,'lat2,lon2=',lat2,lon2 + + ier=1 ! index of obs error + ilon=2 ! index of grid relative obs location (x) + ilat=3 ! index of grid relative obs location (y) + ipres=4 ! index of pressure + ifedob=5 ! index of fed observation + id=6 ! index of station id + itime=7 ! index of observation time in data array + ikxx=8 ! index of ob type + iqmax=9 ! index of max error + itemp=10 ! index of dry temperature + iqc=11 ! index of quality mark + ier2=12 ! index of original-original obs error ratio + iuse=13 ! index of use parameter + idomsfc=14 ! index of dominant surface type + iskint=15 ! index of surface skin temperature + iff10=16 ! index of 10 meter wind factor + isfcr=17 ! index of surface roughness + ilone=18 ! index of longitude (degrees) + ilate=19 ! index of latitude (degrees) + istnelv=20 ! index of station elevation (m) + iobshgt=21 ! index of observation height (m) + izz=22 ! index of surface height + iprvd=23 ! index of observation provider + isprvd=24 ! index of observation subprovider + icat =25 ! index of data level category + iptrb=26 ! index of fed perturbation + do i=1,nobs + muse(i)=nint(data(iuse,i)) <= jiter + end do + + numequal=0 + numnotequal=0 + +! +! If requested, save select data for output to diagnostic file + if(fed_diagsave)then + ii=0 + nchar=1_i_kind + ioff0=26_i_kind ! 21 + 5 (22->Zr; 23->Zs; 24->Zg; 25->tsenges;26->RHO;) + nreal=ioff0 + if (lobsdiagsave) nreal=nreal+4*miter+1 + allocate(cdiagbuf(nobs),rdiagbuf(nreal,nobs)) + rdiagbuf=zero + if(netcdf_diag) call init_netcdf_diag_ + end if + mm1=mype+1 + scale=one + rsig=nsig + + + + !============================================================================================ +! +! Check to see if required guess fields are available +! vars. list: ps, z, q +! vars. list: qr, qs, qg + !============================================================================================ + + call check_vars_(proceed) + if(.not.proceed) then + write(6,*) myname,': some or all necessary variables are not available for fed obs operator. Quit!' + return ! not all vars available, simply return + end if + +! If require guess vars available, extract from bundle ... + call init_vars_ +! qscalar=zero + + !============================================================================================ + ! 1) Calculate the graupel-mass and graupel-volume based flash extent density + ! (FED) on model space, added by R. Kong, 07/05/2018 + !============================================================================================ + ges_qg_mask=ges_qg + where(ges_qg>0.0005) !Count the volume where qg > 0.5/kg + ges_qg_mask=1.0 + elsewhere + ges_qg_mask=0.0 + endwhere + + ! Operator start here + ! set ngx and ngy =2, so the integrated domain is 15kmx15km + ngx = 2 + ngy = 2 + dx_m = 3000. + dy_m = 3000. + print*,'Operator start here!,ngx=',ngx,'ngy=',ngy + rp=zero + + print*, 'mype = ', mype + print*, 'nfldsig = ', nfldsig + print*, 'nsig = ', nsig + print*, 'lon2 = ', lon2 + print*, 'lat2 = ', lat2 + +! compute graupel mass, in kg per 15 km x 15 km column + do jj=1,nfldsig + do k=1,nsig + do i=1,lon2 + do j=1,lat2 !How to handle MPI???? + do igx=1,2*ngx+1 !horizontal integration of qg around point to calculate FED + do jgy=1,2*ngy+1 + itmp = i-ngx+igx-1 + jtmp = j-ngy+jgy-1 + itmp = min(max(1,itmp),lon2) + jtmp = min(max(1,jtmp),lat2) + rp(j,i,jj)=rp(j,i,jj) + ges_qg(jtmp,itmp,k,jj)* & + dx_m*dy_m*(ges_prsi(jtmp,itmp,k,jj)-ges_prsi(jtmp,itmp,k+1,jj))*& + tpwcon * r10 + enddo !igx + enddo !jgy + end do !j + end do !i + end do !k + end do !jj + +! compute FED, in flashes/min + do jj=1,nfldsig + do i=1,lon2 + do j=1,lat2 + if (fed_obs_ob_shape .eq. 1) then + rp(j,i,jj) = CM * glmcoeff * rp(j,i,jj) + else if (fed_obs_ob_shape .eq. 2) then + rp(j,i,jj) = a_coeff + b_coeff & + * tanh(c_coeff * (rp(j,i,jj) - d_coeff)) + else + write(6,*) ' unknown fed_obs_ob_shape: ', fed_obs_ob_shape + write(6,*) ' aborting setupfed' + call stop2(999) + end if + if (rp(j,i,jj) .gt. fed_highbnd) rp(j,i,jj) = fed_highbnd + end do !j + end do !i + end do !jj + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + write(6,*) 'fed_obs_ob_shape=',fed_obs_ob_shape + if (fed_obs_ob_shape .eq. 2) then + write(6,*) 'a_coeff=',a_coeff + write(6,*) 'b_coeff=',b_coeff + write(6,*) 'c_coeff=',c_coeff + write(6,*) 'd_coeff=',d_coeff + end if + write(6,*) 'fed_highbnd=',fed_highbnd + write(6,*) 'maxval(ges_qg)=',maxval(ges_qg),'pe=',mype +! write(6,*) 'maxval(geop_hgtl)=',maxval(geop_hgtl(:,:,:,it)) + write(6,*) 'maxval(ges_tsen)=',maxval(ges_tsen(:,:,:,it)) + write(6,*) 'maxval(FED)=',maxval(rp) + write(6,*) 'ges_prsi',ges_prsi(100,100,1,1),ges_prsi(100,100,nsig,1) + + + !============================================================================================ + + nlat_ll=size(ges_qg,1) + nlon_ll=size(ges_qg,2) + nsig_ll=size(ges_qg,3) + nfld_ll=size(ges_qg,4) + + +! - Observation times are checked in read routine - comment out for now + +! call dtime_setup() + +!print*,"maxval(data(ifedob,:)),mmaxval(data(ilat,:))=",minval(data(ifedob,:)),maxval(data(ifedob,:)),maxval(data(ilat,:)) +!write(6,*) "OKOKOKOKOK, nobs=", nobs + do i=1,nobs + dtime=data(itime,i) + dlat=data(ilat,i) + dlon=data(ilon,i) + + dlon8km=data(iprvd,i) !iprvd=23 + dlat8km=data(isprvd,i) !isprvd=24 + + dpres=data(ipres,i) ! from rdararef_mosaic2: this height abv MSL + ikx = nint(data(ikxx,i)) + error=data(ier2,i) + slat=data(ilate,i)*deg2rad ! needed when converting geophgt to + dlon_earth = data(ilone,i) !the lontitude and latitude on the obs pts. + dlat_earth = data(ilate,i) + ! geometric hgh (hges --> zges below) +! print*,'i,mype,dlat,dlon,dlon8km,dlat8km',i,mype,dlat,dlon,dlon8km,dlat8km,& +! dlon_earth,dlat_earth,dpres,data(ifedob,i) + + if (nobs_bins>1) then + ibin = NINT( dtime/hr_obsbin ) + 1 + else + ibin = 1 + endif + + IF (ibin<1.OR.ibin>nobs_bins) write(6,*)mype,'Error nobs_bins,ibin= ',nobs_bins,ibin + + if (luse_obsdiag) my_diagLL => odiagLL(ibin) + +! Link obs to diagnostics structure + if(luse_obsdiag)then + my_diag => obsdiagLList_nextNode(my_diagLL ,& + create = .not.lobsdiag_allocated ,& + idv = is ,& + iob = ioid(i) ,& + ich = 1 ,& + elat = data(ilate,i) ,& + elon = data(ilone,i) ,& + luse = luse(i) ,& + miter = miter ) + if(.not.associated(my_diag)) call die(myname, & + 'obsdiagLList_nextNode(), create =',.not.lobsdiag_allocated) + endif + +! Interpolate terrain height(model elevation) to obs location. + call tintrp2a11(ges_z,zsges,dlat,dlon,dtime,hrdifsig,& + mype,nfldsig) +! print*,'i,after tintrp2all',i,mype,dlat,zsges +! 1. dpres (MRMS obs height is height above MSL) is adjusted by zsges, so it +! is changed to height relative to model elevation (terrain). +! because in GSI, geop_hgtl is the height relative to terrain (ges_z) +! (subroutine guess_grids) + dpres=dpres-zsges + if (dpres rsig)ratio_errors = zero + +!----------------------------------------------------------------------------! +! ! +! Implementation of forward operator for flash extend densit ----------------! +! ! +!----------------------------------------------------------------------------! + + !============================================================================================ + ! 3) H(x), interpolate the FED from model space on the local domain to obs space (FEDMdiag) + !============================================================================================ + + npt = 0 + FEDMdiag(i) = 0. + call tintrp2a11(rp,FEDMdiag(i),dlat,dlon,dtime,hrdifsig,mype,nfldsig) + dlonobs(i) = dlon_earth + dlatobs(i) = dlat_earth + + ! also Jacobian used for TLM and ADM + !FEDMdiagTL, used for gsi-3dvar,will be implemented in future...... + FEDMdiagTL(i) = 0. + jqg_num = FEDMdiagTL(i) !=dFED/Dqg + jqg = jqg_num + + + !end select + + if(FEDMdiag(i)==data(ifedob,i)) then + numequal=numequal+1 + else + numnotequal=numnotequal+1 + end if + +!!!!!!!!!!!!!!!!!END H(x)!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ! Compute innovations + !--------------Calculate departure from observation----------------! + + ddiff = data(ifedob,i) - FEDMdiag(i) + +! If requested, setup for single obs test. +! Note: do not use this way to run single obs test for fed in the current version. (g.zhao) + if (oneobtest) then + ddiff=maginnov +! if (trim(adjustl(oneob_type))=='fed') then +! data(ifedob,i) = maginnov +! ddiff = data(ifedob,i) - FEDMdiag(i) +! end if + error=one/(magoberr) + ratio_errors=one + end if + + +! Gross error checks + obserror = one/max(ratio_errors*error,tiny_r_kind) + obserrlm = max(cermin(ikx),min(cermax(ikx),obserror)) + + residual = abs(ddiff) != y-H(xb) + ratio = residual/obserrlm != y-H(xb)/sqrt(R) + + if (l_set_oerr_ratio_fed) then + if (ratio > cgross(ikx) .or. ratio_errors < tiny_r_kind) then + if (luse(i)) awork(4) = awork(4)+one + error = zero + ratio_errors = zero + end if + else + ratio_errors = one + end if + + if (ratio_errors*error <=tiny_r_kind) muse(i)=.false. +! if (nobskeep>0 .and. luse_obsdiag) muse(i)=obsdiags(i_fed_ob_type,ibin)%tail%muse(nobskeep) + if (nobskeep>0.and.luse_obsdiag) call obsdiagNode_get(my_diag, jiter=nobskeep, muse=muse(i)) + + val = error*ddiff !=y-H(xb)/sqrt(R) + +! Compute penalty terms (linear & nonlinear qc). + if(luse(i))then + exp_arg = -half*val**2 + rat_err2 = ratio_errors**2 + val2=val*val !(o-g)**2/R, would be saved in awork + if (cvar_pg(ikx) > tiny_r_kind .and. error > tiny_r_kind) then + arg = exp(exp_arg) + wnotgross= one-cvar_pg(ikx) + cg_w=cvar_b(ikx) + wgross = cg_term*cvar_pg(ikx)/(cg_w*wnotgross) + term = log((arg+wgross)/(one+wgross)) + wgt = one-wgross/(arg+wgross) + rwgt = wgt/wgtlim + else + term = exp_arg + wgt = wgtlim + rwgt = wgt/wgtlim + endif + valqc = -two*rat_err2*term + +! print*,'Compute penalty terms' +! Accumulate statistics for obs belonging to this task + if (muse(i)) then + if(rwgt < one) awork(21) = awork(21)+one + jsig = dpres + jsig=max(1,min(jsig,nsig)) + awork(6*nsig+jsig+100)=awork(6*nsig+jsig+100)+val2*rat_err2 + awork(5*nsig+jsig+100)=awork(5*nsig+jsig+100)+one + awork(3*nsig+jsig+100)=awork(3*nsig+jsig+100)+valqc + end if +! Loop over pressure level groupings and obs to accumulate +! statistics as a function of observation type. + ress = scale*ddiff + ressw = ress*ress + nn=1 + if (.not. muse(i)) then + nn=2 + if(ratio_errors*error >=tiny_r_kind)nn=3 + end if + do k = 1,npres_print +! if(presw >=ptop(k) .and. presw<=pbot(k))then + if(presq >=ptopq(k) .and. presq<=pbotq(k))then + bwork(k,ikx,1,nn) = bwork(k,ikx,1,nn)+one ! count + bwork(k,ikx,2,nn) = bwork(k,ikx,2,nn)+ddiff ! bias + bwork(k,ikx,3,nn) = bwork(k,ikx,3,nn)+ressw ! (o-g)**2 + bwork(k,ikx,4,nn) = bwork(k,ikx,4,nn)+val2*rat_err2 ! penalty + bwork(k,ikx,5,nn) = bwork(k,ikx,5,nn)+valqc ! nonlin qc penalty + + end if + end do + end if + + if(luse_obsdiag)then + call obsdiagNode_set(my_diag, luse=luse(i), wgtjo=(error*ratio_errors)**2, & + jiter=jiter, muse=muse(i), nldepart=ddiff) + end if + +! If obs is "acceptable", load array with obs info for use +! in inner loop minimization (int* and stp* routines) + if ( .not. last .and. muse(i)) then + + allocate(my_head) ! YPW added + call fedNode_appendto(my_head,fedhead(ibin)) + + my_head%idv=is + my_head%iob=i + my_head%elat= data(ilate,i) + my_head%elon= data(ilone,i) + + if (istat/=0) write(6,*)'MAKECOBS: allocate error for fedtail_dzg,istat=',istat + + my_head%dlev= dpres + call get_ijk(mm1,dlat,dlon,dpres,my_head%ij,my_head%wij) + my_head%res = ddiff ! Observation - ges + my_head%err2 = error**2 + my_head%raterr2 = ratio_errors**2 +! my_head%jqg = jqg ! for TL and ADJ + my_head%time = dtime + my_head%b = cvar_b(ikx) + my_head%pg = cvar_pg(ikx) + my_head%luse = luse(i) + + if(oberror_tune) then + ! my_head%fedpertb=data(iptrb,i)/error/ratio_errors + my_head%kx=ikx + if(presq > ptabl(2))then + my_head%k1=1 + else if( presq <= ptabl(33)) then + my_head%k1=33 + else + k_loop: do k=2,32 + if(presq > ptabl(k+1) .and. presq <= ptabl(k)) then + my_head%k1=k + exit k_loop + endif + enddo k_loop + endif + endif +!------------------------------------------------- + + if(luse_obsdiag)then + call obsdiagNode_assert(my_diag,my_head%idv,my_head%iob,1,myname,'my_diag:my_head') + my_head%diags => my_diag + endif + + my_head => null() + endif + +! Save select output for diagnostic file + if(.not.luse(i))write(6,*)' luse, mype',luse(i),mype + if(fed_diagsave .and. luse(i))then + ii=ii+1 + rstation_id = data(id,i) + err_input = data(ier2,i) + err_adjst = data(ier,i) + if (ratio_errors*error>tiny_r_kind) then + err_final = one/(ratio_errors*error) + else + err_final = huge_single + endif + errinv_input = huge_single + errinv_adjst = huge_single + errinv_final = huge_single + if (err_input>tiny_r_kind) errinv_input = one/err_input + if (err_adjst>tiny_r_kind) errinv_adjst = one/err_adjst + if (err_final>tiny_r_kind) errinv_final = one/err_final + + if(binary_diag) call contents_binary_diag_(my_diag) + if(netcdf_diag) call contents_netcdf_diag_(my_diag) + + end if + end do + +! Release memory of local guess arrays + call final_vars_ + +! Write information to diagnostic file + if(fed_diagsave .and. netcdf_diag) call nc_diag_write + if(fed_diagsave .and. binary_diag .and. ii>0)then + + write(string,600) jiter +600 format('fed_',i2.2) + diag_file=trim(dirname) // trim(string) + if(init_pass) then + open(66,file=trim(diag_file),form='unformatted',status='unknown',position='rewind') + else + inquire(file=trim(diag_file),exist=diagexist) + if (diagexist) then + open(66,file=trim(diag_file),form='unformatted',status='old',position='append') + else + open(66,file=trim(diag_file),form='unformatted',status='unknown',position='rewind') + endif + endif + if(init_pass .and. mype == 0) then + write(66) ianldate + write(6,*)'SETUPFED: write time record to file ',& + trim(diag_file), ' ',ianldate + endif + +! call dtime_show(myname,'diagsave:fed',i_fed_ob_type) + write(66)'fed',nchar,nreal,ii,mype,ioff0 + write(66)cdiagbuf(1:ii),rdiagbuf(:,1:ii) + ! write(6,*)'fed,nchar,nreal,ii,mype',nchar,nreal,ii,mype + deallocate(cdiagbuf,rdiagbuf) + close(66) + end if + +! End of routine + + +! return + + contains + + subroutine check_vars_ (proceed) + + + logical,intent(inout) :: proceed + integer(i_kind) ivar, istatus +! Check to see if required guess fields are available + call gsi_metguess_get ('var::ps', ivar, istatus ) + proceed=ivar>0 + print*,'For ps, proceed=',proceed + call gsi_metguess_get ('var::z' , ivar, istatus ) + proceed=proceed.and.ivar>0 + print*,'For z, proceed=',proceed + call gsi_metguess_get ('var::q' , ivar, istatus ) + proceed=proceed.and.ivar>0 + print*,'For q, proceed=',proceed +! call gsi_metguess_get ('var::tv' , ivar, istatus ) +! proceed=proceed.and.ivar>0 + call gsi_metguess_get ('var::qs', ivar, istatus ) + proceed=proceed.and.ivar>0 + print*,'For qs, proceed=',proceed + call gsi_metguess_get ('var::qg', ivar, istatus ) + proceed=proceed.and.ivar>0 + print*,'For qg, proceed=',proceed + call gsi_metguess_get ('var::qr', ivar, istatus ) + proceed=proceed.and.ivar>0 + print*,'For qr, proceed=',proceed +! if ( mphyopt == 108 ) then ! comment out by YPW +! proceed=proceed.and.ivar>0 ! comment out by YPW +! print*,'For qnr, proceed=',proceed ! comment out by YPW +! end if ! comment out by YPW + end subroutine check_vars_ + + + subroutine init_vars_ + +! use radaremul_cst, only: mphyopt + + real(r_kind),dimension(:,: ),pointer:: rank2=>NULL() + real(r_kind),dimension(:,:,:),pointer:: rank3=>NULL() + character(len=5) :: varname + integer(i_kind) ifld, istatus + +! If require guess vars available, extract from bundle ... + if(size(gsi_metguess_bundle)==nfldsig) then +! get ps ... + varname='ps' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus) + if (istatus==0) then + if(allocated(ges_ps))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_ps(size(rank2,1),size(rank2,2),nfldsig)) + ges_ps(:,:,1)=rank2 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank2,istatus) + ges_ps(:,:,ifld)=rank2 + enddo + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif +! get z ... + varname='z' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus) + if (istatus==0) then + if(allocated(ges_z))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_z(size(rank2,1),size(rank2,2),nfldsig)) + ges_z(:,:,1)=rank2 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank2,istatus) + ges_z(:,:,ifld)=rank2 + enddo + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif +! get q ... + varname='q' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) + if (istatus==0) then + if(allocated(ges_q))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_q(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) + ges_q(:,:,:,1)=rank3 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) + ges_q(:,:,:,ifld)=rank3 + enddo + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif +! get tv ... +! varname='tv' +! call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) +! if (istatus==0) then +! if(allocated(ges_tv))then +! write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' +! call stop2(999) +! endif +! allocate(ges_tv(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) +! ges_tv(:,:,:,1)=rank3 +! do ifld=2,nfldsig +! call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) +! ges_tv(:,:,:,ifld)=rank3 +! ges_tv(:,:,:,ifld)=rank3 +! enddo +! else +! write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus +! call stop2(999) +! endif +! get qr ... +! get qg ... + varname='qg' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) + if (istatus==0) then + if(allocated(ges_qg))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_qg(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) + if(.not. allocated(ges_qg_mask))then + allocate(ges_qg_mask(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) + endif + + ges_qg(:,:,:,1)=rank3 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) + ges_qg(:,:,:,ifld)=rank3 + enddo + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif + else + write(6,*) trim(myname), ': inconsistent vector sizes (nfldsig,size(metguess_bundle) ',& + nfldsig,size(gsi_metguess_bundle) + call stop2(999) + endif + end subroutine init_vars_ + + subroutine init_netcdf_diag_ + character(len=80) string + character(len=128) diag_conv_file + integer(i_kind) ncd_fileid,ncd_nobs + logical append_diag + logical,parameter::verbose=.false. + write(string,900) jiter +900 format('conv_fed_',i2.2,'.nc4') + diag_conv_file=trim(dirname) // trim(string) + + inquire(file=diag_conv_file, exist=append_diag) + + if (append_diag) then + call nc_diag_read_init(diag_conv_file,ncd_fileid) + ncd_nobs = nc_diag_read_get_dim(ncd_fileid,'nobs') + call nc_diag_read_close(diag_conv_file) + + if (ncd_nobs > 0) then + if(verbose) print *,'file ' // trim(diag_conv_file) // ' exists. Appending. nobs,mype=',ncd_nobs,mype + else + if(verbose) print *,'file ' // trim(diag_conv_file) // ' exists but contains no obs. Not appending. nobs,mype=',ncd_nobs,mype + append_diag = .false. ! if there are no obs in existing file, then do not try to append + endif + end if + + call nc_diag_init(diag_conv_file, append=append_diag) + + if (.not. append_diag) then ! don't write headers on append - the module will break? + call nc_diag_header("date_time",ianldate ) + call nc_diag_header("Number_of_state_vars", nsdim ) + endif + end subroutine init_netcdf_diag_ + subroutine contents_binary_diag_(odiag) + type(obs_diag),pointer,intent(in):: odiag + + cdiagbuf(ii) = station_id ! station id + + rdiagbuf(1,ii) = ictype(ikx) ! observation type + rdiagbuf(2,ii) = icsubtype(ikx) ! observation subtype + + rdiagbuf(3,ii) = data(ilate,i) ! observation latitude (degrees) + rdiagbuf(4,ii) = data(ilone,i) ! observation longitude (degrees) + + rdiagbuf(5,ii) = data(istnelv,i) ! station elevation (meters) + rdiagbuf(6,ii) = presq ! observation pressure (hPa) + rdiagbuf(7,ii) = fed_height ! observation height (meters) + rdiagbuf(8,ii) = (dtime*r60)-time_offset ! obs time (sec relative to analysis time) + rdiagbuf(9,ii) = data(iqc,i) ! input prepbufr qc or event mark + rdiagbuf(10,ii) = rmiss_single ! setup qc or event mark + rdiagbuf(11,ii) = data(iuse,i) ! read_prepbufr data usage flag + if(muse(i)) then + rdiagbuf(12,ii) = one ! analysis usage flag (1=use, -1=not used) + else + rdiagbuf(12,ii) = -one + endif + + rdiagbuf(13,ii) = rwgt ! nonlinear qc relative weight + rdiagbuf(14,ii) = errinv_input ! prepbufr inverse obs error (dBZ)**-1 + rdiagbuf(15,ii) = errinv_adjst ! read_prepbufr inverse obs error (dBZ)**-1 + rdiagbuf(16,ii) = errinv_final ! final inverse observation error (dBZ)**-1 + rdiagbuf(17,ii) = data(ifedob,i) ! radar reflectivity observation (dBZ) + rdiagbuf(18,ii) = ddiff ! obs-ges (dBZ) + rdiagbuf(19,ii) = data(ifedob,i)-FEDMdiag(i) ! obs-ges w/o bias correction (dBZ) (future slot) + rdiagbuf(20,ii) = dlat8km ! j-index on 8km bufr obs grid + rdiagbuf(21,ii) = dlon8km ! i-index on 8km bufr obs grid + +! print*,'data(ilat,i)=',data(ilat,i),'data(ilon,i)=',data(ilon,i) + + rdiagbuf(22,ii) = FEDMdiag(i) ! dBZ from rain water + + rdiagbuf(23,ii) = T1D ! temperature (sensible, K) + rdiagbuf(24,ii) = RHO ! air density (kg/m**3) + + if (lobsdiagsave) then + write(6,*)'wrong here, stop in setupfed.f90 ' + stop + ioff=nreal + do jj=1,miter + ioff=ioff+1 + if (odiag%muse(jj)) then + rdiagbuf(ioff,ii) = one + else + rdiagbuf(ioff,ii) = -one + endif + enddo + do jj=1,miter+1 + ioff=ioff+1 + rdiagbuf(ioff,ii) = odiag%nldepart(jj) + enddo + do jj=1,miter + ioff=ioff+1 + rdiagbuf(ioff,ii) = odiag%tldepart(jj) + enddo + do jj=1,miter + ioff=ioff+1 + rdiagbuf(ioff,ii) = odiag%obssen(jj) + enddo + endif + + end subroutine contents_binary_diag_ + subroutine contents_netcdf_diag_(odiag) + type(obs_diag),pointer,intent(in):: odiag +! Observation class + character(7),parameter :: obsclass = ' fed' + real(r_kind),dimension(miter) :: obsdiag_iuse + call nc_diag_metadata("Station_ID", station_id ) + call nc_diag_metadata("Observation_Class", obsclass ) + call nc_diag_metadata("Observation_Type", ictype(ikx) ) + call nc_diag_metadata("Observation_Subtype", icsubtype(ikx) ) + call nc_diag_metadata("Latitude", sngl(data(ilate,i)) ) + call nc_diag_metadata("Longitude", sngl(data(ilone,i)) ) + call nc_diag_metadata("Station_Elevation", sngl(data(istnelv,i)) ) + call nc_diag_metadata("Pressure", sngl(presq) ) + call nc_diag_metadata("Height", sngl(fed_height) ) + call nc_diag_metadata("Time", sngl(dtime*r60-time_offset)) + call nc_diag_metadata("Prep_QC_Mark", sngl(data(iqc,i)) ) + call nc_diag_metadata("Prep_Use_Flag", sngl(data(iuse,i)) ) +! call nc_diag_metadata("Nonlinear_QC_Var_Jb", var_jb ! ) + call nc_diag_metadata("Nonlinear_QC_Rel_Wgt", sngl(rwgt) ) + if(muse(i)) then + call nc_diag_metadata("Analysis_Use_Flag", sngl(one) ) + else + call nc_diag_metadata("Analysis_Use_Flag", sngl(-one) ) + endif + + call nc_diag_metadata("Errinv_Input", sngl(errinv_input) ) + call nc_diag_metadata("Errinv_Adjust", sngl(errinv_adjst) ) + call nc_diag_metadata("Errinv_Final", sngl(errinv_final) ) + + call nc_diag_metadata("Observation", sngl(data(ifedob,i)) ) + call nc_diag_metadata("Obs_Minus_Forecast_adjusted", sngl(ddiff) ) + call nc_diag_metadata("Obs_Minus_Forecast_unadjusted", sngl(data(ifedob,i)-FEDMdiag(i)) ) + if (lobsdiagsave) then + do jj=1,miter + if (odiag%muse(jj)) then + obsdiag_iuse(jj) = one + else + obsdiag_iuse(jj) = -one + endif + enddo + + call nc_diag_data2d("ObsDiagSave_iuse", obsdiag_iuse ) + call nc_diag_data2d("ObsDiagSave_nldepart", odiag%nldepart ) + call nc_diag_data2d("ObsDiagSave_tldepart", odiag%tldepart ) + call nc_diag_data2d("ObsDiagSave_obssen" , odiag%obssen ) + endif + + end subroutine contents_netcdf_diag_ + + subroutine final_vars_ + if(allocated(ges_z )) deallocate(ges_z ) + if(allocated(ges_q )) deallocate(ges_q ) +! if(allocated(ges_tv)) deallocate(ges_tv) + if(allocated(ges_ps)) deallocate(ges_ps) + if(allocated(ges_qg)) deallocate(ges_qg) + end subroutine final_vars_ + + subroutine init_qcld(t_cld, qxmin_cld, icat_cld, t_dpnd) + use kinds, only: r_kind,r_single,r_double,i_kind + implicit none + real(r_kind), intent(in ) :: t_cld + real(r_kind), intent(inout) :: qxmin_cld + integer, intent(in ) :: icat_cld + logical, intent(in ) :: t_dpnd +! +! local variables + real :: tr_ll, qrmin_ll, tr_hl, qrmin_hl + real :: ts_ll, qsmin_ll, ts_hl, qsmin_hl + real :: tg_ll, qgmin_ll, tg_hl, qgmin_hl + real :: qr_min, qs_min, qg_min +!------------------------------------------------------ + + qr_min = 5.0E-6_r_kind + qs_min = 5.0E-6_r_kind + qg_min = 5.0E-6_r_kind + tr_ll = 275.65; qrmin_ll = 5.0E-6_r_kind; + tr_hl = 270.65; qrmin_hl = 1.0E-8_r_kind; + ts_ll = 275.65; qsmin_ll = 1.0E-8_r_kind; + ts_hl = 270.65; qsmin_hl = 5.0E-6_r_kind; + tg_ll = 275.65; qgmin_ll = 1.0E-6_r_kind; + tg_hl = 270.65; qgmin_hl = 5.0E-6_r_kind; + + select case (icat_cld) + case (1) + if ( t_dpnd ) then + if (t_cld <= tr_hl) then + qxmin_cld = qrmin_hl + else if (t_cld >= tr_ll) then + qxmin_cld = qrmin_ll + else + qxmin_cld = (qrmin_hl + qrmin_ll) * 0.5 + end if + else + qxmin_cld = qr_min + end if + case default + write(6,*) 'wrong cloud hydrometer category ID',icat_cld + end select + + return + + end subroutine init_qcld + +end subroutine setupfed +end module fed_setup From 0c2780e25a0423e4e6fb38285020162c4e15575d Mon Sep 17 00:00:00 2001 From: David Dowell Date: Tue, 11 Jul 2023 01:33:12 +0000 Subject: [PATCH 2/7] updates for adding FED observations and observation operator to GSI observer --- src/gsi/gsi_fedOper.F90 | 185 ++++++++++++++ src/gsi/m_fedNode.F90 | 265 ++++++++++++++++++++ src/gsi/read_fed.f90 | 527 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 977 insertions(+) create mode 100644 src/gsi/gsi_fedOper.F90 create mode 100644 src/gsi/m_fedNode.F90 create mode 100644 src/gsi/read_fed.f90 diff --git a/src/gsi/gsi_fedOper.F90 b/src/gsi/gsi_fedOper.F90 new file mode 100644 index 0000000000..a306868f77 --- /dev/null +++ b/src/gsi/gsi_fedOper.F90 @@ -0,0 +1,185 @@ +module gsi_fedOper +!$$$ subprogram documentation block +! . . . . +! subprogram: module gsi_fedOper +! prgmmr: j guo +! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.3 +! date: 2018-08-10 +! +! abstract: an obOper extension for fedNode type +! +! program history log: +! 2023-04-10 D. Dowell - moved diag_fed and its description here from +! obsmod. +! +! input argument list: see Fortran 90 style document below +! +! output argument list: see Fortran 90 style document below +! +! attributes: +! language: Fortran 90 and/or above +! machine: +! +!$$$ end subprogram documentation block + +! module interface: + + use gsi_obOper, only: obOper + use m_fedNode , only: fedNode + implicit none + public:: fedOper ! data structure + public:: diag_fed + + type,extends(obOper):: fedOper + contains + procedure,nopass:: mytype + procedure,nopass:: nodeMold + procedure:: setup_ + procedure:: intjo1_ + procedure:: stpjo1_ + end type fedOper + +! def diag_fed- namelist logical to compute/write (=true) FED diag files + logical,save:: diag_fed=.false. + +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + character(len=*),parameter :: myname='gsi_fedOper' + type(fedNode),save,target:: myNodeMold_ + +contains + function mytype(nodetype) + implicit none + character(len=:),allocatable:: mytype + logical,optional, intent(in):: nodetype + mytype="[fedOper]" + if(present(nodetype)) then + if(nodetype) mytype=myNodeMold_%mytype() + endif + end function mytype + + function nodeMold() + !> %nodeMold() returns a mold of its corresponding obsNode + use m_obsNode, only: obsNode + implicit none + class(obsNode),pointer:: nodeMold + nodeMold => myNodeMold_ + end function nodeMold + + subroutine setup_(self, lunin, mype, is, nobs, init_pass, last_pass) + use fed_setup, only: setup + use kinds, only: i_kind + use gsi_obOper, only: len_obstype + use gsi_obOper, only: len_isis + + use m_rhs , only: awork => rhs_awork + use m_rhs , only: bwork => rhs_bwork + use m_rhs , only: iwork => i_fed + + use obsmod , only: write_diag + use jfunc , only: jiter + + use mpeu_util, only: die + + use obsmod, only: dirname, ianldate + + implicit none + class(fedOper ), intent(inout):: self + integer(i_kind), intent(in):: lunin + integer(i_kind), intent(in):: mype + integer(i_kind), intent(in):: is + integer(i_kind), intent(in):: nobs + logical , intent(in):: init_pass ! supporting multi-pass setup() + logical , intent(in):: last_pass ! with incremental backgrounds. + + !---------------------------------------- + character(len=*),parameter:: myname_=myname//"::setup_" + + character(len=len_obstype):: obstype + character(len=len_isis ):: isis + integer(i_kind):: nreal,nchanl,ier,nele + logical:: diagsave + integer(i_kind):: lu_diag + character(128):: diag_file + character(80):: string + + if(nobs == 0) then + + if( (mype == 0) .and. init_pass ) then + write(string,600) jiter +600 format('fed_',i2.2) + diag_file=trim(dirname) // trim(string) + write(6,*) 'write ianldate to ', diag_file + open(newunit=lu_diag,file=trim(diag_file),form='unformatted',status='unknown',position='rewind') + write(lu_diag) ianldate + close(lu_diag) + endif + + return + + endif + + read(lunin,iostat=ier) obstype,isis,nreal,nchanl + if(ier/=0) call die(myname_,'read(obstype,...), iostat =',ier) + nele = nreal+nchanl + + diagsave = write_diag(jiter) .and. diag_fed + + call setup(self%obsLL(:), self%odiagLL(:), & + lunin,mype,bwork,awork(:,iwork),nele,nobs,is,diagsave,init_pass) + + end subroutine setup_ + + subroutine intjo1_(self, ibin, rval,sval, qpred,sbias) +! use intfedmod, only: intjo => intfed + use gsi_bundlemod , only: gsi_bundle + use bias_predictors, only: predictors + use m_obsNode , only: obsNode + use m_obsLList, only: obsLList_headNode + use kinds , only: i_kind, r_quad + implicit none + class(fedOper ),intent(in ):: self + integer(i_kind ),intent(in ):: ibin + type(gsi_bundle),intent(inout):: rval ! (ibin) + type(gsi_bundle),intent(in ):: sval ! (ibin) + real(r_quad ),target,dimension(:),intent(inout):: qpred ! (ibin) + type(predictors),target, intent(in ):: sbias + + !---------------------------------------- + character(len=*),parameter:: myname_=myname//"::intjo1_" + class(obsNode),pointer:: headNode + +! headNode => obsLList_headNode(self%obsLL(ibin)) +! call intjo(headNode, rval,sval) +! headNode => null() + + end subroutine intjo1_ + + subroutine stpjo1_(self, ibin, dval,xval,pbcjo,sges,nstep,dbias,xbias) +! use stpfedmod, only: stpjo => stpfed + use gsi_bundlemod, only: gsi_bundle + use bias_predictors, only: predictors + use m_obsNode , only: obsNode + use m_obsLList, only: obsLList_headNode + use kinds, only: r_quad,r_kind,i_kind + implicit none + class(fedOper ),intent(in):: self + integer(i_kind ),intent(in):: ibin + type(gsi_bundle),intent(in):: dval + type(gsi_bundle),intent(in):: xval + real(r_quad ),dimension(:),intent(inout):: pbcjo ! (1:4) + real(r_kind ),dimension(:),intent(in ):: sges + integer(i_kind),intent(in):: nstep + + type(predictors),target, intent(in):: dbias + type(predictors),target, intent(in):: xbias + + !---------------------------------------- + character(len=*),parameter:: myname_=myname//"::stpjo1_" + class(obsNode),pointer:: headNode + +! headNode => obsLList_headNode(self%obsLL(ibin)) +! call stpjo(headNode,dval,xval,pbcjo(:),sges,nstep) +! headNode => null() + end subroutine stpjo1_ + +end module gsi_fedOper diff --git a/src/gsi/m_fedNode.F90 b/src/gsi/m_fedNode.F90 new file mode 100644 index 0000000000..f42503df19 --- /dev/null +++ b/src/gsi/m_fedNode.F90 @@ -0,0 +1,265 @@ +module m_fedNode +!$$$ subprogram documentation block +! . . . . +! subprogram: module m_fedNode +! prgmmr: YPW +! org: CIMMS +! date: 2019-09-24 +! +! abstract: class-module of obs-type fedNode (GLM flash extent density) +! Modified based on m_tdNode.f90 +! +! program history log: +! 2019-09-24 YPW - added this document block for the initial polymorphic +! implementation. +! +! input argument list: see Fortran 90 style document below +! +! output argument list: see Fortran 90 style document below +! +! attributes: +! language: Fortran 90 and/or above +! machine: +! +!$$$ end subprogram documentation block + +! module interface: + use m_obsdiagNode, only: obs_diag + use m_obsdiagNode, only: obs_diags + use kinds , only: i_kind,r_kind + use mpeu_util, only: assert_,die,perr,warn,tell + use m_obsNode, only: obsNode + + implicit none + private + + public:: fedNode + + type,extends(obsNode):: fedNode + !type(td_ob_type),pointer :: llpoint => NULL() + type(obs_diag), pointer :: diags => NULL() + real(r_kind) :: res ! flash extent density residual + real(r_kind) :: err2 ! flash extent density 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) :: jb ! variational quality control parameter + real(r_kind) :: wij(8) ! horizontal interpolation weights + real(r_kind) :: fedpertb ! random number adding to the obs +! logical :: luse ! flag indicating if ob is used in pen. + integer(i_kind) :: k1 ! level of errtable 1-33 + integer(i_kind) :: kx ! ob type + integer(i_kind) :: ij(8) ! horizontal locations + +! integer(i_kind) :: idv,iob ! device id and obs index for sorting + real (r_kind) :: dlev ! reference to the vertical grid + !real (r_kind) :: elat, elon ! earth lat-lon for redistribution + !real (r_kind) :: dlat, dlon ! earth lat-lon for redistribution + contains + procedure,nopass:: mytype + procedure:: setHop => obsNode_setHop_ + procedure:: xread => obsNode_xread_ + procedure:: xwrite => obsNode_xwrite_ + procedure:: isvalid => obsNode_isvalid_ + procedure:: gettlddp => gettlddp_ + + ! procedure, nopass:: headerRead => obsHeader_read_ + ! procedure, nopass:: headerWrite => obsHeader_write_ + ! procedure:: init => obsNode_init_ + ! procedure:: clean => obsNode_clean_ + end type fedNode + + public:: fedNode_typecast + public:: fedNode_nextcast + interface fedNode_typecast; module procedure typecast_ ; end interface + interface fedNode_nextcast; module procedure nextcast_ ; end interface + + public:: fedNode_appendto + interface fedNode_appendto; module procedure appendto_ ; end interface + + character(len=*),parameter:: MYNAME="m_fedNode" + +!#define CHECKSUM_VERBOSE +!#define DEBUG_TRACE +#include "myassert.H" +#include "mytrace.H" +contains +function typecast_(aNode) result(ptr_) +!-- cast a class(obsNode) to a type(tdNode) + use m_obsNode, only: obsNode + implicit none + type(fedNode),pointer:: ptr_ + class(obsNode),pointer,intent(in):: aNode +! character(len=*),parameter:: myname_=MYNAME//"::typecast_" + ptr_ => null() + if(.not.associated(aNode)) return + select type(aNode) + type is(fedNode) + ptr_ => aNode +! class default +! call die(myname_,'unexpected type, aNode%mytype() =',aNode%mytype()) + end select +return +end function typecast_ + +function nextcast_(aNode) result(ptr_) +!-- cast an obsNode_next(obsNode) to a type(fedNode) + use m_obsNode, only: obsNode,obsNode_next + implicit none + type(fedNode),pointer:: ptr_ + class(obsNode),target,intent(in):: aNode + + class(obsNode),pointer:: inode_ + inode_ => obsNode_next(aNode) + ptr_ => typecast_(inode_) +return +end function nextcast_ + +subroutine appendto_(aNode,oll) +!-- append aNode to linked-list oLL + use m_obsNode , only: obsNode + use m_obsLList, only: obsLList,obsLList_appendNode + implicit none + type(fedNode),pointer,intent(in):: aNode + type(obsLList),intent(inout):: oLL + + class(obsNode),pointer:: inode_ + inode_ => aNode + call obsLList_appendNode(oLL,inode_) + inode_ => null() +end subroutine appendto_ + +! obsNode implementations + +function mytype() + implicit none + character(len=:),allocatable:: mytype + mytype="[fedNode]" +end function mytype + + +subroutine obsNode_xread_(aNode,iunit,istat,diagLookup,skip) + use m_obsdiagNode, only: obsdiagLookup_locate + implicit none + class(fedNode),intent(inout):: aNode + integer(i_kind),intent(in ):: iunit + integer(i_kind),intent( out):: istat + type(obs_diags),intent(in ):: diagLookup + logical,optional,intent(in ):: skip + + character(len=*),parameter:: myname_=MYNAME//'.obsNode_xread_' + logical:: skip_ +_ENTRY_(myname_) + skip_=.false. + if(present(skip)) skip_=skip + + istat=0 + if(skip_) then + read(iunit,iostat=istat) + if(istat/=0) then + call perr(myname_,'skipping read(%(res,err2,...)), iostat =',istat) + _EXIT_(myname_) + return + endif + + else + read(iunit,iostat=istat) aNode%res , & + aNode%err2 , & + aNode%raterr2 , & + aNode%b , & + aNode%pg , & + aNode%jb , & + ! aNode%fedpertb , & + aNode%k1 , & + aNode%kx , & + aNode%dlev , & + aNode%wij , & + aNode%ij + if(istat/=0) then + call perr(myname_,'read(%(res,err2,...), iostat =',istat) + _EXIT_(myname_) + return + endif + + aNode%diags => obsdiagLookup_locate(diagLookup,aNode%idv,aNode%iob,1_i_kind) + if(.not.associated(aNode%diags)) then + call perr(myname_,'obsdiagLookup_locate(), %idv =',aNode%idv) + call perr(myname_,' %iob =',aNode%iob) + call die(myname_) + endif + endif +_EXIT_(myname_) +return +end subroutine obsNode_xread_ + +subroutine obsNode_xwrite_(aNode,junit,jstat) + implicit none + class(fedNode),intent(in):: aNode + integer(i_kind),intent(in ):: junit + integer(i_kind),intent( out):: jstat + + character(len=*),parameter:: myname_=MYNAME//'.obsNode_xwrite_' +_ENTRY_(myname_) + + jstat=0 + write(junit,iostat=jstat) aNode%res , & + aNode%err2 , & + aNode%raterr2 , & + aNode%b , & + aNode%pg , & + aNode%jb , & + ! aNode%fedpertb , & + aNode%k1 , & + aNode%kx , & + aNode%dlev , & + aNode%wij , & + aNode%ij + if(jstat/=0) then + call perr(myname_,'write(%res,err2,...), iostat =',jstat) + _EXIT_(myname_) + return + endif +_EXIT_(myname_) +return +end subroutine obsNode_xwrite_ + +subroutine obsNode_setHop_(aNode) + use m_cvgridLookup, only: cvgridLookup_getiw + implicit none + class(fedNode),intent(inout):: aNode + + character(len=*),parameter:: myname_=MYNAME//'::obsNode_setHop_' +_ENTRY_(myname_) + call cvgridLookup_getiw(aNode%elat,aNode%elon,aNode%dlev,aNode%ij,aNode%wij) +_EXIT_(myname_) +return +end subroutine obsNode_setHop_ + +function obsNode_isvalid_(aNode) result(isvalid_) + implicit none + logical:: isvalid_ + class(fedNode),intent(in):: aNode + + character(len=*),parameter:: myname_=MYNAME//'::obsNode_isvalid_' +_ENTRY_(myname_) + isvalid_=associated(aNode%diags) +_EXIT_(myname_) +return +end function obsNode_isvalid_ + +pure subroutine gettlddp_(aNode,jiter,tlddp,nob) + use kinds, only: r_kind + implicit none + class(fedNode), intent(in):: aNode + integer(kind=i_kind),intent(in):: jiter + real(kind=r_kind),intent(inout):: tlddp + integer(kind=i_kind),optional,intent(inout):: nob + + tlddp = tlddp + aNode%diags%tldepart(jiter)*aNode%diags%tldepart(jiter) + if(present(nob)) nob=nob+1 +return +end subroutine gettlddp_ + +end module m_fedNode diff --git a/src/gsi/read_fed.f90 b/src/gsi/read_fed.f90 new file mode 100644 index 0000000000..79c6e2a726 --- /dev/null +++ b/src/gsi/read_fed.f90 @@ -0,0 +1,527 @@ +subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) +!$$$ subprogram documentation block +! . . . . +! ABSTRACT: +! This routine reads in netcdf or prepbufr flash-extent density (FED) data. +! +! PROGRAM HISTORY LOG: +! 2018-07-25 Rong Kong (CAPS/OU) - modified based on read_radarref_mosaic.f90 +! 2019-09-20 Yaping Wang (CIMMS/OU) +! 2021-07-01 David Dowell (DCD; NOAA GSL) - added maximum flashes/min for observed FED +! +! input argument list: +! infile - unit from which to read observation information file +! obstype - observation type to process +! lunout - unit to which to write data for further processing +! twind - input group time window (hours) +! sis - observation variable name +! +! output argument list: +! nread - number of type "obstype" observations read +! ndata - number of type "obstype" observations retained for further processing +! nobs - array of observations on each subdomain for each processor +! +! ATTRIBUTES: +! LANGUAGE: FORTRAN 90 +! +!$$$ +! +!_____________________________________________________________________ +! + use kinds, only: r_kind,r_double,i_kind + use constants, only: zero,one,rad2deg,deg2rad + use convinfo, only: nconvtype,ctwind,cgross,cermax,cermin,cvar_b,cvar_pg, & + ncmiter,ncgroup,ncnumgrp,icuse,ictype,icsubtype,ioctype + use gsi_4dvar, only: l4dvar,l4densvar,winlen + use gridmod, only: tll2xy + use mod_wrfmass_to_a, only: wrfmass_obs_to_a8 + use mpimod, only: npe + use obsmod, only: perturb_obs,iadatemn,time_window + + use netcdf + implicit none + + include 'netcdf.inc' +! + character(len=*), intent(in) :: infile,obstype + integer(i_kind), intent(in) :: lunout + integer(i_kind), intent(inout) :: nread,ndata + integer(i_kind), intent(inout) :: nodata + integer(i_kind), dimension(npe) ,intent(inout) :: nobs + real(r_kind), intent(in ) :: twind + character(len=*), intent(in) :: sis + +! Declare local parameters + real(r_kind),parameter:: r90 = 90.0_r_kind + real(r_kind),parameter:: r360 = 360.0_r_kind + real(r_kind),parameter:: oe_fed = 1.0_r_kind + real(r_kind),parameter:: fed_lowbnd = 0.1_r_kind ! use fed == fed_lowbnd + real(r_kind),parameter:: fed_lowbnd2 = 0.1_r_kind ! use fed >= fed_lowbnd2 +! real(r_kind),parameter:: fed_highbnd = 18.0_r_kind ! 18 flashes/min from Sebok and Back (2021, unpublished) + real(r_kind),parameter:: fed_highbnd = 8.0_r_kind ! 8 flashes/min from Back (2023) for regional FV3 tests + +! +! For fed observations +! + integer(i_kind) nreal,nchanl + + integer(i_kind) ifn,i + + real(r_kind) :: maxfed + integer(i_kind) :: ilon,ilat + + logical :: fedobs, fedob + real(r_kind),allocatable,dimension(:,:):: cdata_out + real(r_kind) :: federr, thiserr + real(r_kind) :: hgt_fed(1) + data hgt_fed / 6500.0 / + + real(r_kind) :: i_maxloc,j_maxloc,k_maxloc + integer(i_kind) :: kint_maxloc + real(r_kind) :: fed_max + integer(i_kind) :: ndata2 + integer(i_kind) :: ppp + +! +! for read in bufr +! + real(r_kind) :: hdr(5),obs(1,3) + character(80):: hdrstr='SID XOB YOB DHR TYP' + character(80):: obsstr='FED' + + character(8) subset + character(8) station_id + real(r_double) :: rstation_id + equivalence(rstation_id,station_id) + integer(i_kind) :: lunin,idate + integer(i_kind) :: ireadmg,ireadsb + + integer(i_kind) :: maxlvl + integer(i_kind) :: numlvl,numfed,numobsa,nmsgmax,maxobs + integer(i_kind) :: k,iret + integer(i_kind) :: nmsg,ntb + + real(r_kind),allocatable,dimension(:,:) :: fed3d_column ! 3D fed in column + real(r_kind),allocatable,dimension(:) :: utime ! time + + integer(i_kind) :: ikx + real(r_kind) :: timeo,t4dv + + character*128 :: myname='read_fed' + + real(r_kind) :: dlat, dlon ! rotated corrdinate + real(r_kind) :: dlat_earth, dlon_earth ! in unit of degree + real(r_kind) :: rlat00, rlon00 ! in unit of rad + + logical :: l_psot_fed + logical :: l_latlon_fedobs + logical :: outside + integer :: unit_table + +! for read netcdf + integer(i_kind) :: idate5(5), sec70,mins_an,mins_ob + integer(i_kind) :: varID, ncdfID, status + character(4) :: idate5s(5) + real(r_kind) :: timeb,twindm,rmins_an,rmins_ob + + + unit_table = 23 +!********************************************************************** +! +! END OF DECLARATIONS....start of program +! + write(6,*) "r_kind=",r_kind + l_psot_fed = .FALSE. + l_latlon_fedobs = .TRUE. + + fedob = obstype == 'fed' + if(fedob) then + nreal=25 + else + write(6,*) ' illegal obs type in read_fed : obstype=',obstype + call stop2(94) + end if + if(perturb_obs .and. fedob)nreal=nreal+1 + write(6,*)'read_fed: nreal=',nreal + + fedobs = .false. + ikx=0 + do i=1,nconvtype + if(trim(obstype) == trim(ioctype(i)) .and. abs(icuse(i))== 1) then + fedobs=.true. + ikx=i + federr = oe_fed ! Obs error (flashes per minute) + thiserr = federr + exit ! Exit loop when finished with initial convinfo fields + else if (i == nconvtype ) then + write(6,*) 'read_fed: Obs Type for fed is not in CONVINFO !' + write(6,*) 'read_fed: PLEASE modify the CONVINFO file !' + write(6,*) 'read_fed: abort read_fed !' + return + endif + end do + write(6,'(1x,A,A30,I4,A15,F7.3,A7)') & + trim(myname),': fed in convinfo-->ikx=',ikx,' fed ob err:',thiserr," (fed)" + + nread=0 + ndata=0 + nchanl=0 + ifn = 15 + + if(fedobs) then + maxlvl= 1 ! fed only has one level + + if(trim(infile) .eq. "fedbufr") then ! prebufr or netcdf format + !! get message and subset counts + ! nmsgmax and maxobs are read in from BUFR data file, not pre-set. + call getcount_bufr(infile,nmsgmax,maxobs) + write(6,*)'read_fed: nmsgmax=',nmsgmax,' maxobs=',maxobs + +! read in fed obs in bufr code format + lunin = 10 + allocate(fed3d_column(maxlvl+2+2,maxobs)) + + open ( unit = lunin, file = trim(infile),form='unformatted',err=200) + call openbf ( lunin, 'IN', lunin ) + open(unit_table,file='prepobs_kr.bufrtable') !temporily dump the bufr table, which is already saved in file + call dxdump(lunin,unit_table) + call datelen ( 10 ) + + nmsg=0 + ntb = 0 + + ndata =0 + ppp = 0 + msg_report: do while (ireadmg(lunin,subset,idate) == 0) + nmsg=nmsg+1 + if (nmsg>nmsgmax) then + write(6,*)'read_fed: messages exceed maximum ',nmsgmax + call stop2(50) + endif + loop_report: do while (ireadsb(lunin) == 0) + ntb = ntb+1 + if (ntb>maxobs) then + write(6,*)'read_fed: reports exceed maximum ',maxobs + call stop2(50) + endif + + ! Extract type, date, and location information from BUFR file + call ufbint(lunin,hdr,5,1,iret,hdrstr) + if(hdr(3) .gt. 90 ) write(6,*) "Inside read_fed.f90, hdr(2)=",hdr(2),"hdr(3)=",hdr(3) + if ( l_latlon_fedobs ) then + if(abs(hdr(3))>r90 .or. abs(hdr(2))>r360) cycle loop_report + if(hdr(2)== r360)hdr(2)=hdr(2)-r360 + if(hdr(2) < zero)hdr(2)=hdr(2)+r360 + end if + +! check time window in subset + if (l4dvar.or.l4densvar) then + t4dv=hdr(4) + if (t4dvwinlen) then + write(6,*)'read_fed: time outside window ',& + t4dv,' skip this report' + cycle loop_report + endif + else + timeo=hdr(4) + if (abs(timeo)>ctwind(ikx) .or. abs(timeo) > twind) then + write(6,*)'read_fed: time outside window ',& + timeo,' skip this report' + cycle loop_report + endif + endif +! read in observations + call ufbint(lunin,obs,1,3,iret,obsstr) !Single level bufr data, Rong Kong + if(obs(1,1) .gt. 5 ) write(6,*) "Inside read_fed.f90, obs(1,1)=",obs(1,1) + numlvl=min(iret,maxlvl) + if (numlvl .ne. maxlvl) then + write(6,*)' read_fed: numlvl is not equalt to maxlvl:',numlvl,maxlvl + end if + if(hdr(3) .gt. 90) write(6,*) "hdr(3)=",hdr(3) + if ( l_latlon_fedobs ) then + if(hdr(2)>= r360)hdr(2)=hdr(2)-r360 + if(hdr(2) < zero)hdr(2)=hdr(2)+r360 + fed3d_column(1,ntb)=hdr(2) ! observation location, earth lon + fed3d_column(2,ntb)=hdr(3) ! observation location, earth lat +! write(6,*) "Inside read_fed.f90, fed3d_column(1,ntb)=",fed3d_column(1,ntb),"fed3d_column(2,ntb)=",fed3d_column(2,ntb) + else + fed3d_column(1,ntb)=hdr(2)*10.0_r_kind ! observation location, grid index i + fed3d_column(2,ntb)=hdr(3)*10.0_r_kind ! observation location, grid index j + end if + + if (l_psot_fed .and. .NOT. l_latlon_fedobs ) then + do k=1,numlvl + if (NINT(fed3d_column(1,ntb)) .eq. 175 .and. NINT(fed3d_column(2,ntb)) .eq. 105 .and. & + NINT(hgt_fed(k)) .ge. 100 ) then + write(6,*) 'read_fed: single point/column obs run on grid: 175 105' + write(6,*) 'read_fed: found the pseudo single(column) fed obs:',fed3d_column(1:2,ntb),hgt_fed(k) + else + obs(1,1) = -999.0 + end if + end do + end if + + fed3d_column(3,ntb)=obs(1,1) + fed3d_column(4,ntb)=obs(1,2) + fed3d_column(5,ntb)=obs(1,3) + if (obs(1,1) == fed_lowbnd .or. obs(1,1) >= fed_lowbnd2 ) then + if (obs(1,1) == 0.0) then + ppp = ppp + 1 + endif + ndata = ndata + 1 + endif + + enddo loop_report + enddo msg_report + + write(6,*)'read_fed: messages/reports = ',nmsg,'/',ntb + print*,'number of Z that is less than 0 is ppp = ', ppp + numfed=ntb + +! - Finished reading fed observations from BUFR format data file +! + call closbf(lunin) + close(lunin) + + else ! NETCDF format +!!!! Start reading fed observations from NETCDF format data file + ! CHECK IF DATA FILE EXISTS + + ! OPEN NETCDF FILE + status = nf90_open(TRIM(infile), NF90_NOWRITE, ncdfID) + print*, '*** OPENING GOES FED OBS NETCDF FILE: ', infile, status + + + !------------------------ + ! Get date information + !------------------------- + ! status = nf90_get_att( ncdfID, nf90_global, 'year', idate5s(1) ) + ! print*, 'year ',status + ! status = nf90_get_att( ncdfID, nf90_global, 'month', idate5s(2) ) + ! status = nf90_get_att( ncdfID, nf90_global, 'day', idate5s(3) ) + ! status = nf90_get_att( ncdfID, nf90_global, 'hour', idate5s(4) ) + ! status = nf90_get_att( ncdfID, nf90_global, 'minute', idate5s(5) ) + ! read(idate5s(:) , *) idate5(:) + ! print*, idate5 + + !------------------------ + ! Get Dimension Info (1-D) + !------------------------- + status = nf90_inq_varid( ncdfID, 'numobs', varID ) + status = nf90_get_var( ncdfID, varID, maxobs ) + + !------------------------ + ! Allocate data arrays + !------------------------- + ALLOCATE( fed3d_column( 5, maxobs ) ) + ALLOCATE( utime( 1 ) ) ! seconds since from 2000-01-01 12:00 + + !------------------------ + ! Get useful data arrays + !------------------------- + ! LON + status = nf90_inq_varid( ncdfID, 'lon', varID ) + status = nf90_get_var( ncdfID, varID, fed3d_column(1, :) ) + ! LAT + status = nf90_inq_varid( ncdfID, 'lat', varID ) + status = nf90_get_var( ncdfID, varID, fed3d_column(2, :) ) + ! FED value + status = nf90_inq_varid( ncdfID, 'value', varID ) + status = nf90_get_var( ncdfID, varID, fed3d_column(3, :) ) + ! TIME + status = nf90_inq_varid( ncdfID, 'time', varID ) + status = nf90_get_var( ncdfID, varID, utime ) + + ! CLOSE NETCDF FILE + status = nf90_close( ncdfID ) + + + !-Obtain analysis time in minutes since reference date + sec70 = 694267200.0 ! seconds since from 1978-01-01 00:00 to 2000-01-01 12:00 + ! because the official GOES prescribed epoch time for GLM data is 2000-01-01 12:00:00 + + call w3fs21(iadatemn,mins_an) !mins_an -integer number of mins snce 01/01/1978 + rmins_an=mins_an !convert to real number + + ! SINCE ALL OBS WILL HAVE THE SAME TIME, CHECK TIME HERE: + rmins_ob = ( utime(1) + sec70 )/60 !Convert to Minutes from seconds + twindm = twind*60. !Convert to Minutes from hours + timeb = rmins_ob-rmins_an + + if(abs(timeb) > abs(twindm)) then + print*, 'WARNING: ALL FED OBSERVATIONS OUTSIDE ASSIMILATION TIME WINDOW: ', timeb, twindm + ! goto 314 + endif + numfed = maxobs + do i=1,numfed + if (fed3d_column( 3, i ) >= fed_lowbnd2 .or. fed3d_column( 3, i ) == fed_lowbnd ) then + ndata = ndata + 1 + end if + end do + end if ! end if prebufr or netcdf format + + write(6,*)'read_fed: total no. of obs = ',ndata + nread=ndata + nodata=ndata +!!! - Finished reading fed observations from NETCDF format data file + + + + allocate(cdata_out(nreal,ndata)) +! +! + DO i=1,numfed + DO k=1,maxlvl + +! DCD 1 July 2021 + if (fed3d_column(k+2,i) .gt. fed_highbnd) fed3d_column(k+2,i) = fed_highbnd + + end do + end do + + write(6,*) ' ------- check max and min value of OBS: bufr fed -------' + write(6,*) ' level maxval(fed) minval(fed)' + DO k=1,maxlvl + write(6,*) k,maxval(fed3d_column(k+2,:)),minval(fed3d_column(k+2,:)) + end do + + + i_maxloc=-1.0 + j_maxloc=-1.0 + k_maxloc=-1.0 + kint_maxloc=-1 + fed_max=-999.99 + ndata2=0 + DO i=1,numfed + DO k=1,maxlvl + if( fed3d_column(k+2,i) >= fed_lowbnd2 .or. fed3d_column(k+2,i) == fed_lowbnd) then !Rong Kong + dlon_earth = fed3d_column(1,i) ! longitude (degrees) of observation + ! ilone=18 ! index of longitude (degrees) + dlat_earth = fed3d_column(2,i) ! latitude (degrees) of observation + ! ilate=19 ! index of latitude (degrees) + !-Check format of longitude and correct if necessary + if(dlon_earth>=r360) dlon_earth=dlon_earth-r360 + if(dlon_earth 0 ) then + call count_obs(ndata,nreal,ilat,ilon,cdata_out,nobs) + write(lunout) obstype,sis,nreal,nchanl,ilat,ilon,ndata + write(lunout) ((cdata_out(k,i),k=1,nreal),i=1,ndata) + ! print*,'cdata_out',cdata_out +! endif + + deallocate(cdata_out) + if (allocated(fed3d_column)) deallocate(fed3d_column) + + write(6,'(1x,A,F12.5,1x,A,3(1x,F8.3),1x,I4)') & + 'read_fed: max fed =',fed_max, '@ i j k =', & + i_maxloc,j_maxloc,k_maxloc,kint_maxloc + + end if +! close(lunout) ! ???? + return + +200 continue + write(6,*) 'read_fed, Warning : cannot find or open bufr fed data file: ', trim(infile) + +314 continue +print* ,'FINISHED WITH READ_FED' +end subroutine read_fed +! +! From a54b475030431523a443a5028c523d3784615cf5 Mon Sep 17 00:00:00 2001 From: David Dowell Date: Tue, 11 Jul 2023 01:36:05 +0000 Subject: [PATCH 3/7] updates to add FED observations and observation operator to GSI observer --- src/gsi/gsi_files.cmake | 4 ++ src/gsi/gsi_obOperTypeManager.F90 | 7 +++ src/gsi/gsimod.F90 | 7 ++- src/gsi/intjo.f90 | 4 +- src/gsi/m_obsNodeTypeManager.F90 | 7 +++ src/gsi/m_rhs.F90 | 2 + src/gsi/obsmod.F90 | 15 +++++-- src/gsi/read_obs.F90 | 14 +++++- src/gsi/setuprhsall.f90 | 3 +- src/gsi/statsconv.f90 | 72 +++++++++++++++++++++++++++++-- 10 files changed, 122 insertions(+), 13 deletions(-) diff --git a/src/gsi/gsi_files.cmake b/src/gsi/gsi_files.cmake index b98cd2d0da..b514e11c1e 100644 --- a/src/gsi/gsi_files.cmake +++ b/src/gsi/gsi_files.cmake @@ -217,6 +217,7 @@ gsi_colvkOper.F90 gsi_dbzOper.F90 gsi_dwOper.F90 gsi_enscouplermod.f90 +gsi_fedOper.F90 gsi_gpsbendOper.F90 gsi_gpsrefOper.F90 gsi_gustOper.F90 @@ -338,6 +339,7 @@ m_distance.f90 m_dtime.F90 m_dwNode.F90 m_extOzone.F90 +m_fedNode.F90 m_find.f90 m_gpsNode.F90 m_gpsrhs.F90 @@ -478,6 +480,7 @@ read_cris.f90 read_dbz_nc.f90 read_dbz_netcdf.f90 read_diag.f90 +read_fed.f90 read_files.f90 read_fl_hdob.f90 read_gfs_ozone_for_regional.f90 @@ -532,6 +535,7 @@ setupco.f90 setupdbz.f90 setupdbz_lib.f90 setupdw.f90 +setupfed.f90 setupgust.f90 setuphowv.f90 setuplag.f90 diff --git a/src/gsi/gsi_obOperTypeManager.F90 b/src/gsi/gsi_obOperTypeManager.F90 index ea306953c4..f7aef2a026 100644 --- a/src/gsi/gsi_obOperTypeManager.F90 +++ b/src/gsi/gsi_obOperTypeManager.F90 @@ -66,6 +66,7 @@ module gsi_obOperTypeManager use gsi_lightOper , only: lightOper use gsi_dbzOper , only: dbzOper + use gsi_fedOper , only: fedOper use gsi_cldtotOper , only: cldtotOper use kinds , only: i_kind @@ -136,6 +137,7 @@ module gsi_obOperTypeManager public:: iobOper_lwcp public:: iobOper_light public:: iobOper_dbz + public:: iobOper_fed public:: iobOper_cldtot enum, bind(C) @@ -181,6 +183,7 @@ module gsi_obOperTypeManager enumerator:: iobOper_lwcp enumerator:: iobOper_light enumerator:: iobOper_dbz + enumerator:: iobOper_fed enumerator:: iobOper_cldtot enumerator:: iobOper_extra_ @@ -242,6 +245,7 @@ module gsi_obOperTypeManager type( lwcpOper), target, save:: lwcpOper_mold type( lightOper), target, save:: lightOper_mold type( dbzOper), target, save:: dbzOper_mold + type( fedOper), target, save:: fedOper_mold type( cldtotOper), target, save:: cldtotOper_mold contains @@ -388,6 +392,7 @@ function dtype2index_(dtype) result(index_) case("goes_glm" ); index_= iobOper_light case("dbz" ,"[dbzoper]" ); index_= iobOper_dbz + case("fed" ,"[fedoper]" ); index_= iobOper_fed case("cldtot" ,"[cldtotoper]" ); index_= iobOper_cldtot case("mta_cld" ); index_= iobOper_cldtot @@ -485,6 +490,7 @@ function index2vmold_(iobOper) result(vmold_) case(iobOper_lwcp ); vmold_ => lwcpOper_mold case(iobOper_light ); vmold_ => lightOper_mold case(iobOper_dbz ); vmold_ => dbzOper_mold + case(iobOper_fed ); vmold_ => fedOper_mold case(iobOper_cldtot ); vmold_ => cldtotOper_mold case( obOper_undef ); vmold_ => null() @@ -600,6 +606,7 @@ subroutine cobstype_config_() cobstype(iobOper_lwcp ) ="lwcp " ! lwcp_ob_type cobstype(iobOper_light ) ="light " ! light_ob_type cobstype(iobOper_dbz ) ="dbz " ! dbz_ob_type + cobstype(iobOper_fed ) ="fed " ! fed_ob_type cobstype(iobOper_cldtot ) ="cldtot " ! using q_ob_type cobstype_configured_=.true. diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index cf885c2b64..de19c85fab 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -21,6 +21,7 @@ module gsimod lread_obs_save,lread_obs_skip,time_window_rad,tcp_posmatch,tcp_box, & neutral_stability_windfact_2dvar,use_similarity_2dvar,ta2tb use gsi_dbzOper, only: diag_radardbz + use gsi_fedOper, only: diag_fed use obsmod, only: doradaroneob,oneoblat,oneoblon,oneobheight,oneobvalue,oneobddiff,oneobradid,& radar_no_thinning,ens_hx_dbz_cut,static_gsi_nopcp_dbz,rmesh_dbz,& @@ -560,6 +561,7 @@ module gsimod ! diag_co - logical to turn off or on the diagnostic carbon monoxide file (true=on) ! diag_light - logical to turn off or on the diagnostic lightning file (true=on) ! diag_radardbz - logical to turn off or on the diagnostic radar reflectivity file (true=on) +! diag_fed - logical to turn off or on the diagnostic flash extent density file (true=on) ! write_diag - logical to write out diagnostic files on outer iteration ! lobsdiagsave - write out additional observation diagnostics ! ltlint - linearize inner loop @@ -738,8 +740,8 @@ module gsimod min_offset,pseudo_q2,& iout_iter,npredp,retrieval,& tzr_qc,tzr_bufrsave,& - diag_rad,diag_pcp,diag_conv,diag_ozone,diag_aero,diag_co,diag_light,diag_radardbz,iguess, & - write_diag,reduce_diag, & + diag_rad,diag_pcp,diag_conv,diag_ozone,diag_aero,diag_co,diag_light,diag_radardbz,diag_fed, & + iguess,write_diag,reduce_diag, & oneobtest,sfcmodel,dtbduv_on,ifact10,l_foto,offtime_data,& use_pbl,use_compress,nsig_ext,gpstop,commgpstop, commgpserrinf, & perturb_obs,perturb_fact,oberror_tune,preserve_restart_date, & @@ -1977,6 +1979,7 @@ subroutine gsimain_initialize diag_pcp=.false. diag_light=.false. diag_radardbz=.false. + diag_fed=.false. use_limit = 0 end if if(reduce_diag) use_limit = 0 diff --git a/src/gsi/intjo.f90 b/src/gsi/intjo.f90 index 91b811147e..a68355471b 100644 --- a/src/gsi/intjo.f90 +++ b/src/gsi/intjo.f90 @@ -31,7 +31,7 @@ module intjomod use gsi_obOperTypeManager, only: & iobOper_t, iobOper_pw, iobOper_q, & iobOper_cldtot, iobOper_w, iobOper_dw, & - iobOper_rw, iobOper_dbz, & + iobOper_rw, iobOper_dbz, iobOper_fed, & iobOper_spd, iobOper_oz, iobOper_o3l, iobOper_colvk, & iobOper_pm2_5, iobOper_pm10, iobOper_ps, iobOper_tcp, iobOper_sst, & iobOper_gpsbend, iobOper_gpsref, & @@ -60,7 +60,7 @@ module intjomod integer(i_kind),parameter,dimension(obOper_count):: ix_obtype = (/ & iobOper_t, iobOper_pw, iobOper_q, & iobOper_cldtot, iobOper_w, iobOper_dw, & - iobOper_rw, iobOper_dbz, & + iobOper_rw, iobOper_dbz, iobOper_fed, & iobOper_spd, iobOper_oz, iobOper_o3l, iobOper_colvk, & iobOper_pm2_5, iobOper_pm10, iobOper_ps, iobOper_tcp, iobOper_sst, & iobOper_gpsbend, iobOper_gpsref, & diff --git a/src/gsi/m_obsNodeTypeManager.F90 b/src/gsi/m_obsNodeTypeManager.F90 index b5ecc6e1ba..43b42e4bf2 100644 --- a/src/gsi/m_obsNodeTypeManager.F90 +++ b/src/gsi/m_obsNodeTypeManager.F90 @@ -70,6 +70,7 @@ module m_obsNodeTypeManager use m_lightNode, only: lightNode use m_dbzNode , only: dbzNode + use m_fedNode, only: fedNode use kinds, only: i_kind use m_obsNode, only: obsNode @@ -124,6 +125,7 @@ module m_obsNodeTypeManager public:: iobsNode_light public:: iobsNode_dbz + public:: iobsNode_fed public :: obsNode_typeMold public :: obsNode_typeIndex @@ -179,6 +181,7 @@ module m_obsNodeTypeManager type( lwcpNode), target, save:: lwcp_mold type( lightNode), target, save:: light_mold type( dbzNode), target, save:: dbz_mold + type( fedNode), target, save:: fed_mold !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ character(len=*),parameter :: myname='m_obsNodeTypeManager' @@ -245,6 +248,7 @@ module m_obsNodeTypeManager enumerator:: iobsNode_lwcp enumerator:: iobsNode_light enumerator:: iobsNode_dbz + enumerator:: iobsNode_fed enumerator:: iobsNode_extra_ end enum @@ -314,6 +318,7 @@ function vname2index_(vname) result(index_) case("light","[lightnode]"); index_ = iobsNode_light case("dbz" , "[dbznode]"); index_ = iobsNode_dbz + case("fed" , "[fednode]"); index_ = iobsNode_fed end select end function vname2index_ @@ -377,6 +382,7 @@ function vmold2index_select_(mold) result(index_) type is(lightNode); index_ = iobsNode_light type is( dbzNode); index_ = iobsNode_dbz + type is( fedNode); index_ = iobsNode_fed end select end function vmold2index_select_ @@ -434,6 +440,7 @@ function index2vmold_(i_obType) result(obsmold_) case(iobsNode_light); obsmold_ => light_mold case(iobsNode_dbz); obsmold_ => dbz_mold + case(iobsNode_fed); obsmold_ => fed_mold end select end function index2vmold_ diff --git a/src/gsi/m_rhs.F90 b/src/gsi/m_rhs.F90 index baee074688..aea417fe27 100644 --- a/src/gsi/m_rhs.F90 +++ b/src/gsi/m_rhs.F90 @@ -80,6 +80,7 @@ module m_rhs public:: i_lwcp public:: i_light public:: i_dbz + public:: i_fed public:: i_cldtot public:: awork_size @@ -146,6 +147,7 @@ module m_rhs enumerator:: i_lwcp enumerator:: i_light enumerator:: i_dbz + enumerator:: i_fed enumerator:: i_cldtot enumerator:: i_outbound diff --git a/src/gsi/obsmod.F90 b/src/gsi/obsmod.F90 index 3066cdb5ca..a059586e67 100644 --- a/src/gsi/obsmod.F90 +++ b/src/gsi/obsmod.F90 @@ -160,6 +160,7 @@ module obsmod ! 2021-11-16 Zhao - add option l_obsprvdiag (if true) to trigger the output of ! observation provider and sub-provider information into ! obsdiags files (used for AutoObsQC) +! 2023-07-10 Y. Wang, D. Dowell - add variables for flash extent density ! ! Subroutines Included: ! sub init_obsmod_dflts - initialize obs related variables to default values @@ -186,6 +187,7 @@ module obsmod ! def write_diag - namelist logical array to compute/write (=true) diag files ! def diag_radardbz- namelist logical to compute/write (=true) radar ! reflectiivty diag files +! def diag_fed - namelist logical to compute/write (=true) flash extent density diag files ! def reduce_diag - namelist logical to produce reduced radiance diagnostic files ! def use_limit - parameter set equal to -1 if diag files produced or 0 if not diag files or reduce_diag ! def obs_setup - prefix for files passing pe relative obs data to setup routines @@ -434,6 +436,7 @@ module obsmod public :: ran01dom,dval_use public :: iout_pcp,iout_rad,iadate,iadatemn,write_diag,reduce_diag,oberrflg,bflag,ndat,dthin,dmesh,l_do_adjoint public :: diag_radardbz + public :: diag_fed public :: lsaveobsens public :: iout_cldch, mype_cldch public :: nprof_gps,time_offset,ianldate,tcp_box @@ -483,7 +486,9 @@ module obsmod public :: iout_dbz, mype_dbz ! --- DBZ DA --- - + + public :: iout_fed, mype_fed + public :: obsmod_init_instr_table public :: obsmod_final_instr_table public :: nobs_sub @@ -583,12 +588,12 @@ module obsmod integer(i_kind) iout_co,iout_gust,iout_vis,iout_pblh,iout_tcamt,iout_lcbas integer(i_kind) iout_cldch integer(i_kind) iout_wspd10m,iout_td2m,iout_mxtm,iout_mitm,iout_pmsl,iout_howv - integer(i_kind) iout_uwnd10m,iout_vwnd10m + integer(i_kind) iout_uwnd10m,iout_vwnd10m,iout_fed integer(i_kind) mype_t,mype_q,mype_uv,mype_ps,mype_pw, & mype_rw,mype_dw,mype_gps,mype_sst, & mype_tcp,mype_lag,mype_co,mype_gust,mype_vis,mype_pblh, & mype_wspd10m,mype_td2m,mype_mxtm,mype_mitm,mype_pmsl,mype_howv,& - mype_uwnd10m,mype_vwnd10m, mype_tcamt,mype_lcbas, mype_dbz + mype_uwnd10m,mype_vwnd10m, mype_tcamt,mype_lcbas, mype_dbz, mype_fed integer(i_kind) mype_cldch integer(i_kind) iout_swcp, iout_lwcp integer(i_kind) mype_swcp, mype_lwcp @@ -638,6 +643,7 @@ module obsmod logical lobserver,l_do_adjoint, lobsdiag_forenkf logical,dimension(0:50):: write_diag logical diag_radardbz + logical diag_fed logical reduce_diag logical offtime_data logical hilbert_curve @@ -789,6 +795,7 @@ subroutine init_obsmod_dflts end do write_diag(1)=.true. diag_radardbz = .false. + diag_fed = .false. reduce_diag = .false. use_limit = -1 lobsdiagsave=.false. @@ -853,6 +860,7 @@ subroutine init_obsmod_dflts iout_lwcp=236 ! liquid-water content path iout_light=237 ! lightning iout_dbz=238 ! radar reflectivity + iout_fed=239 ! flash extent density mype_ps = npe-1 ! surface pressure mype_t = max(0,npe-2) ! temperature @@ -887,6 +895,7 @@ subroutine init_obsmod_dflts mype_lwcp=max(0,npe-31) ! liquid-water content path mype_light=max(0,npe-32)! GOES/GLM lightning mype_dbz=max(0,npe-33) ! radar reflectivity + mype_fed= max(0,npe-34) ! flash extent density ! Initialize arrays used in namelist obs_input diff --git a/src/gsi/read_obs.F90 b/src/gsi/read_obs.F90 index 9017c498c2..ee0209639b 100644 --- a/src/gsi/read_obs.F90 +++ b/src/gsi/read_obs.F90 @@ -192,6 +192,7 @@ subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread) if ( .not. l_use_dbz_directDA) then if(trim(dtype) == 'dbz' )return end if + if(trim(dtype) == 'fed' )return ! Use routine as usual @@ -910,7 +911,8 @@ subroutine read_obs(ndata,mype) obstype == 'mitm' .or. obstype=='pmsl' .or. & obstype == 'howv' .or. obstype=='tcamt' .or. & obstype=='lcbas' .or. obstype=='cldch' .or. obstype == 'larcglb' .or. & - obstype=='uwnd10m' .or. obstype=='vwnd10m' .or. obstype=='dbz' ) then + obstype=='uwnd10m' .or. obstype=='vwnd10m' .or. obstype=='dbz' .or. & + obstype=='fed') then ditype(i) = 'conv' else if (obstype == 'swcp' .or. obstype == 'lwcp') then ditype(i) = 'wcp' @@ -1295,6 +1297,10 @@ subroutine read_obs(ndata,mype) use_hgtl_full=.true. if(belong(i))use_hgtl_full_proc=.true. end if + if(obstype == 'fed')then + use_hgtl_full=.true. + if(belong(i))use_hgtl_full_proc=.true. + end if if(obstype == 'sst')then if(belong(i))use_sfc=.true. endif @@ -1632,6 +1638,12 @@ subroutine read_obs(ndata,mype) endif end if +! Process flash extent density + else if (obstype == 'fed' ) then + print *, "calling read_fed" + call read_fed(nread,npuse,nouse,infile,obstype,lunout,twind,sis,nobs_sub1(1,i)) + string='READ_FED' + ! Process lagrangian data else if (obstype == 'lag') then call read_lag(nread,npuse,nouse,infile,lunout,obstype,& diff --git a/src/gsi/setuprhsall.f90 b/src/gsi/setuprhsall.f90 index 3efcb69859..8075956431 100644 --- a/src/gsi/setuprhsall.f90 +++ b/src/gsi/setuprhsall.f90 @@ -168,6 +168,7 @@ subroutine setuprhsall(ndata,mype,init_pass,last_pass) i_gust,i_vis,i_pblh,i_wspd10m,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv, & i_tcamt,i_lcbas,i_cldch,i_uwnd10m,i_vwnd10m,i_swcp,i_lwcp use m_rhs, only: i_dbz + use m_rhs, only: i_fed use m_rhs, only: i_light use m_gpsStats, only: gpsStats_genstats ! was genstats_gps() @@ -625,7 +626,7 @@ subroutine setuprhsall(ndata,mype,init_pass,last_pass) call statsconv(mype,& i_ps,i_uv,i_t,i_q,i_pw,i_rw,i_dw,i_gps,i_sst,i_tcp,i_lag, & i_gust,i_vis,i_pblh,i_wspd10m,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv, & - i_tcamt,i_lcbas,i_cldch,i_uwnd10m,i_vwnd10m,i_swcp,i_lwcp,i_dbz, & + i_tcamt,i_lcbas,i_cldch,i_uwnd10m,i_vwnd10m,i_swcp,i_lwcp,i_fed,i_dbz, & size(awork1,2),bwork1,awork1,ndata) ! Compute and print statistics for "lightning" data diff --git a/src/gsi/statsconv.f90 b/src/gsi/statsconv.f90 index a01675d8d0..c72adf34e7 100644 --- a/src/gsi/statsconv.f90 +++ b/src/gsi/statsconv.f90 @@ -2,7 +2,7 @@ subroutine statsconv(mype,& i_ps,i_uv,i_t,i_q,i_pw,i_rw,i_dw,i_gps,i_sst,i_tcp,i_lag, & i_gust,i_vis,i_pblh,i_wspd10m,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv, & i_tcamt,i_lcbas,i_cldch,i_uwnd10m,i_vwnd10m,& - i_swcp,i_lwcp,i_dbz,i_ref,bwork,awork,ndata) + i_swcp,i_lwcp,i_fed,i_dbz,i_ref,bwork,awork,ndata) !$$$ subprogram documentation block ! . . . . ! subprogram: statconv prints statistics for conventional data @@ -74,6 +74,7 @@ subroutine statsconv(mype,& ! i_vwnd10m- index in awork array holding vwnd10m info ! i_swcp - index in awork array holding swcp info ! i_lwcp - index in awork array holding lwcp info +! i_fed - index in awork array holding fed info ! i_dbz - index in awork array holding dbz info ! i_ref - size of second dimension of awork array ! bwork - array containing information for statistics @@ -96,12 +97,12 @@ subroutine statsconv(mype,& iout_gust,iout_vis,iout_pblh,iout_wspd10m,iout_td2m,& iout_mxtm,iout_mitm,iout_pmsl,iout_howv,iout_tcamt,iout_lcbas,iout_cldch,& iout_uwnd10m,iout_vwnd10m,& - iout_dbz,iout_swcp,iout_lwcp,& + iout_fed,iout_dbz,iout_swcp,iout_lwcp,& mype_dw,mype_rw,mype_sst,mype_gps,mype_uv,mype_ps,& mype_t,mype_pw,mype_q,mype_tcp,ndat,dtype,mype_lag,mype_gust,& mype_vis,mype_pblh,mype_wspd10m,mype_td2m,mype_mxtm,mype_mitm,& mype_pmsl,mype_howv,mype_tcamt,mype_lcbas,mype_cldch,mype_uwnd10m,mype_vwnd10m,& - mype_dbz,mype_swcp,mype_lwcp + mype_fed,mype_dbz,mype_swcp,mype_lwcp use qcmod, only: npres_print,ptop,pbot,ptopq,pbotq use jfunc, only: first,jiter use gridmod, only: nsig @@ -112,7 +113,7 @@ subroutine statsconv(mype,& integer(i_kind) ,intent(in ) :: mype,i_ps,i_uv,& i_t,i_q,i_pw,i_rw,i_dw,i_gps,i_sst,i_tcp,i_lag,i_gust,i_vis,i_pblh,& i_wspd10m,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv,i_tcamt,i_lcbas,& - i_cldch,i_uwnd10m,i_vwnd10m,i_swcp,i_lwcp,i_dbz,i_ref + i_cldch,i_uwnd10m,i_vwnd10m,i_swcp,i_lwcp,i_fed,i_dbz,i_ref real(r_kind),dimension(7*nsig+100,i_ref) ,intent(in ) :: awork real(r_kind),dimension(npres_print,nconvtype,5,3),intent(in ) :: bwork integer(i_kind),dimension(ndat,3) ,intent(in ) :: ndata @@ -136,6 +137,7 @@ subroutine statsconv(mype,& real(r_kind) dwqcplty,tqcplty,qctt,qctrw,rwqcplty,qctdw,qqcplty,qctgps real(r_kind) gpsqcplty,tpw3,pw3,qctq real(r_kind) tswcp3,tlwcp3,qctdbz,dbzqcplty + real(r_kind) fedmplty,tfed,qctfed,fedqcplty real(r_kind),dimension(1):: pbotall,ptopall logical,dimension(nconvtype):: pflag @@ -1325,6 +1327,68 @@ subroutine statsconv(mype,& end if end if +! Summary report for flash extent density + if(mype==mype_fed) then + nread=0 + nkeep=0 + do i=1,ndat + if(dtype(i)== 'fed')then + nread=nread+ndata(i,2) + nkeep=nkeep+ndata(i,3) + end if + end do + if(nread > 0)then + if(first)then + open(iout_fed) + else + open(iout_fed,position='append') + end if + + fedmplty=zero; fedqcplty=zero ; ntot=0 + tfed=zero ; qctfed=zero + if(nkeep > 0)then + mesage='current vfit of flash extent density, ranges in flashes per minute$' + do j=1,nconvtype + pflag(j)=trim(ioctype(j)) == 'fed' + end do + call dtast(bwork,npres_print,pbot,ptop,mesage,jiter,iout_fed,pflag) + + numgross=nint(awork(4,i_fed)) + numfailqc=nint(awork(21,i_fed)) + do k=1,nsig + num(k)=nint(awork(k+5*nsig+100,i_fed)) + rat=zero + rat3=zero + if(num(k) > 0) then + rat=awork(6*nsig+k+100,i_fed)/float(num(k)) + rat3=awork(3*nsig+k+100,i_fed)/float(num(k)) + end if + ntot=ntot+num(k) + fedmplty=fedmplty+awork(6*nsig+k+100,i_fed) + fedqcplty=fedqcplty+awork(3*nsig+k+100,i_fed) + write(iout_fed,240) 'r',num(k),k,awork(6*nsig+k+100,i_fed), & + awork(3*nsig+k+100,i_fed),rat,rat3 + end do + if(ntot > 0) then + tfed=fedmplty/float(ntot) + qctfed=fedqcplty/float(ntot) + end if + write(iout_fed,925) 'fed',numgross,numfailqc + numlow = nint(awork(2,i_fed)) + numhgh = nint(awork(3,i_fed)) + nhitopo = nint(awork(5,i_fed)) + ntoodif = nint(awork(6,i_fed)) + write(iout_fed,900) 'fed',numhgh,numlow + write(iout_fed,905) 'fed',nhitopo,ntoodif + end if + write(iout_fed,950) 'fed',jiter,nread,nkeep,ntot + write(iout_fed,951) 'fed',fedmplty,fedqcplty,tfed,qctfed + + close(iout_fed) + end if + end if + + if(mype==mype_tcp) then nread=0 nkeep=0 From f9433d8bc69fa35e09a26b6325b05205fff73754 Mon Sep 17 00:00:00 2001 From: David Dowell Date: Tue, 22 Aug 2023 14:26:08 +0000 Subject: [PATCH 4/7] Updated code, with changes suggested by reviewers. --- src/gsi/gsi_fedOper.F90 | 17 ++----- src/gsi/m_fedNode.F90 | 21 +------- src/gsi/read_fed.f90 | 10 ++-- src/gsi/setupfed.f90 | 110 ++++++++++++++++++++-------------------- src/gsi/statsconv.f90 | 2 +- 5 files changed, 66 insertions(+), 94 deletions(-) diff --git a/src/gsi/gsi_fedOper.F90 b/src/gsi/gsi_fedOper.F90 index a306868f77..b2b2400ff0 100644 --- a/src/gsi/gsi_fedOper.F90 +++ b/src/gsi/gsi_fedOper.F90 @@ -2,15 +2,13 @@ module gsi_fedOper !$$$ subprogram documentation block ! . . . . ! subprogram: module gsi_fedOper -! prgmmr: j guo -! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.3 -! date: 2018-08-10 ! ! abstract: an obOper extension for fedNode type ! ! program history log: -! 2023-04-10 D. Dowell - moved diag_fed and its description here from -! obsmod. +! 2023-07-10 D. Dowell - created new module for FED (flash extent +! density); gsi_dbzOper.F90 code used as a +! starting point for developing this new module ! ! input argument list: see Fortran 90 style document below ! @@ -130,7 +128,6 @@ subroutine setup_(self, lunin, mype, is, nobs, init_pass, last_pass) end subroutine setup_ subroutine intjo1_(self, ibin, rval,sval, qpred,sbias) -! use intfedmod, only: intjo => intfed use gsi_bundlemod , only: gsi_bundle use bias_predictors, only: predictors use m_obsNode , only: obsNode @@ -148,14 +145,9 @@ subroutine intjo1_(self, ibin, rval,sval, qpred,sbias) character(len=*),parameter:: myname_=myname//"::intjo1_" class(obsNode),pointer:: headNode -! headNode => obsLList_headNode(self%obsLL(ibin)) -! call intjo(headNode, rval,sval) -! headNode => null() - end subroutine intjo1_ subroutine stpjo1_(self, ibin, dval,xval,pbcjo,sges,nstep,dbias,xbias) -! use stpfedmod, only: stpjo => stpfed use gsi_bundlemod, only: gsi_bundle use bias_predictors, only: predictors use m_obsNode , only: obsNode @@ -177,9 +169,6 @@ subroutine stpjo1_(self, ibin, dval,xval,pbcjo,sges,nstep,dbias,xbias) character(len=*),parameter:: myname_=myname//"::stpjo1_" class(obsNode),pointer:: headNode -! headNode => obsLList_headNode(self%obsLL(ibin)) -! call stpjo(headNode,dval,xval,pbcjo(:),sges,nstep) -! headNode => null() end subroutine stpjo1_ end module gsi_fedOper diff --git a/src/gsi/m_fedNode.F90 b/src/gsi/m_fedNode.F90 index f42503df19..84a319cd12 100644 --- a/src/gsi/m_fedNode.F90 +++ b/src/gsi/m_fedNode.F90 @@ -36,27 +36,21 @@ module m_fedNode public:: fedNode type,extends(obsNode):: fedNode - !type(td_ob_type),pointer :: llpoint => NULL() type(obs_diag), pointer :: diags => NULL() real(r_kind) :: res ! flash extent density residual real(r_kind) :: err2 ! flash extent density 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) :: jb ! variational quality control parameter real(r_kind) :: wij(8) ! horizontal interpolation weights real(r_kind) :: fedpertb ! random number adding to the obs -! logical :: luse ! flag indicating if ob is used in pen. integer(i_kind) :: k1 ! level of errtable 1-33 integer(i_kind) :: kx ! ob type integer(i_kind) :: ij(8) ! horizontal locations -! integer(i_kind) :: idv,iob ! device id and obs index for sorting real (r_kind) :: dlev ! reference to the vertical grid - !real (r_kind) :: elat, elon ! earth lat-lon for redistribution - !real (r_kind) :: dlat, dlon ! earth lat-lon for redistribution contains procedure,nopass:: mytype procedure:: setHop => obsNode_setHop_ @@ -65,10 +59,6 @@ module m_fedNode procedure:: isvalid => obsNode_isvalid_ procedure:: gettlddp => gettlddp_ - ! procedure, nopass:: headerRead => obsHeader_read_ - ! procedure, nopass:: headerWrite => obsHeader_write_ - ! procedure:: init => obsNode_init_ - ! procedure:: clean => obsNode_clean_ end type fedNode public:: fedNode_typecast @@ -81,8 +71,6 @@ module m_fedNode character(len=*),parameter:: MYNAME="m_fedNode" -!#define CHECKSUM_VERBOSE -!#define DEBUG_TRACE #include "myassert.H" #include "mytrace.H" contains @@ -92,14 +80,11 @@ function typecast_(aNode) result(ptr_) implicit none type(fedNode),pointer:: ptr_ class(obsNode),pointer,intent(in):: aNode -! character(len=*),parameter:: myname_=MYNAME//"::typecast_" ptr_ => null() if(.not.associated(aNode)) return select type(aNode) - type is(fedNode) - ptr_ => aNode -! class default -! call die(myname_,'unexpected type, aNode%mytype() =',aNode%mytype()) + type is(fedNode) + ptr_ => aNode end select return end function typecast_ @@ -171,7 +156,6 @@ subroutine obsNode_xread_(aNode,iunit,istat,diagLookup,skip) aNode%b , & aNode%pg , & aNode%jb , & - ! aNode%fedpertb , & aNode%k1 , & aNode%kx , & aNode%dlev , & @@ -210,7 +194,6 @@ subroutine obsNode_xwrite_(aNode,junit,jstat) aNode%b , & aNode%pg , & aNode%jb , & - ! aNode%fedpertb , & aNode%k1 , & aNode%kx , & aNode%dlev , & diff --git a/src/gsi/read_fed.f90 b/src/gsi/read_fed.f90 index 79c6e2a726..849f1b603d 100644 --- a/src/gsi/read_fed.f90 +++ b/src/gsi/read_fed.f90 @@ -370,8 +370,8 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) allocate(cdata_out(nreal,ndata)) ! ! - DO i=1,numfed - DO k=1,maxlvl + do i=1,numfed + do k=1,maxlvl ! DCD 1 July 2021 if (fed3d_column(k+2,i) .gt. fed_highbnd) fed3d_column(k+2,i) = fed_highbnd @@ -381,7 +381,7 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) write(6,*) ' ------- check max and min value of OBS: bufr fed -------' write(6,*) ' level maxval(fed) minval(fed)' - DO k=1,maxlvl + do k=1,maxlvl write(6,*) k,maxval(fed3d_column(k+2,:)),minval(fed3d_column(k+2,:)) end do @@ -392,8 +392,8 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) kint_maxloc=-1 fed_max=-999.99 ndata2=0 - DO i=1,numfed - DO k=1,maxlvl + do i=1,numfed + do k=1,maxlvl if( fed3d_column(k+2,i) >= fed_lowbnd2 .or. fed3d_column(k+2,i) == fed_lowbnd) then !Rong Kong dlon_earth = fed3d_column(1,i) ! longitude (degrees) of observation ! ilone=18 ! index of longitude (degrees) diff --git a/src/gsi/setupfed.f90 b/src/gsi/setupfed.f90 index 6432791443..350a2a68ec 100644 --- a/src/gsi/setupfed.f90 +++ b/src/gsi/setupfed.f90 @@ -193,26 +193,26 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa !------------------------------------------------! - integer(i_kind) :: ibgn, iend, jbgn, jend, ips,ipe,jps,jpe,itmp,jtmp,ktmp + integer(i_kind) :: ibgn, iend, jbgn, jend, ips,ipe,jps,jpe,itmp,jtmp,ktmp character(256) :: binfilename integer(i_kind), parameter :: ntimesfed=1 - character(256) ::fedfilename - integer(i_kind),parameter :: nxfed=99, nyfed=99, nzfed=1, nfldfed=3 + character(256) :: fedfilename + integer(i_kind),parameter :: nxfed=99, nyfed=99, nzfed=1, nfldfed=3 real(4) :: a(ntimesfed,nfldfed,nzfed,nyfed,nxfed) real(4) :: gga(ntimesfed,nfldfed,nzfed,nyfed,nxfed) integer(i_kind) irec1, irec2, irec3, irec4, itot - integer(i_kind) :: la, iobs, lt, nnnnn + integer(i_kind) :: la, iobs, lt, nnnnn real(r_kind),dimension(nobs) :: FEDMdiag,FEDMdiagTL real(r_kind),dimension(nobs) :: FEDMdiag2D - integer(i_kind) :: npt - integer(i_kind) :: nobsfed - real(r_kind) :: dlat_earth,dlon_earth - logical :: outside + integer(i_kind) :: npt + integer(i_kind) :: nobsfed + real(r_kind) :: dlat_earth,dlon_earth + logical :: outside ! YPW added the next lines - logical :: l_set_oerr_ratio_fed=.False. - logical :: l_gpht2gmht = .True. + logical :: l_set_oerr_ratio_fed=.False. + logical :: l_gpht2gmht = .True. integer(i_kind) :: ncid,status,x_dimid,y_dimid,z_dimid,varid,x_varid,y_varid integer(i_kind),dimension(3):: dimids integer(i_kind),dimension(2):: dimids_2d @@ -341,8 +341,8 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa rp(j,i,jj)=rp(j,i,jj) + ges_qg(jtmp,itmp,k,jj)* & dx_m*dy_m*(ges_prsi(jtmp,itmp,k,jj)-ges_prsi(jtmp,itmp,k+1,jj))*& tpwcon * r10 - enddo !igx - enddo !jgy + end do !igx + end do !jgy end do !j end do !i end do !k @@ -419,7 +419,7 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa ibin = NINT( dtime/hr_obsbin ) + 1 else ibin = 1 - endif + end if IF (ibin<1.OR.ibin>nobs_bins) write(6,*)mype,'Error nobs_bins,ibin= ',nobs_bins,ibin @@ -438,7 +438,7 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa miter = miter ) if(.not.associated(my_diag)) call die(myname, & 'obsdiagLList_nextNode(), create =',.not.lobsdiag_allocated) - endif + end if ! Interpolate terrain height(model elevation) to obs location. call tintrp2a11(ges_z,zsges,dlat,dlon,dtime,hrdifsig,& @@ -641,7 +641,7 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa term = exp_arg wgt = wgtlim rwgt = wgt/wgtlim - endif + end if valqc = -two*rat_err2*term ! print*,'Compute penalty terms' @@ -718,19 +718,19 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa if(presq > ptabl(k+1) .and. presq <= ptabl(k)) then my_head%k1=k exit k_loop - endif - enddo k_loop - endif - endif + end if + end do k_loop + end if + end if !------------------------------------------------- if(luse_obsdiag)then call obsdiagNode_assert(my_diag,my_head%idv,my_head%iob,1,myname,'my_diag:my_head') my_head%diags => my_diag - endif + end if my_head => null() - endif + end if ! Save select output for diagnostic file if(.not.luse(i))write(6,*)' luse, mype',luse(i),mype @@ -743,7 +743,7 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa err_final = one/(ratio_errors*error) else err_final = huge_single - endif + end if errinv_input = huge_single errinv_adjst = huge_single errinv_final = huge_single @@ -775,13 +775,13 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa open(66,file=trim(diag_file),form='unformatted',status='old',position='append') else open(66,file=trim(diag_file),form='unformatted',status='unknown',position='rewind') - endif - endif + end if + end if if(init_pass .and. mype == 0) then write(66) ianldate write(6,*)'SETUPFED: write time record to file ',& trim(diag_file), ' ',ianldate - endif + end if ! call dtime_show(myname,'diagsave:fed',i_fed_ob_type) write(66)'fed',nchar,nreal,ii,mype,ioff0 @@ -849,17 +849,17 @@ subroutine init_vars_ if(allocated(ges_ps))then write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' call stop2(999) - endif + end if allocate(ges_ps(size(rank2,1),size(rank2,2),nfldsig)) ges_ps(:,:,1)=rank2 do ifld=2,nfldsig call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank2,istatus) ges_ps(:,:,ifld)=rank2 - enddo + end do else write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus call stop2(999) - endif + end if ! get z ... varname='z' call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus) @@ -867,17 +867,17 @@ subroutine init_vars_ if(allocated(ges_z))then write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' call stop2(999) - endif + end if allocate(ges_z(size(rank2,1),size(rank2,2),nfldsig)) ges_z(:,:,1)=rank2 do ifld=2,nfldsig call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank2,istatus) ges_z(:,:,ifld)=rank2 - enddo + end do else write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus call stop2(999) - endif + end if ! get q ... varname='q' call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) @@ -885,17 +885,17 @@ subroutine init_vars_ if(allocated(ges_q))then write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' call stop2(999) - endif + end if allocate(ges_q(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) ges_q(:,:,:,1)=rank3 do ifld=2,nfldsig call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) ges_q(:,:,:,ifld)=rank3 - enddo + end do else write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus call stop2(999) - endif + end if ! get tv ... ! varname='tv' ! call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) @@ -903,18 +903,18 @@ subroutine init_vars_ ! if(allocated(ges_tv))then ! write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' ! call stop2(999) -! endif +! end if ! allocate(ges_tv(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) ! ges_tv(:,:,:,1)=rank3 ! do ifld=2,nfldsig ! call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) ! ges_tv(:,:,:,ifld)=rank3 ! ges_tv(:,:,:,ifld)=rank3 -! enddo +! end do ! else ! write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus ! call stop2(999) -! endif +! end if ! get qr ... ! get qg ... varname='qg' @@ -923,26 +923,26 @@ subroutine init_vars_ if(allocated(ges_qg))then write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' call stop2(999) - endif + end if allocate(ges_qg(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) if(.not. allocated(ges_qg_mask))then allocate(ges_qg_mask(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) - endif + end if ges_qg(:,:,:,1)=rank3 do ifld=2,nfldsig call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) ges_qg(:,:,:,ifld)=rank3 - enddo + end do else write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus call stop2(999) - endif + end if else write(6,*) trim(myname), ': inconsistent vector sizes (nfldsig,size(metguess_bundle) ',& nfldsig,size(gsi_metguess_bundle) call stop2(999) - endif + end if end subroutine init_vars_ subroutine init_netcdf_diag_ @@ -967,7 +967,7 @@ subroutine init_netcdf_diag_ else if(verbose) print *,'file ' // trim(diag_conv_file) // ' exists but contains no obs. Not appending. nobs,mype=',ncd_nobs,mype append_diag = .false. ! if there are no obs in existing file, then do not try to append - endif + end if end if call nc_diag_init(diag_conv_file, append=append_diag) @@ -975,7 +975,7 @@ subroutine init_netcdf_diag_ if (.not. append_diag) then ! don't write headers on append - the module will break? call nc_diag_header("date_time",ianldate ) call nc_diag_header("Number_of_state_vars", nsdim ) - endif + end if end subroutine init_netcdf_diag_ subroutine contents_binary_diag_(odiag) type(obs_diag),pointer,intent(in):: odiag @@ -999,7 +999,7 @@ subroutine contents_binary_diag_(odiag) rdiagbuf(12,ii) = one ! analysis usage flag (1=use, -1=not used) else rdiagbuf(12,ii) = -one - endif + end if rdiagbuf(13,ii) = rwgt ! nonlinear qc relative weight rdiagbuf(14,ii) = errinv_input ! prepbufr inverse obs error (dBZ)**-1 @@ -1028,21 +1028,21 @@ subroutine contents_binary_diag_(odiag) rdiagbuf(ioff,ii) = one else rdiagbuf(ioff,ii) = -one - endif - enddo + end if + end do do jj=1,miter+1 ioff=ioff+1 rdiagbuf(ioff,ii) = odiag%nldepart(jj) - enddo + end do do jj=1,miter ioff=ioff+1 rdiagbuf(ioff,ii) = odiag%tldepart(jj) - enddo + end do do jj=1,miter ioff=ioff+1 rdiagbuf(ioff,ii) = odiag%obssen(jj) - enddo - endif + end do + end if end subroutine contents_binary_diag_ subroutine contents_netcdf_diag_(odiag) @@ -1068,7 +1068,7 @@ subroutine contents_netcdf_diag_(odiag) call nc_diag_metadata("Analysis_Use_Flag", sngl(one) ) else call nc_diag_metadata("Analysis_Use_Flag", sngl(-one) ) - endif + end if call nc_diag_metadata("Errinv_Input", sngl(errinv_input) ) call nc_diag_metadata("Errinv_Adjust", sngl(errinv_adjst) ) @@ -1083,14 +1083,14 @@ subroutine contents_netcdf_diag_(odiag) obsdiag_iuse(jj) = one else obsdiag_iuse(jj) = -one - endif - enddo + end if + end do call nc_diag_data2d("ObsDiagSave_iuse", obsdiag_iuse ) call nc_diag_data2d("ObsDiagSave_nldepart", odiag%nldepart ) call nc_diag_data2d("ObsDiagSave_tldepart", odiag%tldepart ) call nc_diag_data2d("ObsDiagSave_obssen" , odiag%obssen ) - endif + end if end subroutine contents_netcdf_diag_ diff --git a/src/gsi/statsconv.f90 b/src/gsi/statsconv.f90 index c72adf34e7..0da8606f24 100644 --- a/src/gsi/statsconv.f90 +++ b/src/gsi/statsconv.f90 @@ -1340,7 +1340,7 @@ subroutine statsconv(mype,& if(nread > 0)then if(first)then open(iout_fed) - else + else open(iout_fed,position='append') end if From d3e64d34e17733acaf22f4ba5d2e79558358e5f7 Mon Sep 17 00:00:00 2001 From: David Dowell Date: Wed, 23 Aug 2023 18:00:16 +0000 Subject: [PATCH 5/7] Bug fix requested by Guoqing Ge and Chunhua Zhou. --- src/gsi/constants.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gsi/constants.f90 b/src/gsi/constants.f90 index b4cf775068..2d0d53a8ad 100644 --- a/src/gsi/constants.f90 +++ b/src/gsi/constants.f90 @@ -90,7 +90,7 @@ module constants ! Declare derived constants integer(i_kind):: huge_i_kind - integer(i_kind), parameter :: max_varname_length=20 + integer(i_kind), parameter :: max_varname_length=60 real(r_single):: tiny_single, huge_single real(r_kind):: xai, xa, xbi, xb, dldt, rozcon,ozcon,fv, tpwcon,eps, rd_over_g real(r_kind):: el2orc, g_over_rd, rd_over_cp, cpr, omeps, epsm1, factor2 From 278ee2f6d3026a1bcba6226df94b933c53a4e8da Mon Sep 17 00:00:00 2001 From: "Ming.Hu" Date: Fri, 1 Sep 2023 22:53:43 +0000 Subject: [PATCH 6/7] Clean read_fed.f90 and setupfed.f90: remove unused variables found by DEBUG mode. --- src/gsi/read_fed.f90 | 13 +++---- src/gsi/setupfed.f90 | 88 +++++++++----------------------------------- 2 files changed, 23 insertions(+), 78 deletions(-) diff --git a/src/gsi/read_fed.f90 b/src/gsi/read_fed.f90 index 849f1b603d..8aaafeba92 100644 --- a/src/gsi/read_fed.f90 +++ b/src/gsi/read_fed.f90 @@ -29,14 +29,13 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) !_____________________________________________________________________ ! use kinds, only: r_kind,r_double,i_kind - use constants, only: zero,one,rad2deg,deg2rad - use convinfo, only: nconvtype,ctwind,cgross,cermax,cermin,cvar_b,cvar_pg, & - ncmiter,ncgroup,ncnumgrp,icuse,ictype,icsubtype,ioctype + use constants, only: zero,one,deg2rad + use convinfo, only: nconvtype,ctwind,icuse,ioctype use gsi_4dvar, only: l4dvar,l4densvar,winlen use gridmod, only: tll2xy use mod_wrfmass_to_a, only: wrfmass_obs_to_a8 use mpimod, only: npe - use obsmod, only: perturb_obs,iadatemn,time_window + use obsmod, only: perturb_obs,iadatemn use netcdf implicit none @@ -67,7 +66,6 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) integer(i_kind) ifn,i - real(r_kind) :: maxfed integer(i_kind) :: ilon,ilat logical :: fedobs, fedob @@ -97,7 +95,7 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) integer(i_kind) :: ireadmg,ireadsb integer(i_kind) :: maxlvl - integer(i_kind) :: numlvl,numfed,numobsa,nmsgmax,maxobs + integer(i_kind) :: numlvl,numfed,nmsgmax,maxobs integer(i_kind) :: k,iret integer(i_kind) :: nmsg,ntb @@ -119,9 +117,8 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) integer :: unit_table ! for read netcdf - integer(i_kind) :: idate5(5), sec70,mins_an,mins_ob + integer(i_kind) :: sec70,mins_an integer(i_kind) :: varID, ncdfID, status - character(4) :: idate5s(5) real(r_kind) :: timeb,twindm,rmins_an,rmins_ob diff --git a/src/gsi/setupfed.f90 b/src/gsi/setupfed.f90 index 350a2a68ec..cf6334e567 100644 --- a/src/gsi/setupfed.f90 +++ b/src/gsi/setupfed.f90 @@ -41,16 +41,16 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa use gsi_4dvar, only: nobs_bins,hr_obsbin use oneobmod, only: oneobtest,maginnov,magoberr use guess_grids, only: hrdifsig,nfldsig,ges_prsi - use guess_grids, only: ges_lnprsl, ges_prsl, ges_tsen, geop_hgtl + use guess_grids, only: ges_lnprsl, geop_hgtl use gridmod, only: lat2, lon2 - use gridmod, only: nsig, get_ij,get_ijk,regional,tll2xy + use gridmod, only: nsig, get_ij,get_ijk,tll2xy use constants, only: flattening,semi_major_axis,grav_ratio,zero,grav,wgtlim use constants, only: half,one,two,grav_equator,eccentricity,somigliana - use constants, only: rad2deg,deg2rad,r60,tiny_r_kind,cg_term,huge_single + use constants, only: deg2rad,r60,tiny_r_kind,cg_term,huge_single use constants, only: r10,r100,r1000 - use constants, only: rd,grav,tpwcon - use qcmod, only: npres_print,ptop,pbot,ptopq,pbotq,dfact,dfact1 - use jfunc, only: jiter,last,jiterstart,miter + use constants, only: grav,tpwcon + use qcmod, only: npres_print,ptopq,pbotq + use jfunc, only: jiter,last,miter use convinfo, only: nconvtype,cermin,cermax,cgross,cvar_b,cvar_pg,ictype use convinfo, only: icsubtype use converr, only: ptabl @@ -100,8 +100,6 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa real(r_kind),parameter:: D608=0.608_r_kind character(len=*),parameter:: myname='setupfed' - real(r_kind) :: Cs_tmp, Cg_tmp ! temporary coefficients for check-up - ! Declare external calls for code analysis external:: tintrp2a1 external:: tintrp2a11 @@ -112,10 +110,8 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa ! Declare local variables real(r_kind) rlow,rhgh,rsig - real(r_kind) dz,denom + real(r_kind) dz real(r_kind) jqg_num,jqg - real(r_kind) wgt_dry, wgt_wet - real(r_kind) jqg_num_dry, jqg_num_wet real(r_kind) dlnp,pobl,zob real(r_kind) sin2,termg,termr,termrg real(r_kind) psges,zsges @@ -129,13 +125,9 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa real(r_double) rstation_id real(r_kind) dlat,dlon,dtime,dpres,ddiff,error,slat,dlat8km,dlon8km real(r_kind) ratio_errors - real(r_kind) qgges,rhoges - real(r_kind) Ze,rdBZ,presw,fednoise,fednoise_runits - real(r_kind) Ze_orig, Zer, Zes, Zeg - real(r_kind) Zeg_dry, Zeg_wet + real(r_kind) presw real(r_kind) errinv_input,errinv_adjst,errinv_final real(r_kind) err_input,err_adjst,err_final - real(r_kind) qgexp real(r_kind),dimension(nele,nobs):: data real(r_kind),dimension(lat2,lon2,nfldsig)::rp real(r_single),allocatable,dimension(:,:)::rdiagbuf @@ -145,29 +137,22 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa real(r_kind),allocatable,dimension(:,:,:,:) :: ges_qg,ges_qg_mask real(r_kind) :: presq - real(r_kind) :: P1D,T1D,Q1D,RHO - real(r_kind) :: qges,tsenges ! used to calculate tv - virtual temperature - real(r_kind) :: lnprslges ! use log(p) for vertical interpolation - real(r_kind) :: qg_min + real(r_kind) :: T1D,RHO real(r_kind) :: glmcoeff = 2.088_r_kind*10.0**(-8.0) ! Allen et al. (2016,MWR) real(r_kind) :: CM = 0.5_r_kind ! tuning factor in eq. 14 of Kong et al. 2020 - integer(i_kind) i,nchar,nreal,k,j,k1,ii,nii,jj,im,jm,km - integer(i_kind) mm1,k2,isli + integer(i_kind) i,nchar,nreal,k,j,k1,ii,jj + integer(i_kind) mm1,k2 integer(i_kind) jsig,ikxx,nn,ibin,ioff,ioff0 - integer(i_kind) ier,ilat,ilon,ihgt,ifedob,ikx,itime,iuse - integer(i_kind) ielev,id,itilt,iazm,ilone,ilate,irange - integer(i_kind) ier2,ifednoise,it,istatus - integer(i_kind) ier_b - integer(i_kind) ijk + integer(i_kind) ier,ilat,ilon,ifedob,ikx,itime,iuse + integer(i_kind) id,ilone,ilate + integer(i_kind) ier2 - integer(i_kind) i4,j4,k4,n4 integer(i_kind) nlat_ll,nlon_ll,nsig_ll,nfld_ll integer(i_kind) ipres,iqmax,iqc,icat,itemp integer(i_kind) istnelv,iobshgt,izz,iprvd,isprvd,iptrb integer(i_kind) idomsfc,iskint,isfcr,iff10 - integer(i_kind) nguess character(8) station_id character(8),allocatable,dimension(:):: cdiagbuf @@ -179,11 +164,8 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa logical proceed equivalence(rstation_id,station_id) - real(r_kind) wrange - integer(i_kind) numequal,numnotequal,kminmin,kmaxmax,istat + integer(i_kind) numequal,numnotequal - logical:: in_curbin, in_anybin - type(fedNode),pointer:: my_head type(obs_diag),pointer:: my_diag type(obs_diags),pointer:: my_diagLL @@ -193,33 +175,18 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa !------------------------------------------------! - integer(i_kind) :: ibgn, iend, jbgn, jend, ips,ipe,jps,jpe,itmp,jtmp,ktmp - character(256) :: binfilename + integer(i_kind) :: itmp,jtmp integer(i_kind), parameter :: ntimesfed=1 - character(256) :: fedfilename integer(i_kind),parameter :: nxfed=99, nyfed=99, nzfed=1, nfldfed=3 - real(4) :: a(ntimesfed,nfldfed,nzfed,nyfed,nxfed) - real(4) :: gga(ntimesfed,nfldfed,nzfed,nyfed,nxfed) - integer(i_kind) irec1, irec2, irec3, irec4, itot - integer(i_kind) :: la, iobs, lt, nnnnn real(r_kind),dimension(nobs) :: FEDMdiag,FEDMdiagTL - real(r_kind),dimension(nobs) :: FEDMdiag2D integer(i_kind) :: npt - integer(i_kind) :: nobsfed real(r_kind) :: dlat_earth,dlon_earth - logical :: outside ! YPW added the next lines logical :: l_set_oerr_ratio_fed=.False. logical :: l_gpht2gmht = .True. - integer(i_kind) :: ncid,status,x_dimid,y_dimid,z_dimid,varid,x_varid,y_varid - integer(i_kind),dimension(3):: dimids - integer(i_kind),dimension(2):: dimids_2d - character(256) :: outfile real(r_kind),dimension(nobs) :: dlatobs,dlonobs - real(r_kind),dimension(4):: wgrd - integer(i_kind),dimension(4):: jgrd integer(i_kind):: ngx,ngy,igx,jgy real(r_kind):: dx_m, dy_m @@ -318,7 +285,7 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa ngy = 2 dx_m = 3000. dy_m = 3000. - print*,'Operator start here!,ngx=',ngx,'ngy=',ngy + print*,'FED Operator start here!,ngx=',ngx,'ngy=',ngy rp=zero print*, 'mype = ', mype @@ -377,10 +344,6 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa end if write(6,*) 'fed_highbnd=',fed_highbnd write(6,*) 'maxval(ges_qg)=',maxval(ges_qg),'pe=',mype -! write(6,*) 'maxval(geop_hgtl)=',maxval(geop_hgtl(:,:,:,it)) - write(6,*) 'maxval(ges_tsen)=',maxval(ges_tsen(:,:,:,it)) - write(6,*) 'maxval(FED)=',maxval(rp) - write(6,*) 'ges_prsi',ges_prsi(100,100,1,1),ges_prsi(100,100,nsig,1) !============================================================================================ @@ -412,8 +375,6 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa dlon_earth = data(ilone,i) !the lontitude and latitude on the obs pts. dlat_earth = data(ilate,i) ! geometric hgh (hges --> zges below) -! print*,'i,mype,dlat,dlon,dlon8km,dlat8km',i,mype,dlat,dlon,dlon8km,dlat8km,& -! dlon_earth,dlat_earth,dpres,data(ifedob,i) if (nobs_bins>1) then ibin = NINT( dtime/hr_obsbin ) + 1 @@ -693,14 +654,11 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa my_head%elat= data(ilate,i) my_head%elon= data(ilone,i) - if (istat/=0) write(6,*)'MAKECOBS: allocate error for fedtail_dzg,istat=',istat - my_head%dlev= dpres - call get_ijk(mm1,dlat,dlon,dpres,my_head%ij,my_head%wij) + call get_ijk(mm1,dlat,dlon,dpres,my_head%ij,my_head%wij) my_head%res = ddiff ! Observation - ges my_head%err2 = error**2 my_head%raterr2 = ratio_errors**2 -! my_head%jqg = jqg ! for TL and ADJ my_head%time = dtime my_head%b = cvar_b(ikx) my_head%pg = cvar_pg(ikx) @@ -806,28 +764,18 @@ subroutine check_vars_ (proceed) ! Check to see if required guess fields are available call gsi_metguess_get ('var::ps', ivar, istatus ) proceed=ivar>0 - print*,'For ps, proceed=',proceed call gsi_metguess_get ('var::z' , ivar, istatus ) proceed=proceed.and.ivar>0 - print*,'For z, proceed=',proceed call gsi_metguess_get ('var::q' , ivar, istatus ) proceed=proceed.and.ivar>0 - print*,'For q, proceed=',proceed ! call gsi_metguess_get ('var::tv' , ivar, istatus ) ! proceed=proceed.and.ivar>0 call gsi_metguess_get ('var::qs', ivar, istatus ) proceed=proceed.and.ivar>0 - print*,'For qs, proceed=',proceed call gsi_metguess_get ('var::qg', ivar, istatus ) proceed=proceed.and.ivar>0 - print*,'For qg, proceed=',proceed call gsi_metguess_get ('var::qr', ivar, istatus ) proceed=proceed.and.ivar>0 - print*,'For qr, proceed=',proceed -! if ( mphyopt == 108 ) then ! comment out by YPW -! proceed=proceed.and.ivar>0 ! comment out by YPW -! print*,'For qnr, proceed=',proceed ! comment out by YPW -! end if ! comment out by YPW end subroutine check_vars_ From 57ae9acd34478e04a06dcfa972cdf96fb026e261 Mon Sep 17 00:00:00 2001 From: "Ming.Hu" Date: Tue, 5 Sep 2023 17:28:52 +0000 Subject: [PATCH 7/7] Add boundary check for lat and lon. --- src/gsi/read_fed.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/gsi/read_fed.f90 b/src/gsi/read_fed.f90 index 8aaafeba92..c478b3d93f 100644 --- a/src/gsi/read_fed.f90 +++ b/src/gsi/read_fed.f90 @@ -204,7 +204,7 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) ! Extract type, date, and location information from BUFR file call ufbint(lunin,hdr,5,1,iret,hdrstr) - if(hdr(3) .gt. 90 ) write(6,*) "Inside read_fed.f90, hdr(2)=",hdr(2),"hdr(3)=",hdr(3) + if(hdr(3) .gt. r90 ) write(6,*) "Inside read_fed.f90, hdr(2)=",hdr(2),"hdr(3)=",hdr(3) if ( l_latlon_fedobs ) then if(abs(hdr(3))>r90 .or. abs(hdr(2))>r360) cycle loop_report if(hdr(2)== r360)hdr(2)=hdr(2)-r360 @@ -399,6 +399,7 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) !-Check format of longitude and correct if necessary if(dlon_earth>=r360) dlon_earth=dlon_earth-r360 if(dlon_earth=r360 .or. dlat_earth >90.0_r_kind) cycle !-Convert back to radians rlon00 = dlon_earth*deg2rad