!------------------------------------------------------------------------- ! NOAA/NCEP, National Centers for Environmental Prediction GSI ! !------------------------------------------------------------------------- !BOP ! ! !ROUTINE: setuppcp --- Compute rhs of oi equation for precipitation ! ! !INTERFACE: ! subroutine setuppcp(lunin,mype,aivals,nele,nobs,& obstype,isis,is,pcp_diagsave,init_pass) ! !USES: use mpeu_util, only: die,perr use kinds, only: r_kind,r_single,i_kind use pcpinfo, only: deltim use pcpinfo, only: dtphys use pcpinfo, only: npcptype use pcpinfo, only: iusep use pcpinfo, only: xkt2d use pcpinfo, only: npredp use pcpinfo, only: predxp use pcpinfo, only: varchp use pcpinfo, only: gross_pcp use pcpinfo, only: b_pcp use pcpinfo, only: pg_pcp use pcpinfo, only: ibias use pcpinfo, only: nupcp use gridmod, only: nlat ! no. lat use gridmod, only: nsig,nsig2,nsig3,nsig4,nsig5 ! no. levels use gridmod, only: rbs2 ! 1./sin(grid latitudes))**2 use gridmod, only: istart ! start lat of the whole array on each pe use gridmod, only: get_ij use guess_grids, only: geop_hgtl,hrdifsig,nfldsig use guess_grids, only: ges_prsl,ges_prsi,ges_tsen use guess_grids, only: isli2 use derivsmod, only: drv_initialized use derivsmod, only: gsi_xderivative_bundle use derivsmod, only: gsi_yderivative_bundle use tendsmod, only: tnd_initialized use tendsmod, only: gsi_tendency_bundle use obsmod, only: ndat,dplat,time_offset 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_obsdiags, only: pcphead use obsmod, only: time_offset use obsmod, only: i_pcp_ob_type,obsdiags,lobsdiagsave,ianldate use obsmod, only: mype_diaghdr,nobskeep,lobsdiag_allocated,dirname use m_obsNode, only: obsNode use m_pcpNode, only: pcpNode use m_obsLList, only: obsLList_appendNode use obsmod, only: obs_diag,luse_obsdiag use gsi_4dvar, only: nobs_bins,hr_obsbin,l4dvar,l4densvar use gsi_metguess_mod, only: gsi_metguess_bundle use gsi_bundlemod, only: gsi_bundlegetpointer use constants, only: rd,cp,pi,zero,quarter,r60, & half,one,two,three,tiny_r_kind,one_tenth,cg_term,r1000,wgtlim,fv,r3600,r10 use jfunc, only: jiter,miter use m_dtime, only: dtime_setup, dtime_check, dtime_show use gsi_bundlemod, only : gsi_bundlegetpointer use gsi_metguess_mod, only : gsi_metguess_get,gsi_metguess_bundle implicit none ! Turn off implicit typing ! !INPUT PARAMETERS: integer(i_kind) , intent(in ) :: lunin ! unit from which to read ! precpitation observations integer(i_kind) , intent(in ) :: mype ! mpi task id integer(i_kind) , intent(in ) :: nele ! number of pieces of information ! per precipitation observation integer(i_kind) , intent(in ) :: nobs ! number of precipitation obs to process character(len=20) , intent(in ) :: isis ! sensor/instrument/satellite id integer(i_kind) , intent(in ) :: is ! counter for number of obs types to process character(10) , intent(in ) :: obstype ! type of precipitation observation logical , intent(in ) :: pcp_diagsave ! switch diagnostic output on/off ! (.false.=no output) logical , intent(in ) :: init_pass ! state of "setup" processing ! !INPUT/OUTPUT PARAMETERS: real(r_kind),dimension(40,ndat), intent(inout) :: aivals ! array holding sums for ! various statistical output ! !DESCRIPTION: For precipitation rate observations, this routine ! \begin{enumerate} ! \item reads obs assigned to given mpi task (geographic region), ! \item simulates obs from guess, ! \item apply some quality control to obs, ! \item load weight and innovation arrays used in minimization ! \item collects statistics for runtime diagnostic output ! \item writes additional diagnostic information to output file ! \end{enumerate} ! ! !REVISION HISTORY: ! ! 1998-04-30 treadon ! 1999-08-24 derber, j., treadon, r., yang, w., first frozen mpp version ! 2004-05-28 kleist, subroutine call update ! 2004-06-15 treadon - update documentation ! 2004-07-15 todling - protex-compliant prologue; added intent/only's; removed unused mods ! 2004-10-06 parrish - modifications for nonlinear qc ! 2004-10-28 treadon - replace tiny with tiny_r_kind ! 2004-11-22 derber - remove weight, add logical for boundary point ! 2004-12-22 treadon - add outer loop number to name of diagnostic file ! 2005-03-09 parrish - nonlinear qc change to account for inflated obs error ! 2005-03-16 derber - save observation time ! 2005-04-11 treadon - add logical to toggle on/off nonlinear qc code ! 2005-08-05 treadon - increase size of diag_pcp_sat character variable; ! remove mype and add idate to diag_pcp_sat header ! 2005-09-28 derber - modify for surface info from read routines and ! additional input info from pcpinfo ! 2005-10-14 derber - input grid location and fix regional lat/lon ! 2005-10-17 treadon - generalize accessing of elements from obs array ! 2005-11-29 parrish - remove call to deter_sfc_reg (earlier patch for regional mode) ! 2006-01-09 treadon - introduce get_ij ! 2006-02-03 derber - modify for new obs control and obs count ! 2006-03-21 treadon - add option to perturb observation ! 2006-04-12 treadon - replace 1d del,sl with 2d del0,sl0 ! 2006-05-03 derber - modify to write data in jppf chunks ! 2006-06-06 su - move to wgtlim to constants module ! 2006-07-28 derber - modify to use new inner loop obs data structure ! - unify NL qc ! - process one ob at a time ! 2006-07-31 kleist - use separate derivative arrays, change to use ges_ps ! 2006-10-12 treadon - replace virtual temperature with sensible ! 2007-01-19 derber - load log(one+satpcp) instead of satpcp into structure ! 2007-03-19 tremolet - binning of observations ! 2007-03-21 su - remove option to perturb observation ! 2007-06-05 tremolet - add observation diagnostics structure ! 2007-06-08 kleist/treadon - add prefix (task id or path) to diag_pcp_file ! 2008-05-23 safford - rm unused vars and uses ! 2008-12-03 todling - changed handle of tail%time ! 2009-08-19 guo - changed for multi-pass setup with dtime_check(), and ! new arguments init_pass and last_pass. ! - fixed a bug of using "i" instead of "n" when setting the ! vericiation index of obsdiag. ! 2009-12-08 guo - cleaned diag output rewind with open(position='rewind') ! - fixed a bug in diag header output while is not init_pass. ! 2011-05-01 todling - add metguess-bundle; cwmr no longer in guess-grids ! 2013-10-19 todling - metguess now holds background ! tendencies now in bundle ! 2014-12-30 derber - Modify for possibility of not using obsdiag ! 2015-08-31 thomas - Added 4DEnVar to logical to exclude downweighting of ! obs away from center of window ! 2015-10-01 guo - full res obvsr: index to allow redistribution of obsdiags ! 2016-05-18 guo - replaced ob_type with polymorphic obsNode through type casting ! 2016-06-24 guo - fixed the default value of obsdiags(:,:)%tail%luse to luse(n) ! . removed (%dlat,%dlon) debris. ! ! ! !REMARKS: This routine is NOT correctly set up if running ! with the hybrid vertical coordinate, but will ! essentially assign a sigma-like vertical structure ! based on the 'last' observation location ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP; SGI Origin 2000; Compaq/HP ! ! !AUTHOR: ! treadon org: np23 date: 1998-04-30 ! !EOP !------------------------------------------------------------------------- ! Declare local parameters integer(i_kind),parameter:: iint=8 integer(i_kind),parameter:: ncld=1 ! Declare external calls for code analysis external:: stop2 external:: pcp_k ! Declare local variables integer(i_kind) isatid,itime,ilon,ilat,isfcflg,ipcp,isdv integer(i_kind) icnt,ilone,ilate,icnv,itype,iclw,icli integer(i_kind) itim,itimp,istat logical sea logical ssmi,amsu,tmi,stage3,muse logical land,coast,ice logical,allocatable,dimension(:):: luse integer(i_kind),dimension(nobs):: ioid ! initial (pre-distribution) obs ID character(10) filex character(12) string character(128) diag_pcp_file integer(i_kind) km1,mm1,iiflg,iextra,ireal integer(i_kind) ii,i,j,k,m,n,ibin,ioff,ioff0 integer(i_kind) ipt integer(i_kind) nsphys,ixp,iyp,ixx,iyy integer(i_kind) km,ncnt integer(i_kind) kx,jj integer(i_kind) kdt integer(i_kind) ncloud integer(i_kind),dimension(4):: jgrd integer(i_kind),dimension(iint):: idiagbuf integer(i_kind) kbcon,jmin,ktcon,kuo,kbcon4,jmin4,ktcon4 integer(i_kind) ksatid,isflg,istatus real(r_kind) avg,sdv,rterm1,rterm2,rterm real(r_kind) error,a0,a1 real(r_kind) errlog real(r_kind) rdocp,frain,dtp,dtf,sum,sixthpi real(r_kind) drad,vfact,efact,fhour,rtime real(r_kind) xlo,xhi real(r_kind) wgtij real(r_kind) rmiss,pterm,delt,deltp real(r_kind) rmmhr,term,detect_threshold real(r_kind) dzmax,dtdz,dqdz,dudz,dvdz,dcwmdz real(r_kind) cg_pcp,wgross,wnotgross,wgt,arg real(r_single),allocatable,dimension(:):: diagbuf real(r_kind),dimension(nsig):: sl real(r_kind),dimension(nsig+1):: si real(r_kind),dimension(4):: wgrd real(r_kind),dimension(npredp):: pred real(r_kind),dimension(2):: dpcpmag real(r_kind) errf,varinv,error0,dtime,slats,slons,satpcp,pcpsdv, & pcpnum,cenlat,cenlon,satclw,satcli,obserr,pcpnbc,varinvint, & pcpbc,pcpsas,pcplrg4,pcpsas4,xkt4,rbs0,rmask0,psexp4,satcnv, & plon0,plat0,pcplrg,rn1_ad,rn4,psexp,cldwrk,cldwrk4,xkt2 real(r_kind),allocatable,dimension(:,:):: data_p real(r_kind),dimension(nsig)::tsas,qsas,tlrg,qlrg, & t0,q0,cwm0,u0,v0,div0,t1,q1,cwm1,u1,v1,& t0_ad,q0_ad,cwm0_ad,& u0_ad,v0_ad,div0_ad,t1_ad,q1_ad,cwm1_ad,& u1_ad,v1_ad,t4_ad,q4_ad,cwm4_ad,u4_ad,v4_ad,div4_ad, & t_ges,q_ges,u_ges,v_ges,div_ges,cwm_ges,zges,z0,& prsl0,del0,sl0,tsen_ten0,q_ten0,p_ten0 real(r_kind),dimension(nsig+1):: prsi0 real(r_kind),pointer,dimension(:,:,:)::ges_cwmr_im,ges_cwmr_ip real(r_kind),parameter:: zero_7 = 0.7_r_kind real(r_kind),parameter:: r1em6 = 0.000001_r_kind real(r_kind),parameter:: r15 = 15.0_r_kind real(r_kind),parameter:: r25 = 25.0_r_kind ! real(r_kind),parameter:: r40 = 40.0_r_kind real(r_kind),parameter:: r45 = 45.0_r_kind ! real(r_kind),parameter:: r272_16 = 272.16_r_kind real(r_kind),parameter:: r10000 = 10000.0_r_kind character(len=*),parameter:: myname='setuppcp' logical:: in_curbin, in_anybin logical proceed integer(i_kind),dimension(nobs_bins) :: n_alloc integer(i_kind),dimension(nobs_bins) :: m_alloc class(obsNode),pointer:: my_node type(pcpNode),pointer:: my_head type(obs_diag),pointer:: my_diag real(r_kind),allocatable,dimension(:,:,: ) :: ges_ps real(r_kind),allocatable,dimension(:,:,:,:) :: ges_u real(r_kind),allocatable,dimension(:,:,:,:) :: ges_v real(r_kind),allocatable,dimension(:,:,:,:) :: ges_div real(r_kind),allocatable,dimension(:,:,:,:) :: ges_q real(r_kind),pointer,dimension(:,:,:) :: ges_prs_ten=>NULL() real(r_kind),pointer,dimension(:,:,:) :: ges_tv_ten =>NULL() real(r_kind),pointer,dimension(:,:,:) :: ges_q_ten =>NULL() real(r_kind),allocatable,dimension(:,:,:) :: ges_ps_lon real(r_kind),allocatable,dimension(:,:,:) :: ges_ps_lat data rmiss / -999._r_kind / ! Check to see if required guess fields are available call check_vars_(proceed) if(.not.proceed) return ! not all vars available, simply return ! If require guess vars available, extract from bundle ... call init_vars_ !! Verify preconditions if(.not. (drv_initialized.and.tnd_initialized) ) then if(.not.drv_initialized) call perr(myname,'drv_initialized =',drv_initialized) if(.not.tnd_initialized) call perr(myname,'tnd_initialized =',tnd_initialized) call die(myname) endif n_alloc(:)=0 m_alloc(:)=0 !********************************************************************************* ! ONE TIME, INITIAL SETUP PRIOR TO PROCESSING SATELLITE DATA ! ! Initialize variables iiflg = 1 ncloud = ncld nsphys = max(int(two*deltim/dtphys+0.9999_r_kind),1) dtp = two*deltim/nsphys dtf = half*dtp frain = dtf/dtp kdt = 1 fhour = kdt*deltim/r3600 rtime = one/(r3600*fhour) rmmhr = r1000*rtime * r3600 mm1 = mype+1 detect_threshold = one_tenth ! Initialize variables. rdocp = rd/cp ! ! ! Set other constants. sixthpi = half*pi/three ! ! !***** ! INITIALIZE VARIABLE AND DATA ARRAYS !***** ncnt = 0 ! Initialize logical flags for satellite platform ssmi = obstype == 'pcp_ssmi' tmi = obstype == 'pcp_tmi' amsu = obstype == 'pcp_amsu' stage3= obstype == 'pcp_stage3' if(pcp_diagsave)then iextra=0 filex=obstype if(ssmi)filex='pcp_ssmi' write(string,1976)jiter 1976 format('_',i2.2) diag_pcp_file= trim(dirname) // trim(filex) // '_' // trim(dplat(is)) // trim(string) if(init_pass) then open(4,file=trim(diag_pcp_file),form='unformatted',status='unknown',position='rewind') else open(4,file=trim(diag_pcp_file),form='unformatted',status='old',position='append') endif ioff0=22 ireal=ioff0 if (lobsdiagsave) ireal=ireal+4*miter+1 allocate(diagbuf(ireal)) ! ! Initialize/write parameters for precip. diagnostic file on ! first outer iteration. if (init_pass .and. mype==mype_diaghdr(is)) then write(4) isis,dplat(is),obstype,jiter,ianldate,iint,ireal,iextra,ioff0 ! write(6,*)'SETUPPCP: write header record for ',& ! isis,iint,ireal,iextra,' to file ',trim(diag_pcp_file),' ',ianldate endif idiagbuf= 0 diagbuf = rmiss end if ! ! Load data array for current satellite allocate(data_p(nele,nobs)) allocate(luse(nobs)) read(lunin) data_p,luse,ioid ! Index information for data array (see reading routine) isatid = 1 ! index of satellite id itime = 2 ! index of analysis relative obs time ( ilon = 3 ! index of grid relative obs location (x) ilat = 4 ! index of grid relative obs location (y) isfcflg = 5 ! index of surface flag ipcp = 6 ! index of rain rate (mm/hr) isdv = 7 ! index of standard deviation of superob icnt = 8 ! index of number of obs in superob ilone = 9 ! index of earth realtive obs longitude (degrees) ilate = 10 ! index of earth relative obs latitude (degrees) if (tmi) then icnv = 7 iclw = 8 icli = 9 icnt = 10 ilone = 11 ilate = 12 elseif (amsu) then itype = 8 endif rterm1=one/float(nsig) rterm2=one/float(nsig*(nsig-1)) call dtime_setup() do n = 1,nobs ! Extract obs date/time. dtime = data_p(itime,n) call dtime_check(dtime, in_curbin, in_anybin) if(.not.in_anybin) cycle if(in_curbin) then ! ! Determine block of data to use ! ! Initialize variables/arrays. cldwrk = zero psexp = zero xkt2 = zero kbcon = 0 ktcon = 0 jmin = 0 ! INITIAL PROCESSING OF SATELLITE DATA ! ! Extract satellite id ksatid = nint(data_p(isatid,n)) kx = -999 do i = 1,npcptype if (isis == nupcp(i)) kx = i end do if (kx<=0) write(6,*) & 'SETUPPCP: ***WARNING*** Problem with satellite id. ksatid,kx=',& ksatid,kx,isis,(nupcp(i),i=1,npcptype) ! ! Set observation error error0 = varchp(kx) ! ! Extract obs date/time. dtime = data_p(itime,n) ! ! Extract obs lon and lat. slons = data_p(ilon,n) ! grid relative longitude slats = data_p(ilat,n) ! grid relative latitude ! Extract land/sea/ice flag (0=sea, 1=land, 2=ice, 3=snow, 4=mixed) isflg = data_p(isfcflg,n) ! Set logical flags sea = isflg == 0 land = isflg == 1 .or. isflg == 3 ice = isflg == 2 coast = isflg >= 4 ! Check for coastal (transition) points ! ! Extract observations satpcp = data_p(ipcp,n) ! superob rain rate (mm/hr) cenlat = data_p(ilate,n) ! earth relative latitude (degrees) cenlon = data_p(ilone,n) ! earth relative longitude (degrees) satcnv = rmiss if (tmi) then pcpsdv = rmiss satcnv = data_p(icnv,n) satclw = data_p(iclw,n) satcli = data_p(icli,n) pcpnum = data_p(icnt,n) ! number of obs in superob else satclw = rmiss ! cloud liquid water satcli = rmiss ! cloud ice if(.not. amsu)then pcpsdv = data_p(isdv,n) ! standard deviation of superob pcpnum = data_p(icnt,n) ! number of obs in superob else pcpsdv = rmiss pcpnum = rmiss end if endif ! ! Adjust obs error based on precipitation rate. ! ! Bill Olson TMI error model. if (ksatid==211) then if (sea) then a0=0.137_r_kind; a1=0.118_r_kind !tmi,sea else a0=0.335_r_kind; a1=0.045_r_kind !tmi, land endif ! Bob Kuligowski SSMI and AMSU error models elseif (ksatid==264) then if (sea) then a0=0.3835_r_kind; a1=0.1699_r_kind !ssmi, sea (n=50 bin) (scattering) else a0=0.3148_r_kind; a1=0.1781_r_kind !ssmi, land (n=50 bin) endif elseif (ksatid==258) then if (sea) then a0=0.2708_r_kind; a1=0.1997_r_kind !amsu, sea, (n=50 bin) else a0=0.2232_r_kind; a1=0.2619_r_kind !amsu, land (n=50 bin) endif else a0=varchp(kx); a1=zero endif term = log(one+satpcp) errlog = a0 + a1*term obserr = errlog error = obserr errf = error varinv = one/(error*error) ! If negative precipitation, reset to zero if (satpcp < zero) satpcp=zero endif ! (in_curbin) ! Link observation to appropriate observation bin 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 ! Link obs to diagnostics structure if (luse_obsdiag) then if (.not.lobsdiag_allocated) then if (.not.associated(obsdiags(i_pcp_ob_type,ibin)%head)) then obsdiags(i_pcp_ob_type,ibin)%n_alloc = 0 allocate(obsdiags(i_pcp_ob_type,ibin)%head,stat=istat) if (istat/=0) then write(6,*)'setuppcp: failure to allocate obsdiags',istat call stop2(263) end if obsdiags(i_pcp_ob_type,ibin)%tail => obsdiags(i_pcp_ob_type,ibin)%head else allocate(obsdiags(i_pcp_ob_type,ibin)%tail%next,stat=istat) if (istat/=0) then write(6,*)'setuppcp: failure to allocate obsdiags',istat call stop2(264) end if obsdiags(i_pcp_ob_type,ibin)%tail => obsdiags(i_pcp_ob_type,ibin)%tail%next end if obsdiags(i_pcp_ob_type,ibin)%n_alloc = obsdiags(i_pcp_ob_type,ibin)%n_alloc +1 allocate(obsdiags(i_pcp_ob_type,ibin)%tail%muse(miter+1)) allocate(obsdiags(i_pcp_ob_type,ibin)%tail%nldepart(miter+1)) allocate(obsdiags(i_pcp_ob_type,ibin)%tail%tldepart(miter)) allocate(obsdiags(i_pcp_ob_type,ibin)%tail%obssen(miter)) obsdiags(i_pcp_ob_type,ibin)%tail%indxglb=ioid(n) obsdiags(i_pcp_ob_type,ibin)%tail%nchnperobs=-99999 obsdiags(i_pcp_ob_type,ibin)%tail%luse=luse(n) obsdiags(i_pcp_ob_type,ibin)%tail%muse(:)=.false. obsdiags(i_pcp_ob_type,ibin)%tail%nldepart(:)=-huge(zero) obsdiags(i_pcp_ob_type,ibin)%tail%tldepart(:)=zero obsdiags(i_pcp_ob_type,ibin)%tail%wgtjo=-huge(zero) obsdiags(i_pcp_ob_type,ibin)%tail%obssen(:)=zero n_alloc(ibin) = n_alloc(ibin) +1 my_diag => obsdiags(i_pcp_ob_type,ibin)%tail my_diag%idv = is my_diag%iob = ioid(n) my_diag%ich = 1 my_diag%elat= data_p(ilate,n) my_diag%elon= data_p(ilone,n) else if (.not.associated(obsdiags(i_pcp_ob_type,ibin)%tail)) then obsdiags(i_pcp_ob_type,ibin)%tail => obsdiags(i_pcp_ob_type,ibin)%head else obsdiags(i_pcp_ob_type,ibin)%tail => obsdiags(i_pcp_ob_type,ibin)%tail%next end if if (.not.associated(obsdiags(i_pcp_ob_type,ibin)%tail)) then call die(myname,'.not.associated(obsdiags(i_pcp_ob_type,ibin)%tail)') end if if (obsdiags(i_pcp_ob_type,ibin)%tail%indxglb/=ioid(n)) then write(6,*)'setuppcp: index error' call stop2(265) end if endif endif if(.not.in_curbin) cycle ! ! Do some prelimiary qc of the data. Check for points covered ! with ice/snow, land points, and coastal points. ! ! ! !***** ! CALL Precipitation (FORWARD) MODEL !***** ! ! ! Call precipitation model. Set ajm forcing ! for precipitation rate to unity. ! Initialize variables and arrays km = nsig do k = 1,nsig t0_ad(k) = zero q0_ad(k) = zero u0_ad(k) = zero v0_ad(k) = zero div0_ad(k) = zero cwm0_ad(k) = zero t_ges(k) = zero q_ges(k) = zero u_ges(k) = zero v_ges(k) = zero cwm_ges(k) = zero div_ges(k) = zero zges(k) = zero end do pcpnbc = zero pcpsas = zero pcplrg = zero if(dtime > hrdifsig(1) .and. dtime < hrdifsig(nfldsig))then do jj=1,nfldsig-1 if(dtime > hrdifsig(jj) .and. dtime <= hrdifsig(jj+1))then itim=jj itimp=jj+1 delt=((hrdifsig(jj+1)-dtime)/(hrdifsig(jj+1)-hrdifsig(jj))) end if end do else if(dtime <=hrdifsig(1))then itim=1 itimp=1 delt=one else itim=nfldsig itimp=nfldsig delt=one end if deltp=one-delt ! Get pointer to could water mixing ratio call gsi_bundlegetpointer (gsi_metguess_bundle(itim), 'cw',ges_cwmr_im,istatus) if (istatus/=0) call die('setuppcp','cannot get pointer to cwmr(itim), istatus =',istatus) call gsi_bundlegetpointer (gsi_metguess_bundle(itimp),'cw',ges_cwmr_ip,istatus) if (istatus/=0) call die('setuppcp','cannot get pointer to cwmr(itimp), istatus =',istatus) ! Set and save spatial interpolation indices and weights. call get_ij(mm1,slats,slons,jgrd,wgrd,jjlat=ixx,jjlon=iyy) ixp=ixx+1 iyp=iyy+1 ! Loop over surrounding analysis gridpoints do ipt=1,4 ! Set (i,j) index for ipt-th point for observation n. if(ipt == 1)then i=ixx j=iyy else if (ipt == 2)then i=ixp j=iyy else if (ipt == 3)then i=ixx j=iyp else i=ixp j=iyp end if ! Load arrays used by forward/adjoint model xkt4 = xkt2d(i,j) ii = i+istart(mm1)-2 ii = max(1,min(ii,nlat)) rbs0 = rbs2(ii) rmask0= isli2(i,j) ! Perform interpolation. pterm = delt*ges_ps(i,j,itim) + deltp*ges_ps(i,j,itimp) psexp4 = pterm plon0 = delt*ges_ps_lon(i,j,itim) + deltp*ges_ps_lon(i,j,itimp) plat0 = delt*ges_ps_lat(i,j,itim) + deltp*ges_ps_lat(i,j,itimp) do k=1,nsig t0(k) = delt *ges_tsen(i,j,k,itim) + deltp*ges_tsen(i,j,k,itimp) q0(k) = delt *ges_q(i,j,k,itim) + deltp*ges_q(i,j,k,itimp) u0(k) = delt *ges_u(i,j,k,itim) + deltp*ges_u(i,j,k,itimp) v0(k) = delt *ges_v(i,j,k,itim) + deltp*ges_v(i,j,k,itimp) div0(k) = delt *ges_div(i,j,k,itim) + deltp*ges_div(i,j,k,itimp) cwm0(k) = delt *ges_cwmr_im(i,j,k) + deltp*ges_cwmr_ip(i,j,k) z0(k) = delt *geop_hgtl(i,j,k,itim) + deltp*geop_hgtl(i,j,k,itimp) prsl0(k) = delt *ges_prsl(i,j,k,itim) + deltp*ges_prsl(i,j,k,itimp) prsi0(k) = delt *ges_prsi(i,j,k,itim) + deltp*ges_prsi(i,j,k,itimp) q_ten0(k) = ges_q_ten(i,j,k) tsen_ten0(k) = (ges_tv_ten(i,j,k) - fv*t0(k)*ges_q_ten(i,j,k))/(one+fv*q0(k)) ! NOTE: k=1 is surface pressure tendency p_ten0(k) = ges_prs_ten(i,j,k+1) end do k=nsig+1 prsi0(k) = delt *ges_prsi(i,j,k,itim ) + deltp*ges_prsi(i,j,k,itimp) ! Create sigma-like coefficients (used in pcp_k) do k=1,nsig sl(k) = prsl0(k) / psexp4 si(k) = prsi0(k) / psexp4 end do k=nsig+1 si(k) = prsi0(k) / psexp4 do k=1,nsig del0(k) = si(k)-si(k+1) sl0(k) = sl(k) end do ! Initialize arrays below prior to each call to forward/adjoint model do k = 1,nsig t4_ad(k) = zero q4_ad(k) = zero u4_ad(k) = zero v4_ad(k) = zero div4_ad(k) = zero cwm4_ad(k) = zero t1_ad(k) = zero q1_ad(k) = zero u1_ad(k) = zero v1_ad(k) = zero cwm1_ad(k) = zero t1(k) = zero q1(k) = zero u1(k) = zero v1(k) = zero cwm1(k) = zero end do rn4 = zero rn1_ad = one pcpsas4 = zero pcplrg4 = zero ! Call forward/adjoint model for current analysis gridpoint call pcp_k(km,dtp,del0,sl0,rbs0, & rmask0,xkt4,ncloud,frain,rmmhr,& psexp4,plat0,plon0,& t0,q0,u0,v0,div0,cwm0,& tsen_ten0,q_ten0,p_ten0,& tsas,qsas,pcpsas4,cldwrk4,kbcon4,ktcon4,jmin4,kuo, & tlrg,qlrg,pcplrg4, & t1,q1,cwm1,u1,v1,rn4,& t4_ad,q4_ad,cwm4_ad,u4_ad,v4_ad,div4_ad,& t1_ad,q1_ad,cwm1_ad,u1_ad,v1_ad,rn1_ad) ! Based on bilinear interpolation weights, accumulate rain rates ! and sensitivity over four surrounding analysis gridpoints. wgtij=wgrd(ipt) pcpsas = pcpsas + wgtij*pcpsas4 pcplrg = pcplrg + wgtij*pcplrg4 pcpnbc = pcpnbc + wgtij*rn4 kbcon = kbcon + wgtij*kbcon4 jmin = jmin + wgtij*jmin4 ktcon = ktcon + wgtij*ktcon4 psexp = psexp + wgtij*psexp4 cldwrk = cldwrk + wgtij*cldwrk4 xkt2 = xkt2 + wgtij*xkt4 do k=1,nsig t0_ad(k) = t0_ad(k) + wgtij*t4_ad(k) q0_ad(k) = q0_ad(k) + wgtij*q4_ad(k) cwm0_ad(k) = cwm0_ad(k) + wgtij*cwm4_ad(k) u0_ad(k) = u0_ad(k) + wgtij*u4_ad(k) v0_ad(k) = v0_ad(k) + wgtij*v4_ad(k) div0_ad(k) = div0_ad(k) + wgtij*div4_ad(k) t_ges(k) = t_ges(k) + wgtij*t0(k) q_ges(k) = q_ges(k) + wgtij*q0(k) u_ges(k) = u_ges(k) + wgtij*u0(k) v_ges(k) = v_ges(k) + wgtij*v0(k) cwm_ges(k) = cwm_ges(k) + wgtij*cwm0(k) div_ges(k) = div_ges(k) + wgtij*div0(k) zges(k) = zges(k) + wgtij*z0(k) end do ! End of loop over surrounding gridpoints end do ! If magnitude of sensitivity vector is zero, set inverse observation ! error to zero since we don't have gradient information for this ! observation. For now, ommit divergence sensitivity since it is ! not used in pcgsoi. Likewise, if the gradient is too large, do not ! use this observation (this check rarely catches any points). dpcpmag=zero dzmax =-9.99e11_r_kind do k=1,nsig dpcpmag(1) = dpcpmag(1) + t0_ad(k) dpcpmag(1) = dpcpmag(1) + q0_ad(k) dpcpmag(1) = dpcpmag(1) + u0_ad(k) dpcpmag(1) = dpcpmag(1) + v0_ad(k) dpcpmag(1) = dpcpmag(1) + cwm0_ad(k) dpcpmag(2) = dpcpmag(2) + t0_ad(k)*t0_ad(k) dpcpmag(2) = dpcpmag(2) + q0_ad(k)*q0_ad(k) dpcpmag(2) = dpcpmag(2) + u0_ad(k)*u0_ad(k) dpcpmag(2) = dpcpmag(2) + v0_ad(k)*v0_ad(k) dpcpmag(2) = dpcpmag(2) + cwm0_ad(k)*cwm0_ad(k) km1=max(1,k-1) term=zges(k)-zges(km1) ! rterm=one/term rterm = zero if (k>km1) rterm=one/term dtdz = abs((t0_ad(k)-t0_ad(km1))*rterm) dqdz = abs((q0_ad(k)-q0_ad(km1))*rterm) dudz = abs((u0_ad(k)-u0_ad(km1))*rterm) dvdz = abs((v0_ad(k)-v0_ad(km1))*rterm) dcwmdz = abs((cwm0_ad(k)-cwm0_ad(km1))*rterm) dzmax = max(dtdz,dqdz,dudz,dvdz,dcwmdz,dzmax) end do avg = dpcpmag(1)*rterm1 term= (nsig*dpcpmag(2)-dpcpmag(1)**2)*rterm2 sdv = zero if (term>zero) sdv=sqrt(term) dpcpmag(2)=sqrt(dpcpmag(2)) if (dpcpmag(2)r10000) then if(luse(n))aivals(27,is) = aivals(27,is) + one varinv=zero elseif (abs(avg)>zero) then if (abs(sdv/avg)>r25) then if(luse(n))aivals(28,is) = aivals(28,is) + one varinv = zero endif endif if (abs(dzmax)>one) then if(luse(n))aivals(29,is) = aivals(29,is) + one varinv = zero endif if(luse(n))aivals(1,is) = aivals(1,is) + one aivals(2,is) = aivals(2,is) + one if (varinv>r1em6 .and. luse(n)) aivals(3,is) = aivals(3,is) + one ! ! !***** ! COMPUTE AND APPLY BIAS CORRECTION TO SIMULATED VALUES !***** ! ! Prepare for application of bias correction to simulated values. ! ! Construct predictors for bias correction. do k=1,npredp pred(k) = zero end do ! ! Apply bias correction to simulated values. sum = zero do j = 1,npredp sum = sum + predxp(j,kx)*pred(j)*ibias(kx) end do pcpbc = pcpnbc + sum pcpbc = max(pcpbc,zero) ! ! !*** ! QC OBSERVATIONS BASED ON VARIOUS CRITERIA !*** ! ! QC SSMI, TMI, AMSU, and STAGE3 precipitation. if (ssmi .or. tmi .or. amsu .or. stage3) then ! efact = one vfact = quarter ! Reduce qc bounds in extratropics. Between 45 ! and 60 degrees latitude linearly reduce error ! bound. Poleward of 60 degrees latitude do not use if ( abs(cenlat) > r45 ) then if(luse(n))aivals(4,is) = aivals(4,is) + one term = (r60-abs(cenlat))/r15 term = max(zero,term) term = min(term,one) if (abs(cenlat) > r60) term=zero efact = term*efact endif ! ! If land obs, reduce error bounds and weight if (land) then efact = zero_7*efact vfact = zero_7*vfact if(luse(n))aivals(5,is) = aivals(5,is) + one end if ! ! Toss obs if over snow or ice if (ice) then efact = zero vfact = zero if(luse(n))aivals(6,is) = aivals(6,is) + one end if ! If coastal obs, reduce error bounds and weight if (coast) then efact = half*efact vfact = half*vfact if(luse(n))aivals(7,is) = aivals(7,is) + one end if ! ! Reduce weight and tighten obs-ges bound as function of obs time. if(.not.(l4dvar.or.l4densvar))then term = cos(sixthpi*(dtime-half*hr_obsbin)) term = term**3 term = max(zero,min(term,one)) efact = term*efact vfact = term*vfact end if ! If the rain rate is greater than zero but less than the ! detection threshold (currently take to be 0.1 mm/hr), then ! set the observation error to zero. ! That is, do not use this observation in the analysis. if ( (satpcp>zero) .and. (satpcp r1em6) then if (drad < xlo .or. drad > xhi) then varinv = zero if(luse(n))aivals(12,is) = aivals(12,is) + one endif endif muse= (varinv>r1em6.and.iusep(kx)>=1) if (nobskeep>0.and.luse_obsdiag) muse=obsdiags(i_pcp_ob_type,ibin)%tail%muse(nobskeep) if (luse_obsdiag) then obsdiags(i_pcp_ob_type,ibin)%tail%muse(jiter)=muse obsdiags(i_pcp_ob_type,ibin)%tail%nldepart(jiter)= drad obsdiags(i_pcp_ob_type,ibin)%tail%wgtjo= varinv endif ! !*** ! CONSTRUCT SENSITIVITY VECTORS AND SAVE INFO FOR INNER ITERATION !*** ! ! Optionally, set parameter iusep to turn off data. if (muse) then ncnt = ncnt+1 allocate(my_head) m_alloc(ibin) = m_alloc(ibin) +1 my_node => my_head ! this is a workaround call obsLList_appendNode(pcphead(ibin),my_node) my_node => null() my_head%idv = is my_head%iob = ioid(n) my_head%elat= data_p(ilate,n) my_head%elon= data_p(ilone,n) allocate(my_head%predp(npredp),my_head%dpcp_dvar(nsig5), & stat=istat) if(istat /= 0)write(6,*)' failure to write pcphead arrays ' do m=1,4 my_head%ij(m)=jgrd(m) my_head%wij(m)=wgrd(m) end do my_head%err2=one/error0**2 my_head%raterr2=error0**2*varinv do m=1,npredp my_head%predp(m)=pred(m) end do my_head%obs=log(one+satpcp) my_head%ges=pcpbc my_head%time=dtime do m=1,nsig my_head%dpcp_dvar(m)=t0_ad(m) my_head%dpcp_dvar(m+nsig)=q0_ad(m) my_head%dpcp_dvar(m+nsig2)=u0_ad(m) my_head%dpcp_dvar(m+nsig3)=v0_ad(m) my_head%dpcp_dvar(m+nsig4)=cwm0_ad(m) ! my_head%dpcp_dvar(m+nsig5)=div0_ad(m) !need to increase size end do my_head%icxp=kx my_head%luse=luse(n) if (luse_obsdiag) then my_head%diags => obsdiags(i_pcp_ob_type,ibin)%tail my_diag => my_head%diags if(my_head%idv /= my_diag%idv .or. & my_head%iob /= my_diag%iob ) then call perr(myname,'mismatching %[head,diags]%(idv,iob,ibin) =', & (/is,ioid(n),ibin/)) call perr(myname,'my_head%(idv,iob) =',(/my_head%idv,my_head%iob/)) call perr(myname,'my_diag%(idv,iob) =',(/my_diag%idv,my_diag%iob/)) call die(myname) endif endif my_head => null() end if if(luse(n))then ! ! Save data in diagnostic arrays. if (pcp_diagsave) then idiagbuf(1) = ksatid idiagbuf(2) = 0 ! spare idiagbuf(3) = 0 ! spare idiagbuf(4) = isflg idiagbuf(5) = 0 ! spare idiagbuf(6) = kbcon idiagbuf(7) = ktcon idiagbuf(8) = jmin diagbuf(1) = cenlat diagbuf(2) = cenlon diagbuf(3) = dtime-time_offset diagbuf(4) = pcpnum diagbuf(5) = satpcp diagbuf(6) = pcpsdv diagbuf(7) = satcnv diagbuf(8) = satclw diagbuf(9) = satcli diagbuf(10) = one/(error*error) diagbuf(11) = error if (kbcon>=1 .and. kbcon<=nsig) diagbuf(12) = r10*psexp * sl0(kbcon) if (ktcon>=1 .and. ktcon<=nsig) diagbuf(13) = r10*psexp * sl0(ktcon) if (jmin>=1 .and. jmin<=nsig) diagbuf(14) = r10*psexp * sl0(jmin) diagbuf(15) = cldwrk diagbuf(16) = pcpsas diagbuf(17) = pcpnbc diagbuf(18) = pcpbc diagbuf(19) = varinvint diagbuf(20) = varinv diagbuf(21) = errf diagbuf(22) = xkt2 ioff=ioff0 if (lobsdiagsave) then do jj=1,miter ioff=ioff+1 if (obsdiags(i_pcp_ob_type,ibin)%tail%muse(jj)) then diagbuf(ioff) = one else diagbuf(ioff) = -one endif enddo do jj=1,miter+1 ioff=ioff+1 diagbuf(ioff) = obsdiags(i_pcp_ob_type,ibin)%tail%nldepart(jj) enddo do jj=1,miter ioff=ioff+1 diagbuf(ioff) = obsdiags(i_pcp_ob_type,ibin)%tail%tldepart(jj) enddo do jj=1,miter ioff=ioff+1 diagbuf(ioff) = obsdiags(i_pcp_ob_type,ibin)%tail%obssen(jj) enddo endif ! Write diagnostics to output file. write(4) mype,(idiagbuf(m),m=1,iint),(diagbuf(m),m=1,ireal) endif ! ! ! !****** ! DIAGNOSTIC CALCULATIONS. !****** ! ! ! Accumulate data for diagnostic print. if (varinv > r1em6) then drad = log(one+satpcp) - log(one+pcpbc) obserr = sqrt(one/varinv) if (iusep(kx)<1) varinv=zero aivals(11,is) = aivals(11,is) + one aivals(13,is) = aivals(13,is) + drad !bias aivals(14,is) = aivals(14,is) + drad*drad !bias**2 aivals(15,is) = aivals(15,is) + drad*drad*varinv !pen aivals(16,is) = aivals(16,is) + obserr !obs error error aivals(21,is) = aivals(21,is) + satpcp aivals(22,is) = aivals(22,is) + pcpnbc aivals(23,is) = aivals(23,is) + pcpbc aivals(24,is) = aivals(24,is) + pcpsas aivals(25,is) = aivals(25,is) + one term = (drad/error0)**2 if (pg_pcp(kx) > tiny_r_kind .and. error > tiny_r_kind) then arg = exp(-half*term) wnotgross= one-pg_pcp(kx) cg_pcp=b_pcp(kx)*error wgross = cg_term*pg_pcp(kx)/(cg_pcp*wnotgross) term = log((arg+wgross)/(one+wgross)) wgt = one-wgross/(arg+wgross) if(wgt < wgtlim) aivals(40,is)=aivals(40,is)+one endif aivals(39,is) = aivals(39,is) +(error0**2)*varinv*term else aivals(31,is) = aivals(31,is) + satpcp aivals(32,is) = aivals(32,is) + pcpnbc aivals(33,is) = aivals(33,is) + pcpbc aivals(34,is) = aivals(34,is) + pcpsas aivals(35,is) = aivals(35,is) + one endif end if ! End of loop over observations end do ! Release memory of local guess arrays call final_vars_ ! ! Deallocate arrays ! deallocate(data_p) deallocate(luse) if (pcp_diagsave) then close(4) deallocate(diagbuf) call dtime_show(myname,'diagsave:pcp',i_pcp_ob_type) endif ! 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 call gsi_metguess_get ('var::u' , ivar, istatus ) proceed=proceed.and.ivar>0 call gsi_metguess_get ('var::v' , ivar, istatus ) proceed=proceed.and.ivar>0 call gsi_metguess_get ('var::div', ivar, istatus ) proceed=proceed.and.ivar>0 call gsi_metguess_get ('var::q', ivar, istatus ) proceed=proceed.and.ivar>0 end subroutine check_vars_ subroutine init_vars_ 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 u ... varname='u' call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) if (istatus==0) then if(allocated(ges_u))then write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' call stop2(999) endif allocate(ges_u(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) ges_u(:,:,:,1)=rank3 do ifld=2,nfldsig call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) ges_u(:,:,:,ifld)=rank3 enddo else write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus call stop2(999) endif ! get v ... varname='v' call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) if (istatus==0) then if(allocated(ges_v))then write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' call stop2(999) endif allocate(ges_v(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) ges_v(:,:,:,1)=rank3 do ifld=2,nfldsig call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) ges_v(:,:,:,ifld)=rank3 enddo else write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus call stop2(999) endif ! get div ... varname='div' call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) if (istatus==0) then if(allocated(ges_div))then write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' call stop2(999) endif allocate(ges_div(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) ges_div(:,:,:,1)=rank3 do ifld=2,nfldsig call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) ges_div(:,:,:,ifld)=rank3 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 else write(6,*) trim(myname), ': inconsistent vector sizes (nfldsig,size(metguess_bundle) ',& nfldsig,size(gsi_metguess_bundle) call stop2(999) endif ! extract tendencies from tendency bundle if(tnd_initialized) then ! get prs_ten ... varname='prse' call gsi_bundlegetpointer(gsi_tendency_bundle,trim(varname),ges_prs_ten,istatus) if (istatus/=0) then write(6,*) trim(myname),': ', trim(varname), ' not found in tend bundle, ier= ',istatus call stop2(999) endif ! get tv_ten ... varname='tv' call gsi_bundlegetpointer(gsi_tendency_bundle,trim(varname),ges_tv_ten,istatus) if (istatus/=0) then write(6,*) trim(myname),': ', trim(varname), ' not found in tend bundle, ier= ',istatus call stop2(999) endif ! get q_ten ... varname='q' call gsi_bundlegetpointer(gsi_tendency_bundle,trim(varname),ges_q_ten,istatus) if (istatus/=0) then write(6,*) trim(myname),': ', trim(varname), ' not found in tend bundle, ier= ',istatus call stop2(999) endif else write(6,*) trim(myname), ': tendency bundle not properly allocated, aborting ... ' call stop2(999) endif ! extract derivatives from derivatives bundle if(drv_initialized) then ! get ps_lon ... varname='ps' call gsi_bundlegetpointer(gsi_xderivative_bundle(1),trim(varname),rank2,istatus) if (istatus==0) then if(allocated(ges_ps_lon))then write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' call stop2(999) endif allocate(ges_ps_lon(size(rank2,1),size(rank2,2),nfldsig)) ges_ps_lon(:,:,1)=rank2 do ifld=2,nfldsig call gsi_bundlegetpointer(gsi_xderivative_bundle(ifld),trim(varname),rank2,istatus) ges_ps_lon(:,:,ifld)=rank2 enddo else write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus call stop2(999) endif ! get ps_lat ... varname='ps' call gsi_bundlegetpointer(gsi_yderivative_bundle(1),trim(varname),rank2,istatus) if (istatus==0) then if(allocated(ges_ps_lat))then write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' call stop2(999) endif allocate(ges_ps_lat(size(rank2,1),size(rank2,2),nfldsig)) ges_ps_lat(:,:,1)=rank2 do ifld=2,nfldsig call gsi_bundlegetpointer(gsi_xderivative_bundle(ifld),trim(varname),rank2,istatus) ges_ps_lat(:,:,ifld)=rank2 enddo else write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus call stop2(999) endif endif end subroutine init_vars_ subroutine init_netcdf_diag_ end subroutine init_netcdf_diag_ subroutine contents_binary_diag_ end subroutine contents_binary_diag_ subroutine contents_netcdf_diag_ ! Observation class character(7),parameter :: obsclass = ' pcp' end subroutine contents_netcdf_diag_ subroutine final_vars_ if(associated(ges_q_ten )) nullify(ges_q_ten ) if(associated(ges_tv_ten )) nullify(ges_tv_ten ) if(associated(ges_prs_ten)) nullify(ges_prs_ten) ! if(allocated(ges_q )) deallocate(ges_q ) if(allocated(ges_div))deallocate(ges_div) if(allocated(ges_v )) deallocate(ges_v ) if(allocated(ges_u )) deallocate(ges_u ) if(allocated(ges_ps)) deallocate(ges_ps) end subroutine final_vars_ end subroutine setuppcp