subroutine setuprhsall(ndata,mype,init_pass,last_pass) !$$$ subprogram documentation block ! . . . . ! subprogram: setuprhsall sets up rhs of oi ! prgmmr: derber org: np23 date: 2003-05-22 ! ! abstract: This routine sets up the right hand side (rhs) of the ! analysis equation. Functions performed in this routine ! include: ! a) calculate increments between current solutions and obs, ! b) generate statistical summaries of quality control and innovations, ! c) generate diagnostic files (optional), and ! d) prepare/save information for use in inner minimization loop ! ! program history log: ! 2003-05-22 derber ! 2003-12-23 kleist - ozone calculation modified to use guess pressure ! 2004-06-17 treadon - update documentation ! 2004-07-23 derber - modify to include conventional sst ! 2004-07-29 treadon - add only to module use, add intent in/out ! 2004-10-06 parrish - increase dimension of work arrays for nonlin qc ! 2004-12-08 xu li - replace local logical flag retrieval with that in radinfo ! 2004-12-22 treadon - restructure code to compute and write out ! innovation information on select outer iterations ! 2005-01-20 okamoto - add ssmi/amsre/ssmis ! 2005-03-30 lpchang - statsoz call was passing ozmz var unnecessarily ! 2005-04-18 treadon - deallocate fbias ! 2005-05-27 derber - level output change ! 2005-07-06 derber - include mhs and hirs/4 ! 2005-06-14 wu - add OMI oz ! 2005-07-27 derber - add print of monitoring and reject data ! 2005-09-28 derber - simplify data file info handling ! 2005-10-20 kazumori - modify for real AMSR-E data process ! 2005-12-01 cucurull - add GPS bending angle ! 2005-12-21 treadon - modify processing of GPS data ! 2006-01-09 derber - move create/destroy array, compute_derived, q_diag ! from glbsoi outer loop into this routine ! 2006-01-12 treadon - add channelinfo ! 2006-02-03 derber - modify for new obs control and obs count- clean up! ! 2006-02-24 derber - modify to take advantage of convinfo module ! 2006-03-21 treadon - add code to generate optional observation perturbations ! 2006-07-28 derber - modify code for new inner loop obs data structure ! 2006-07-29 treadon - remove create_atm_grids and destroy_atm_grids ! 2006-07-31 kleist - change call to atm arrays routines ! 2007-02-21 sienkiewicz - add MLS ozone changes ! 2007-03-01 treadon - add toss_gps and toss_gps_sub ! 2007-03-10 su - move the observation perturbation to each setup routine ! 2007-03-19 tremolet - Jo table ! 2007-06-05 tremolet - add observation diagnostics structure ! 2007-06-08 kleist/treadon - add prefix (task id or path) to diag_conv_file ! 2007-07-09 tremolet - observation sensitivity ! 2007-06-20 cucurull - changes related to gps diagnostics ! 2007-06-29 jung - change channelinfo to array ! 2007-09-30 todling - add timer ! 2007-10-03 todling - add observer split option ! 2007-12-15 todling - add prefix to diag filenames ! 2008-03-28 wu - move optional randon seed for perturb_obs to read_obs ! 2008-04-14 treadon - remove super_gps, toss_gps (moved into genstats_gps) ! 2008-05-23 safford - rm unused vars and uses ! 2008-12-08 todling - move 3dprs/geop-hght calculation from compute_derivate into here ! 2009-01-17 todling - update interface to intjo ! 2009-03-05 meunier - add call to lagragean operator ! 2009-08-19 guo - moved all rhs related statistics variables to m_rhs ! for multi-pass setuprhsall(); ! - added control arguments init_pass and last_pass for ! multi-pass setuprhsall(). ! 2009-09-14 guo - invoked compute_derived() even under lobserver. This is ! the right way to do it. It trigged moving of statments ! from glbsoi() to observer_init(). ! - cleaned up redandent calls to setupyobs() and inquire_obsdiags(). ! 2009-10-22 shen - add high_gps and high_gps_sub ! 2009-12-08 guo - fixed diag_conv output rewind while is not init_pass, with open(position='rewind') ! 2010-04-09 cucurull - remove high_gps and high_gps_sub ! 2010-04-01 tangborn - start adding call for carbon monoxide data. ! 2010-04-28 zhu - add ostats and rstats for additional precoditioner ! 2010-05-28 todling - obtain variable id's on the fly (add getindex) ! 2010-10-14 pagowski - added pm2_5 conventional obs ! 2010-10-20 hclin - added aod ! 2011-02-16 zhu - add gust,vis,pblh ! 2011-04-07 todling - newpc4pred now in radinfo ! 2011-09-17 todling - automatic sizes definition for mpi-reduce calls ! 2012-01-11 Hu - add load_gsdgeop_hgt to compute 2d subdomain pbl heights from the guess fields ! 2012-04-08 Hu - add code to skip the observations that are not used in minimization ! 2013-02-22 Carley - Add call to load_gsdgeop_hgt for NMMB/WRF-NMM if using ! PBL pseudo obs ! 2013-10-19 todling - metguess now holds background ! 2013-05-24 zhu - add ostats_t and rstats_t for aircraft temperature bias correction ! 2014-03-19 pondeca - add wspd10m ! 2014-04-10 pondeca - add td2m,mxtm,mitm,pmsl ! 2014-05-07 pondeca - add howv ! 2014-12-30 derber - Modify for possibility of not using obsdiag ! 2014-0?-16 carley/zhu - add tcamt and lcbas ! 2015-07-10 pondeca - add cldch ! 2015-10-01 guo - full res obvsr: index to allow redistribution of obsdiags ! 2016-05-05 pondeca - add uwnd10m, vwund10m ! 2017-05-12 Y. Wang and X. Wang - add dbz for reflectivity DA. POC: xuguang.wang@ou.edu ! 2018-02-15 wu - add code for fv3_regional ! 2018-01-01 Apodaca - add GOES/GLM lightning ! 2019-03-15 Ladwig - add option for cloud analysis in observer ! 2019-03-28 Ladwig - add metar cloud obs as pseudo water vapor in var analysis ! ! input argument list: ! ndata(*,1)- number of prefiles retained for further processing ! ndata(*,2)- number of observations read ! ndata(*,3)- number of observations keep after read ! mype - mpi task id ! ! output argument list: ! ! ! comments: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,i_kind,r_quad,r_single use constants, only: zero,one,fv,zero_quad use guess_grids, only: load_prsges,load_geop_hgt,load_gsdpbl_hgt use guess_grids, only: ges_tsen,nfldsig ! use obsmod, only: mype_diaghdr use obsmod, only: nsat1,iadate,nobs_type,obscounts,& ndat,obs_setup,& dirname,write_diag,ditype,obsdiags,lobserver,& destroyobs,inquire_obsdiags,lobskeep,nobskeep,lobsdiag_allocated, & luse_obsdiag use obsmod, only: lobsdiagsave,diag_radardbz use obsmod, only: binary_diag use obs_sensitivity, only: lobsensfc, lsensrecompute use radinfo, only: newpc4pred use radinfo, only: mype_rad,diag_rad,jpch_rad,retrieval,fbias,npred,ostats,rstats ! use aircraftinfo, only: aircraft_t_bc_pof,aircraft_t_bc,ostats_t,rstats_t,npredt,max_tail use aircraftinfo, only: aircraft_t_bc_pof,aircraft_t_bc,ostats_t,rstats_t,npredt,ntail use pcpinfo, only: diag_pcp use ozinfo, only: diag_ozone,mype_oz,jpch_oz,ihave_oz use lightinfo, only: mype_light,diag_light use coinfo, only: diag_co,mype_co,jpch_co,ihave_co use mpimod, only: ierror,mpi_comm_world,mpi_rtype,mpi_sum use gridmod, only: nsig,twodvar_regional,wrf_mass_regional,nems_nmmb_regional use gridmod, only: cmaq_regional,fv3_regional use gsi_4dvar, only: nobs_bins,l4dvar use gsi_4dvar, only: mPEs_observer use jfunc, only: jiter,jiterstart,miter,first,last use qcmod, only: npres_print use convinfo, only: nconvtype,diag_conv use timermod, only: timer_ini,timer_fnl use lag_fields, only: lag_presetup,lag_state_write,lag_state_read,lag_destroy_uv use state_vectors, only: svars2d use mpeu_util, only: getindex use mpl_allreducemod, only: mpl_allreduce use aeroinfo, only: diag_aero use berror, only: reset_predictors_var use rapidrefresh_cldsurf_mod, only: l_PBL_pseudo_SurfobsT,l_PBL_pseudo_SurfobsQ,& l_PBL_pseudo_SurfobsUV,i_gsdcldanal_type,& i_cloud_q_innovation use m_rhs, only: rhs_alloc use m_rhs, only: rhs_dealloc use m_rhs, only: rhs_allocated use m_rhs, only: awork => rhs_awork use m_rhs, only: bwork => rhs_bwork use m_rhs, only: aivals => rhs_aivals use m_rhs, only: stats => rhs_stats use m_rhs, only: stats_co => rhs_stats_co use m_rhs, only: stats_oz => rhs_stats_oz use m_rhs, only: toss_gps_sub => rhs_toss_gps use m_gpsStats, only: gpsStats_genstats ! was genstats_gps() use m_gpsStats, only: gpsStats_destroy ! was done by genstats_gps() use gsi_bundlemod, only: GSI_BundleGetPointer use gsi_metguess_mod, only: GSI_MetGuess_Bundle use m_obsdiags, only: obsdiags_reset use m_obsdiags, only: obsdiags_read use m_obsdiags, only: obsdiags_sort use m_obsdiags, only: obsdiags_write use mpeu_util, only: die,warn,perr use mpeu_util, only: basename implicit none ! Declare passed variables integer(i_kind) ,intent(in ) :: mype integer(i_kind),dimension(ndat,3),intent(in ) :: ndata logical ,intent(in ) :: init_pass, last_pass ! state of "setup" processing ! Declare external calls for code analysis external:: compute_derived external:: evaljo !external:: genstats_gps external:: mpi_allreduce external:: mpi_finalize external:: mpi_reduce !external:: read_obsdiags external:: setupaod external:: setupbend external:: setupdw external:: setuplag external:: setupozlay external:: setupozlev external:: setuppcp external:: setupps external:: setuppw external:: setupq external:: setupcldtot external:: setuprad external:: setupdbz external:: setupref external:: setuprw external:: setupspd external:: setupsst external:: setupt external:: setuptcp external:: setupw external:: setupgust external:: setupvis external:: setuppblh external:: setupwspd10m external:: setuptd2m external:: setupmxtm external:: setupmitm external:: setuppmsl external:: setuphowv external:: setuptcamt external:: setuplcbas external:: setupcldch external:: setupuwnd10m external:: setupvwnd10m external:: setupswcp external:: setuplwcp external:: setuplight external:: statsconv external:: statsoz external:: statspcp external:: statsrad external:: statslight external:: stop2 external:: w3tage ! Delcare local variables logical rad_diagsave,ozone_diagsave,pcp_diagsave,conv_diagsave,llouter,getodiag,co_diagsave logical aero_diagsave logical light_diagsave logical radardbz_diagsave character(80):: string character(10)::obstype character(20)::isis character(128):: diag_conv_file character(128):: diag_light_file character(len=12) :: clfile integer(i_kind) lunin,nobs,nchanl,nreal,nele,& is,idate,i_dw,i_rw,i_sst,i_tcp,i_gps,i_uv,i_ps,i_lag,i_light,& i_t,i_pw,i_q,i_co,i_gust,i_vis,i_ref,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_cldtot,iobs,nprt,ii,jj integer(i_kind) it,ier,istatus real(r_quad):: zjo real(r_kind),dimension(40,ndat):: aivals1 real(r_kind),dimension(7,jpch_rad):: stats1 real(r_kind),dimension(9,jpch_oz):: stats_oz1 real(r_kind),dimension(9,jpch_co):: stats_co1 real(r_kind),dimension(npres_print,nconvtype,5,3):: bwork1 real(r_kind),allocatable,dimension(:,:):: awork1 real(r_kind),dimension(:,:,:),pointer:: ges_tv_it=>NULL() real(r_kind),dimension(:,:,:),pointer:: ges_q_it =>NULL() character(len=*),parameter:: myname='setuprhsall' logical,parameter:: OBSDIAGS_RELOAD = .false. !logical,parameter:: OBSDIAGS_RELOAD = .true. logical:: opened character(len=256):: tmpname,tmpaccess,tmpform if(.not.init_pass .and. .not.lobsdiag_allocated) call die('setuprhsall','multiple lobsdiag_allocated',lobsdiag_allocated) !****************************************************************************** ! Initialize timer call timer_ini('setuprhsall') ! Initialize variables and constants. first = jiter == jiterstart ! .true. on first outer iter last = jiter == miter+1 ! .true. following last outer iter llouter=.true. ! Set diagnostic output flag rad_diagsave = write_diag(jiter) .and. diag_rad pcp_diagsave = write_diag(jiter) .and. diag_pcp conv_diagsave = write_diag(jiter) .and. diag_conv ozone_diagsave= write_diag(jiter) .and. diag_ozone .and. ihave_oz co_diagsave = write_diag(jiter) .and. diag_co .and. ihave_co aero_diagsave = write_diag(jiter) .and. diag_aero light_diagsave= write_diag(jiter) .and. diag_light radardbz_diagsave = write_diag(jiter) .and. diag_radardbz i_ps = 1 i_uv = 2 i_t = 3 i_q = 4 i_pw = 5 i_rw = 6 i_dw = 7 i_gps= 8 i_sst= 9 i_tcp= 10 i_lag= 11 i_co = 12 i_gust=13 i_vis =14 i_pblh=15 i_wspd10m=16 i_td2m=17 i_mxtm=18 i_mitm=19 i_pmsl=20 i_howv=21 i_tcamt=22 i_lcbas=23 i_cldch=24 i_uwnd10m=26 i_vwnd10m=27 i_swcp=28 i_lwcp=29 i_light=30 i_dbz=31 i_cldtot=32 i_ref =i_cldtot allocate(awork1(7*nsig+100,i_ref)) if(.not.rhs_allocated) call rhs_alloc(aworkdim2=size(awork1,2)) ! Reset observation pointers if(init_pass) call destroyobs if(init_pass) call obsdiags_reset(obsdiags_keep=lobsdiagsave) ! replacing destroyobs() ! Read observation diagnostics if available if (l4dvar) then getodiag=(.not.lobserver) .or. (lobserver.and.jiter>1) clfile='obsdiags.ZZZ' if (lobsensfc .and. .not.lsensrecompute) then write(clfile(10:12),'(I3.3)') miter !call read_obsdiags(clfile) call obsdiags_read(clfile,mPEs=mPEs_observer) ! replacing read_obsdiags() call inquire_obsdiags(miter) else if (getodiag) then if (.not.lobserver) then write(clfile(10:12),'(I3.3)') jiter !call read_obsdiags(clfile) call obsdiags_read(clfile,mPEs=mPEs_observer) ! replacing read_obsdiags() call inquire_obsdiags(miter) endif endif endif if(init_pass) then if (allocated(obscounts)) then write(6,*)'setuprhsall: obscounts allocated' call stop2(285) end if allocate(obscounts(nobs_type,nobs_bins)) endif if (jiter>1.and.lobskeep) then nobskeep=1 else nobskeep=0 endif ! The 3d pressure and geopotential grids are initially loaded at ! the end of the call to read_guess. Thus, we don't need to call ! load_prsges and load_geop_hgt on the first outer loop. We need ! to update these 3d pressure arrays on all subsequent outer loops. ! Hence, the conditional call to load_prsges and load_geop_hgt if (lobserver .or. jiter>jiterstart) then ! Get sensible temperature (after bias correction's been applied) do it=1,nfldsig ier=0 call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'tv',ges_tv_it,istatus);ier=ier+istatus call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'q' ,ges_q_it ,istatus);ier=ier+istatus if(ier/=0) exit ges_tsen(:,:,:,it)= ges_tv_it(:,:,:)/(one+fv*max(zero,ges_q_it(:,:,:))) enddo ! Load 3d subdomain pressure arrays from the guess fields call load_prsges ! Compute 3d subdomain geopotential heights from the guess fields call load_geop_hgt ! if (sfcmod_gfs .or. sfcmod_mm5) then ! if (mype==0) write(6,*)'COMPUTE_DERIVED: call load_fact10' ! call load_fact10 ! endif endif ! Compute 2d subdomain pbl heights from the guess fields if (wrf_mass_regional) then call load_gsdpbl_hgt(mype) else if (nems_nmmb_regional .or. fv3_regional) then if (l_PBL_pseudo_SurfobsT .or. l_PBL_pseudo_SurfobsQ .or. l_PBL_pseudo_SurfobsUV) then call load_gsdpbl_hgt(mype) end if endif ! Compute derived quantities on grid if(.not.cmaq_regional) call compute_derived(mype,init_pass) ! ------------------------------------------------------------------------ if ( (l4dvar.and.lobserver) .or. .not.l4dvar ) then ! Init for Lagrangian data assimilation (gather winds and NL integration) call lag_presetup() ! Save state for inner loop if in 4Dvar observer mode if (l4dvar.and.lobserver) then call lag_state_write() end if ! Reset observation pointers. This is assumed by setup*() routines. do ii=1,size(obsdiags,2) do jj=1,size(obsdiags,1) obsdiags(jj,ii)%tail => NULL() enddo enddo lunin=1 open(lunin,file=obs_setup,form='unformatted') rewind lunin ! If requested, create conventional diagnostic files if(conv_diagsave.and.binary_diag)then write(string,900) jiter 900 format('conv_',i2.2) diag_conv_file=trim(dirname) // trim(string) if(init_pass) then open(7,file=trim(diag_conv_file),form='unformatted',status='unknown',position='rewind') else ! open(7,file=trim(diag_conv_file),form='unformatted',status='old',position='append') ! Without a close(7) until the last_pass=.true., the same file ! is expected to remain open "asis", equivalent to an "append" ! position through a re-open(). Therefore, a sequence of ! verification steps are taken to replace the earlier open() ! statement, to avoid re-open() without a close(). inquire(unit=7,opened=opened) if(opened) then inquire(unit=7,name=tmpname,form=tmpform,access=tmpaccess) tmpname=basename(tmpname) if(trim(tmpname)/=trim(diag_conv_file)) then call perr(myname,'unexpectly occupied, unit =',7) call perr(myname,' diag_conv_file =',trim(diag_conv_file)) call perr(myname,' inquire(unit=7, name= )',trim(tmpname)) call perr(myname,' inquire(unit=7, form= )',trim(tmpform)) call perr(myname,' inquire(unit=7,access= )',trim(tmpaccess)) call die(myname) endif else call perr(myname,'unexpectly closed, unit =',7) call perr(myname,' diag_conv_file =',trim(diag_conv_file)) call die(myname) endif endif idate=iadate(4)+iadate(3)*100+iadate(2)*10000+iadate(1)*1000000 if(init_pass .and. mype == 0)write(7)idate end if if (newpc4pred) then ostats=zero rstats=zero_quad end if if (aircraft_t_bc_pof .or. aircraft_t_bc) then ostats_t=zero_quad rstats_t=zero_quad end if ! If requested, create lightning diagnostic files if(light_diagsave)then write(string,500) jiter 500 format('light_',i2.2) diag_light_file=trim(dirname) // trim(string) if(init_pass) then open(55,file=trim(diag_light_file),form='unformatted',status='unknown',position='rewind') else open(55,file=trim(diag_light_file),form='unformatted',status='old',position='append') endif idate=iadate(4)+iadate(3)*100+iadate(2)*10000+iadate(1)*1000000 if(init_pass .and. mype == 0)write(55)idate end if ! Loop over data types to process do is=1,ndat nobs=nsat1(is) if(nobs > 0)then read(lunin,iostat=ier) obstype,isis,nreal,nchanl ! if(mype == mype_diaghdr(is)) then ! write(6,300) obstype,isis,nreal,nchanl !300 format(' SETUPALL:,obstype,isis,nreal,nchanl=',a12,a20,i5,i5) ! endif if(ier/=0) call die('setuprhsall','read(), iostat =',ier) nele=nreal+nchanl ! Set up for radiance data if(ditype(is) == 'rad')then call setuprad(lunin,& mype,aivals,stats,nchanl,nreal,nobs,& obstype,isis,is,rad_diagsave,init_pass,last_pass) ! Set up for aerosol data else if(ditype(is) == 'aero')then call setupaod(lunin,& mype,nchanl,nreal,nobs,& obstype,isis,is,aero_diagsave,init_pass) ! Set up for precipitation data else if(ditype(is) == 'pcp')then call setuppcp(lunin,mype,& aivals,nele,nobs,obstype,isis,is,pcp_diagsave,init_pass) ! Set up conventional data else if(ditype(is) == 'conv')then ! Set up temperature data if(obstype=='t')then call setupt(lunin,mype,bwork,awork(1,i_t),nele,nobs,is,conv_diagsave) ! Set up uv wind data else if(obstype=='uv')then call setupw(lunin,mype,bwork,awork(1,i_uv),nele,nobs,is,conv_diagsave) ! Set up wind speed data else if(obstype=='spd')then call setupspd(lunin,mype,bwork,awork(1,i_uv),nele,nobs,is,conv_diagsave) ! Set up surface pressure data else if(obstype=='ps')then call setupps(lunin,mype,bwork,awork(1,i_ps),nele,nobs,is,conv_diagsave) ! Set up tc-mslp data else if(obstype=='tcp')then call setuptcp(lunin,mype,bwork,awork(1,i_tcp),nele,nobs,is,conv_diagsave) ! Set up moisture data else if(obstype=='q') then call setupq(lunin,mype,bwork,awork(1,i_q),nele,nobs,is,conv_diagsave) ! Set up lidar wind data else if(obstype=='dw')then call setupdw(lunin,mype,bwork,awork(1,i_dw),nele,nobs,is,conv_diagsave) ! Set up radar wind data else if(obstype=='rw')then call setuprw(lunin,mype,bwork,awork(1,i_rw),nele,nobs,is,conv_diagsave) ! Set up radar reflectivity data else if(obstype=='dbz')then call setupdbz(lunin,mype,bwork,awork(1,i_dbz),nele,nobs,is,radardbz_diagsave,init_pass) ! Set up total precipitable water (total column water) data else if(obstype=='pw')then call setuppw(lunin,mype,bwork,awork(1,i_pw),nele,nobs,is,conv_diagsave) ! Set up conventional sst data else if(obstype=='sst' .and. getindex(svars2d,'sst')>0) then call setupsst(lunin,mype,bwork,awork(1,i_sst),nele,nobs,is,conv_diagsave) ! Set up conventional lagrangian data else if(obstype=='lag') then call setuplag(lunin,mype,bwork,awork(1,i_lag),nele,nobs,is,conv_diagsave) else if(obstype == 'pm2_5')then call setuppm2_5(lunin,mype,nele,nobs,isis,is,conv_diagsave) else if(obstype == 'pm10')then call setuppm10(lunin,mype,nele,nobs,isis,is,conv_diagsave) ! Set up conventional wind gust data else if(obstype=='gust' .and. getindex(svars2d,'gust')>0) then call setupgust(lunin,mype,bwork,awork(1,i_gust),nele,nobs,is,conv_diagsave) ! Set up conventional visibility data else if(obstype=='vis' .and. getindex(svars2d,'vis')>0) then call setupvis(lunin,mype,bwork,awork(1,i_vis),nele,nobs,is,conv_diagsave) ! Set up conventional pbl height data else if(obstype=='pblh' .and. getindex(svars2d,'pblh')>0) then call setuppblh(lunin,mype,bwork,awork(1,i_pblh),nele,nobs,is,conv_diagsave) ! Set up conventional wspd10m data else if(obstype=='wspd10m' .and. getindex(svars2d,'wspd10m')>0) then call setupwspd10m(lunin,mype,bwork,awork(1,i_wspd10m),nele,nobs,is,conv_diagsave) ! Set up conventional td2m data else if(obstype=='td2m' .and. getindex(svars2d,'td2m')>0) then call setuptd2m(lunin,mype,bwork,awork(1,i_td2m),nele,nobs,is,conv_diagsave) ! Set up conventional mxtm data else if(obstype=='mxtm' .and. getindex(svars2d,'mxtm')>0) then call setupmxtm(lunin,mype,bwork,awork(1,i_mxtm),nele,nobs,is,conv_diagsave) ! Set up conventional mitm data else if(obstype=='mitm' .and. getindex(svars2d,'mitm')>0) then call setupmitm(lunin,mype,bwork,awork(1,i_mitm),nele,nobs,is,conv_diagsave) ! Set up conventional pmsl data else if(obstype=='pmsl' .and. getindex(svars2d,'pmsl')>0) then call setuppmsl(lunin,mype,bwork,awork(1,i_pmsl),nele,nobs,is,conv_diagsave) ! Set up conventional howv data else if(obstype=='howv' .and. getindex(svars2d,'howv')>0) then call setuphowv(lunin,mype,bwork,awork(1,i_howv),nele,nobs,is,conv_diagsave) ! Set up total cloud amount data else if(obstype=='tcamt' .and. getindex(svars2d,'tcamt')>0) then call setuptcamt(lunin,mype,bwork,awork(1,i_tcamt),nele,nobs,is,conv_diagsave) ! Set up base height of lowest cloud seen else if(obstype=='lcbas' .and. getindex(svars2d,'lcbas')>0) then call setuplcbas(lunin,mype,bwork,awork(1,i_lcbas),nele,nobs,is,conv_diagsave) ! Set up conventional cldch data else if(obstype=='cldch' .and. getindex(svars2d,'cldch')>0) then call setupcldch(lunin,mype,bwork,awork(1,i_cldch),nele,nobs,is,conv_diagsave) ! Set up conventional uwnd10m data else if(obstype=='uwnd10m' .and. getindex(svars2d,'uwnd10m')>0) then call setupuwnd10m(lunin,mype,bwork,awork(1,i_uwnd10m),nele,nobs,is,conv_diagsave) ! Set up conventional vwnd10m data else if(obstype=='vwnd10m' .and. getindex(svars2d,'vwnd10m')>0) then call setupvwnd10m(lunin,mype,bwork,awork(1,i_vwnd10m),nele,nobs,is,conv_diagsave) ! Set up metar cloud pseudo obs else if(obstype=='mta_cld') then if(i_cloud_q_innovation==2) then call setupcldtot(lunin,mype,bwork,awork(1,i_cldtot),nele,nobs,is,conv_diagsave) else read(lunin,iostat=ier) if(ier/=0) call die('setuprhsall','read(), iostat =',ier) endif ! skip this kind of data because they are not used in the var analysis else if(obstype == 'gos_ctp' .or. & obstype == 'rad_ref' .or. obstype=='lghtn' .or. & obstype == 'larccld' .or. obstype == 'larcglb') then read(lunin,iostat=ier) if(ier/=0) call die('setuprhsall','read(), iostat =',ier) ! else write(6,*) 'Warning, unknown data type in setuprhsall,', obstype end if ! Set up solid/liquid-water content path data else if(ditype(is) == 'wcp')then if(obstype=='swcp')then call setupswcp(lunin,mype,bwork,awork(1,i_swcp),nele,nobs,is,conv_diagsave) else if (obstype=='lwcp')then call setuplwcp(lunin,mype,bwork,awork(1,i_lwcp),nele,nobs,is,conv_diagsave) endif ! set up ozone (sbuv/omi/mls) data else if(ditype(is) == 'ozone' .and. ihave_oz)then if (obstype == 'o3lev' .or. index(obstype,'mls')/=0 ) then call setupozlev(lunin,mype,stats_oz,nchanl,nreal,nobs,& obstype,isis,is,ozone_diagsave,init_pass) else call setupozlay(lunin,mype,stats_oz,nchanl,nreal,nobs,& obstype,isis,is,ozone_diagsave,init_pass) end if ! Set up co (mopitt) data else if(ditype(is) == 'co')then call setupco(lunin,mype,stats_co,nchanl,nreal,nobs,& obstype,isis,is,co_diagsave,init_pass) ! Set up GPS local refractivity data else if(ditype(is) == 'gps')then if(obstype=='gps_ref')then call setupref(lunin,mype,awork(1,i_gps),nele,nobs,toss_gps_sub,is,init_pass,last_pass,conv_diagsave) ! Set up GPS local bending angle data else if(obstype=='gps_bnd')then call setupbend(lunin,mype,awork(1,i_gps),nele,nobs,toss_gps_sub,is,init_pass,last_pass,conv_diagsave) end if ! Set up lightning (GOES/GLM) data else if(ditype(is) == 'light')then if(obstype == 'goes_glm')then call setuplight(lunin,mype,bwork,awork,nele,nobs,is,light_diagsave) end if endif ! ditype end if !if(nobs > 0)then end do !loop over all data types close(lunin) ! run cloud analysis in observer if(i_gsdcldanal_type==7) then call gsdcloudanalysis(mype) ! Write output analysis files call write_all(-1,mype) call prt_guess('analysis') endif else ! Init for Lagrangian data assimilation (read saved parameters) call lag_state_read() endif ! < lobserver > lobsdiag_allocated=.true. if(.not.last_pass) then call timer_fnl('setuprhsall') return endif ! Deallocate wind field array for Lagrangian data assimilation call lag_destroy_uv() ! Finalize qc and accumulate statistics for GPSRO data call gpsStats_genstats(bwork,awork(:,i_gps),toss_gps_sub,conv_diagsave,mype) call gpsStats_destroy() ! replacing ... ! -- call genstats_gps(bwork,awork(1,i_gps),toss_gps_sub,conv_diagsave,mype) if (conv_diagsave.and.binary_diag) close(7) if (light_diagsave) close(55) if(.not.(l_PBL_pseudo_SurfobsT .or. l_PBL_pseudo_SurfobsQ .or. & l_PBL_pseudo_SurfobsUV .or. (i_cloud_q_innovation==2)) ) then call obsdiags_sort() endif ! for temporary testing purposes, _write and _read. if(OBSDIAGS_RELOAD) then call obsdiags_write('obsdiags.ttt',force=.true.) ! call Barrier() before obsdiags_read(), to make sure all PEs have ! finished their obsdiags_write(). if(mPEs_observer>0) call MPI_Barrier(mpi_comm_world,ier) call obsdiags_read('obsdiags.ttt',mPEs=mPEs_observer,force=.true., & ignore_iter=.true.) call inquire_obsdiags(miter) endif ! call inquire_obsdiags(miter) ! Collect information for preconditioning if (newpc4pred) then call mpl_allreduce(jpch_rad,rpvals=ostats) call mpl_allreduce(npred,jpch_rad,rstats) end if ! Collect information for aircraft data if (aircraft_t_bc_pof .or. aircraft_t_bc) then ! call mpl_allreduce(npredt,max_tail,ostats_t) ! call mpl_allreduce(npredt,max_tail,rstats_t) call mpl_allreduce(npredt,ntail,ostats_t) call mpl_allreduce(npredt,ntail,rstats_t) end if if (newpc4pred .or. aircraft_t_bc_pof .or. aircraft_t_bc) then call reset_predictors_var end if ! Collect satellite and precip. statistics call mpi_reduce(aivals,aivals1,size(aivals1),mpi_rtype,mpi_sum,mype_rad, & mpi_comm_world,ierror) call mpi_reduce(stats,stats1,size(stats1),mpi_rtype,mpi_sum,mype_rad, & mpi_comm_world,ierror) if (ihave_oz) call mpi_reduce(stats_oz,stats_oz1,size(stats_oz1),mpi_rtype,mpi_sum,mype_oz, & mpi_comm_world,ierror) if (ihave_co) call mpi_reduce(stats_co,stats_co1,size(stats_co1),mpi_rtype,mpi_sum,mype_co, & mpi_comm_world,ierror) ! Collect conventional data statistics call mpi_allreduce(bwork,bwork1,size(bwork1),mpi_rtype,mpi_sum,& mpi_comm_world,ierror) call mpi_allreduce(awork,awork1,size(awork1),mpi_rtype,mpi_sum, & mpi_comm_world,ierror) ! Compute and print statistics for radiance, precipitation, and ozone data. ! These data types are NOT processed when running in 2dvar mode. Hence ! the check on the 2dvar flag below. if ( (l4dvar.and.lobserver) .or. .not.l4dvar ) then if (.not.twodvar_regional) then ! Compute and print statistics for radiance data if(mype==mype_rad) call statsrad(aivals1,stats1,ndata) ! Compute and print statistics for precipitation data if(mype==mype_rad) call statspcp(aivals1,ndata) ! Compute and print statistics for ozone if (mype==mype_oz .and. ihave_oz) call statsoz(stats_oz1,ndata) ! Compute and print statistics for carbon monoxide !???? if (mype==mype_co .and. ihave_co) call statsco(stats_co1,bwork1,awork1(1,i_co),ndata) endif ! Compute and print statistics for "conventional" data 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_ref,bwork1,awork1,ndata) ! Compute and print statistics for "lightning" data if (mype==mype_light) call statslight(mype,i_light,bwork,awork,i_ref,ndata) endif ! < .not. lobserver > deallocate(awork1) call rhs_dealloc() ! destroy the workspace: awork, bwork, etc. ! Print Jo table nprt=2 llouter=.true. if(luse_obsdiag)call evaljo(zjo,iobs,nprt,llouter) ! If only performing sst retrieval, end program execution if(retrieval)then deallocate(fbias) if(mype==0)then write(6,*)'SETUPRHSALL: normal completion for retrieval' call w3tage('GLOBAL_SSI') end if call mpi_finalize(ierror) stop end if ! Finalize timer call timer_fnl('setuprhsall') return end subroutine setuprhsall