PROGRAM arpsdas,175
!
!##################################################################
!##################################################################
!###### ######
!###### PROGRAM ARPSDAS ######
!###### ARPS Data Analysis System ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
! PURPOSE:
!
! This program does three-dimensional data analyses on
! the ARPS grid. It is loosely based on the rewritten
! OLAPS surface analysis.
!
! The driver establishes sizing of variables
! and adjustable parameters.
! It calls a routine to set up the grid, calls the data readers,
! a routine to analyse the data and another to write it out.
!
! AUTHOR:
!
! Keith Brewster, CAPS, March, 1994
! OLAPS sfc analysis Keith Brewster, CAPS, March, 1994
!
! MODIFICATION HISTORY:
! 01/17/96 (KB) Named changed to ADAS, dimensions imported from
! adas.dims and input variables from namelist
! rather than parameter statements in this code.
! 11/27/96 (KB) Cleaned-up to remove unused arrays.
! 01/25/96 (KB) Added tem variables passed to anxiter to accomodate
! new version of anxiter to speed-up execution.
! 06/15/97 (KB) Update for version 4.3.0, new arguments to call
! to INITVAR.
! 07/16/97 (KB) Added nt to calling list of INITVAR and made nt
! equal to one to save memory.
! 05/26/98 (KB) Modified logic in superadiabatic check. (kb)
! 08/31/98 (KB) Added code for nudging to official version (kb)
! 09/22/98 Jason Levit
! Updated call in initgrdvar for arps4.4.0, including
! change of dimension for qcumsrc.
! 12/07/98 Added background-based surface correlations to
! account for cumulus-parameterization small-scale changes.
! 12/20/98 Changed moisture analysis variable to qv.
! Error stats are still in RH, converted to qv based on
! background value.
! 12/30/98 Added processing of qr from retrieval data.
!
! 2000/04/14 (Gene Bassett)
! Merged the ADAS input file back into the ARPS input file.
!
! 2000/05/09 (Gene Bassett)
! Converted to F90, creating allocation and adas main
! subroutines.
!
! 2000/07/28 (Ming Xue)
! Converted to F90 free format. Did away with adaasamin for
! better compilation efficiency. Use ALLOCATABLE instead of
! POINTER allocation to avoid double memory usage.
!
! 07/23/2001 (K. Brewster)
! Added dynamic array allocation for the IAU increment arrays.
! Made corresponding changes to the increment handling
! routines to account for the array allocation.
!
! 08/07/2001 (K. Brewster)
! Added a preliminary calculation of vertical velocity
! before the cloud analyis (before cmpclddrv call).
!
! 11/01/2001 (K. Brewster)
! Renamed radar_ref_3d array to ref_mos_3d since its a mosaic.
! Removed potentially very large 4-d temporary arrays for
! radar by reorganizing logic for radar mosaic creation.
!
! 04/10/2002 (K. Brewster)
! Added code for soil surface temperature adjustment.
! Includes allocation and deallocation of arrays temporarily
! needed for radiation calculations.
!
! 05/16/2002 (K. Brewster)
! Added BMJ arrays to initgrdvar call to conform with
! latest updates to that subroutine. Added deallocate after
! the call to initgrdvar for some arrays not needed after
! that call.
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INCLUDE 'alloc.inc'
INCLUDE 'adas.inc'
INCLUDE 'nudging.inc'
INCLUDE 'exbc.inc'
INCLUDE 'radcst.inc'
REAL :: exbcbuf( 1 ) ! EXBC buffer array (unused)
!
!-----------------------------------------------------------------------
!
! Arrays defining model grid
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: x (:) ! The x-coord. of the physical and
! computational grid. Defined at u-point.
REAL, ALLOCATABLE :: y (:) ! The y-coord. of the physical and
! computational grid. Defined at v-point.
REAL, ALLOCATABLE :: z (:) ! The z-coord. of the computational grid.
! Defined at w-point on the staggered grid.
REAL, ALLOCATABLE :: zp (:,:,:) ! The physical height coordinate defined at
! w-point of the staggered grid.
REAL, ALLOCATABLE :: hterain(:,:) ! The height of the terrain.
REAL, ALLOCATABLE :: mapfct(:,:,:) ! Map factors at scalar, u and v points
REAL, ALLOCATABLE :: j1 (:,:,:) ! Coordinate transformation Jacobian defined
! as - d( zp )/d( x ).
REAL, ALLOCATABLE :: j2 (:,:,:) ! Coordinate transformation Jacobian defined
! as - d( zp )/d( y ).
REAL, ALLOCATABLE :: j3 (:,:,:) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ).
REAL, ALLOCATABLE :: aj3x (:,:,:) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE X-DIR.
REAL, ALLOCATABLE :: aj3y (:,:,:) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE Y-DIR.
REAL, ALLOCATABLE :: aj3z (:,:,:) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE Z-DIR.
REAL, ALLOCATABLE :: j3inv (:,:,:) ! Inverse of j3
!
!-----------------------------------------------------------------------
!
! ARPS Time-dependent variables
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: u (:,:,:) ! Total u-velocity (m/s)
REAL, ALLOCATABLE :: v (:,:,:) ! Total v-velocity (m/s)
REAL, ALLOCATABLE :: w (:,:,:) ! Total w-velocity (m/s)
REAL, ALLOCATABLE :: wcont (:,:,:) ! Contravariant vertical velocity (m/s)
REAL, ALLOCATABLE :: ptprt (:,:,:) ! Perturbation potential temperature (K)
REAL, ALLOCATABLE :: pprt (:,:,:) ! Perturbation pressure (Pascal)
REAL, ALLOCATABLE :: qv (:,:,:) ! Water vapor specific humidity (kg/kg)
REAL, ALLOCATABLE :: qc (:,:,:) ! Cloud water mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qr (:,:,:) ! Rain water mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qi (:,:,:) ! Cloud ice mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qs (:,:,:) ! Snow mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qh (:,:,:) ! Hail mixing ratio (kg/kg)
REAL, ALLOCATABLE :: tke (:,:,:) ! Turbulent Kinetic Energy ((m/s)**2)
REAL, ALLOCATABLE :: ubar (:,:,:) ! Base state u-velocity (m/s)
REAL, ALLOCATABLE :: vbar (:,:,:) ! Base state v-velocity (m/s)
REAL, ALLOCATABLE :: ptbar (:,:,:) ! Base state potential temperature (K)
REAL, ALLOCATABLE :: pbar (:,:,:) ! Base state pressure (Pascal).
REAL, ALLOCATABLE :: rhostr(:,:,:) ! Base state density rhobar times j3.
REAL, ALLOCATABLE :: qvbar (:,:,:) ! Base state water vapor specific
! humidity(kg/kg)
REAL, ALLOCATABLE :: udteb (:,:) ! T-tendency of u at e-boundary (m/s**2)
REAL, ALLOCATABLE :: udtwb (:,:) ! T-tendency of u at w-boundary (m/s**2)
REAL, ALLOCATABLE :: vdtnb (:,:) ! T-tendency of v at n-boundary (m/s**2)
REAL, ALLOCATABLE :: vdtsb (:,:) ! T-tendency of v at s-boundary (m/s**2)
REAL, ALLOCATABLE :: pdteb (:,:) ! T-tendency of pprt at e-boundary (PASCAL/s)
REAL, ALLOCATABLE :: pdtwb (:,:) ! T-tendency of pprt at w-boundary (PASCAL/s)
REAL, ALLOCATABLE :: pdtnb (:,:) ! T-tendency of pprt at n-boundary (PASCAL/s)
REAL, ALLOCATABLE :: pdtsb (:,:) ! T-tendency of pprt at s-boundary (PASCAL/s)
REAL, ALLOCATABLE :: trigs1(:) ! Array containing pre-computed trig
! function for use in fft.
REAL, ALLOCATABLE :: trigs2(:) ! Array containing pre-computed trig
! function for use in fft.
INTEGER, ALLOCATABLE :: ifax1(:) ! Array containing the factors of nx.
INTEGER, ALLOCATABLE :: ifax2(:) ! Array containing the factors of ny.
REAL, ALLOCATABLE :: vwork1 (:,:) ! 2-D work array for fftopt=2.
REAL, ALLOCATABLE :: vwork2 (:,:) ! 2-D work array for fftopt=2.
REAL, ALLOCATABLE :: wsave1 (:) ! Work array for fftopt=2.
REAL, ALLOCATABLE :: wsave2 (:) ! Work array for fftopt=2.
REAL, ALLOCATABLE :: ppi(:,:,:) ! Exner function.
REAL, ALLOCATABLE :: csndsq(:,:,:) ! Speed of sound squared
REAL, ALLOCATABLE :: radfrc(:,:,:) ! Radiation forcing (K/s)
REAL, ALLOCATABLE :: radsw(:,:) ! Solar radiation reaching the surface
REAL, ALLOCATABLE :: rnflx(:,:) ! Net absorbed radiation by the surface
!
!-----------------------------------------------------------------------
!
! ARPS Surface variables:
!
!-----------------------------------------------------------------------
!
INTEGER, ALLOCATABLE :: soiltyp (:,:,:) ! Soil type
REAL, ALLOCATABLE :: stypfrct(:,:,:) ! Soil type fraction
INTEGER, ALLOCATABLE :: vegtyp (:,:) ! Vegetation type
REAL, ALLOCATABLE :: lai (:,:) ! Leaf Area Index
REAL, ALLOCATABLE :: roufns (:,:) ! Surface roughness
REAL, ALLOCATABLE :: veg (:,:) ! Vegetation fraction
REAL, ALLOCATABLE :: tsfc(:,:,:) ! Temperature at ground (K)
! (in top 1 cm layer)
REAL, ALLOCATABLE :: qvsfc(:,:,:) ! Effective qv at sfc.
REAL, ALLOCATABLE :: tsoil(:,:,:) ! Deep soil temperature (K)
! (in deep 1 m layer)
REAL, ALLOCATABLE :: wetsfc(:,:,:) ! Surface soil moisture in the top
! 1 cm layer
REAL, ALLOCATABLE :: wetdp(:,:,:) ! Deep soil moisture in the deep
! 1 m layer
REAL, ALLOCATABLE :: wetcanp(:,:,:) ! Canopy water amount
REAL, ALLOCATABLE :: snowdpth(:,:) ! Snow depth (:)
REAL, ALLOCATABLE :: ptcumsrc(:,:,:)! Source term in pt-equation due
! to cumulus parameterization
REAL, ALLOCATABLE :: qcumsrc(:,:,:,:) ! Source term in Qv equation due
! to cumulus parameterization
REAL, ALLOCATABLE :: w0avg(:,:,:) ! a closing running average vertical
! velocity in 10min for K-F scheme
REAL, ALLOCATABLE :: kfraincv(:,:) ! K-F convective rainfall (:)
INTEGER, ALLOCATABLE :: nca(:,:) ! K-F counter for CAPE release
REAL, ALLOCATABLE :: cldefi(:,:) ! BMJ cloud efficiency
REAL, ALLOCATABLE :: xland(:,:) ! BMJ land mask
! (1.0 = land, 2.0 = sea)
REAL, ALLOCATABLE :: bmjraincv(:,:) ! BMJ convective rainfall (cm)
REAL, ALLOCATABLE :: raing(:,:) ! Grid supersaturation rain
REAL, ALLOCATABLE :: rainc(:,:) ! Cumulus convective rain
REAL, ALLOCATABLE :: prcrate(:,:,:) ! precipitation rate (kg/(m**2*s))
! prcrate(:,:,:) = total precipitation rate
! prcrate(:,:,:) = grid scale precip. rate
! prcrate(:,:,:) = cumulus precip. rate
! prcrate(:,:,:) = microphysics precip. rate
REAL, ALLOCATABLE :: usflx (:,:) ! Surface flux of u-momentum (kg/(m*s**2))
REAL, ALLOCATABLE :: vsflx (:,:) ! Surface flux of v-momentum (kg/(m*s**2))
REAL, ALLOCATABLE :: ptsflx(:,:) ! Surface heat flux (K*kg/(m*s**2))
REAL, ALLOCATABLE :: qvsflx(:,:) ! Surface moisture flux (kg/(m**2*s))
!
!-----------------------------------------------------------------------
!
! Arrays for analysis increment updating (a type of nudging) output
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: uincr(:,:,:)
REAL, ALLOCATABLE :: vincr(:,:,:)
REAL, ALLOCATABLE :: wincr(:,:,:)
REAL, ALLOCATABLE :: pincr(:,:,:)
REAL, ALLOCATABLE :: ptincr(:,:,:)
REAL, ALLOCATABLE :: qvincr(:,:,:)
REAL, ALLOCATABLE :: qcincr(:,:,:)
REAL, ALLOCATABLE :: qrincr(:,:,:)
REAL, ALLOCATABLE :: qiincr(:,:,:)
REAL, ALLOCATABLE :: qsincr(:,:,:)
REAL, ALLOCATABLE :: qhincr(:,:,:)
!
!-----------------------------------------------------------------------
!
! Analysis scratch space
!
!-----------------------------------------------------------------------
!
INTEGER, ALLOCATABLE :: knt(:,:)
REAL, ALLOCATABLE :: rngsqi(:)
REAL, ALLOCATABLE :: wgtsum(:,:)
REAL, ALLOCATABLE :: zsum(:,:)
REAL, ALLOCATABLE :: sqsum(:,:)
!
!-----------------------------------------------------------------------
!
! Additional grid variables
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: xs(:)
REAL, ALLOCATABLE :: ys(:)
REAL, ALLOCATABLE :: zs(:,:,:)
REAL, ALLOCATABLE :: latgr(:,:)
REAL, ALLOCATABLE :: longr(:,:)
!
!-----------------------------------------------------------------------
!
! Temporary arrays
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: tem1 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem2 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem3 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem4 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem5 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem6 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem7 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem8 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem9 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem10 (:,:,:) ! Temporary work array.
INTEGER, ALLOCATABLE :: ibegin(:)
INTEGER, ALLOCATABLE :: iend(:)
!
!-----------------------------------------------------------------------
!
! Radiation variables only needed for adjtsfc
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: rsirbm(:,:)
REAL, ALLOCATABLE :: rsirdf(:,:)
REAL, ALLOCATABLE :: rsuvbm(:,:)
REAL, ALLOCATABLE :: rsuvdf(:,:)
REAL, ALLOCATABLE :: cosz(:,:)
REAL, ALLOCATABLE :: cosss(:,:)
REAL, ALLOCATABLE :: fdirir(:,:)
REAL, ALLOCATABLE :: fdifir(:,:)
REAL, ALLOCATABLE :: fdirpar(:,:)
REAL, ALLOCATABLE :: fdifpar(:,:)
REAL, ALLOCATABLE :: st4(:,:)
REAL, ALLOCATABLE :: radbuf(:)
REAL, ALLOCATABLE :: tem11 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem12 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem13 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem14 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem15 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem16 (:,:,:) ! Temporary work array.
!
!-----------------------------------------------------------------------
!
! Cloud analysis variables other than model's ones:
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: ref_mos_3d (:,:,:)
! 3D remapped radar reflectivity field as input
! may be modified by sate. VIS as output
REAL, ALLOCATABLE :: w_cld (:,:,:) ! vertical velocity (m/s) in clouds
INTEGER(kind=selected_int_kind(4)), ALLOCATABLE :: cldpcp_type_3d(:,:,:)
! cloud/precip type field
INTEGER, ALLOCATABLE :: icing_index_3d(:,:,:) ! icing severity index
LOGICAL, ALLOCATABLE :: l_mask_pcptype(:,:) ! analyze precip type
!
!-----------------------------------------------------------------------
!
! Analysis variables
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: anx(:,:,:,:)
!
!-----------------------------------------------------------------------
!
! Cross-category correlations
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: xcor(:,:)
INTEGER, ALLOCATABLE :: icatg(:,:)
!
!-----------------------------------------------------------------------
!
! ARPS dimensions:
!
! nx, ny, nz: Dimensions of analysis grid.
!
!-----------------------------------------------------------------------
!
INTEGER :: nx ! Number of grid points in the x-direction
INTEGER :: ny ! Number of grid points in the y-direction
INTEGER :: nz ! Number of grid points in the z-direction
INTEGER :: nxndg ! Number of x grid points for IAU
INTEGER :: nyndg ! Number of y grid points for IAU
INTEGER :: nzndg ! Number of z grid points for IAU
!
!-----------------------------------------------------------------------
!
! Soil types.
!
!-----------------------------------------------------------------------
!
INTEGER :: nstyps ! Number of soil types
!
!-----------------------------------------------------------------------
!
! Misc
!
!-----------------------------------------------------------------------
!
INTEGER :: nt ! Number of time levels of data
INTEGER :: ncat ! Number of correlation categories
!-----------------------------------------------------------------------
!
! Analysis parameters
!
!-----------------------------------------------------------------------
!
REAL :: rhmin
INTEGER :: iqspr,klim
PARAMETER (rhmin=0.05, & ! rh safety net value to prevent neg qv
iqspr=3, & ! Use qobs of pstn to combine x,y,elev
klim= 1) ! Min of one other sfc station for QC
!
!-----------------------------------------------------------------------
!
! ARPS cloud cover analysis variables:
!
!-----------------------------------------------------------------------
!
! Remapping parameters of the remapped radar data grid
!
INTEGER :: strhopt_rad ! streching option
INTEGER :: mapproj_rad ! map projection indicator
REAL :: dx_rad,dy_rad,dz_rad,dzmin_rad ! grid spcngs
REAL :: ctrlat_rad,ctrlon_rad ! central lat and lon
REAL :: tlat1_rad,tlat2_rad,tlon_rad ! true lat and lon
REAL :: scale_rad ! map scale factor
!
!-----------------------------------------------------------------------
!
! Common block that stores remapping parameters for the radar
! data file.
!
!-----------------------------------------------------------------------
!
COMMON /remapfactrs_rad/ strhopt_rad,mapproj_rad
COMMON /remapfactrs_rad2/ dx_rad,dy_rad,dz_rad,dzmin_rad, &
ctrlat_rad,ctrlon_rad, &
tlat1_rad,tlat2_rad,tlon_rad,scale_rad
!
!-----------------------------------------------------------------------
!
! Non-dimensional smoothing coefficient
!
!-----------------------------------------------------------------------
!
INTEGER, PARAMETER :: tsfcadj=1
!
!-----------------------------------------------------------------------
!
! ARPS include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
INCLUDE 'bndry.inc'
INCLUDE 'grid.inc' ! Grid parameters
INCLUDE 'adjust.inc'
!
!-----------------------------------------------------------------------
!
! Surface Station variables
!
!-----------------------------------------------------------------------
!
INTEGER :: isrcsng(mx_sng)
INTEGER :: icatsng(mx_sng)
REAL :: latsng(mx_sng,ntime),lonsng(mx_sng,ntime)
REAL :: hgtsng(mx_sng,ntime)
REAL :: xsng(mx_sng),ysng(mx_sng)
INTEGER :: timesng(mx_sng,ntime)
CHARACTER (LEN=5) :: stnsng(mx_sng,ntime)
!
!-----------------------------------------------------------------------
!
! Surface (single-level) read-in observation variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=8) :: wx(mx_sng,ntime)
CHARACTER (LEN=8) :: csrcsng(mx_sng,ntime)
CHARACTER (LEN=1) :: store_emv(mx_sng,5,ntime)
CHARACTER (LEN=4) :: store_amt(mx_sng,5,ntime)
INTEGER :: kloud(mx_sng,ntime),idp3(mx_sng,ntime)
REAL :: store_hgt(mx_sng,5,ntime)
REAL :: obrdsng(mx_sng,nvar_sng,ntime)
REAL :: obsng(nvar_anx,mx_sng)
REAL :: odifsng(nvar_anx,mx_sng)
REAL :: corsng(nvar_anx,mx_sng)
REAL :: oanxsng(nvar_anx,mx_sng)
REAL :: thesng(mx_sng)
INTEGER :: ival(mx_sng)
!
!-----------------------------------------------------------------------
!
! Upper Air Station variables
!
!-----------------------------------------------------------------------
!
INTEGER :: isrcua(mx_ua)
REAL :: xua(mx_ua),yua(mx_ua)
REAL :: elevua(mx_ua)
REAL :: hgtua(nz_ua,mx_ua)
INTEGER :: nlevsua(mx_ua)
CHARACTER (LEN=5) :: stnua(mx_ua)
!
!-----------------------------------------------------------------------
!
! Upper-air observation variables
!
!-----------------------------------------------------------------------
!
REAL :: obsua(nvar_anx,nz_ua,mx_ua)
REAL :: odifua(nvar_anx,nz_ua,mx_ua)
REAL :: corua(nvar_anx,nz_ua,mx_ua)
REAL :: oanxua(nvar_anx,nz_ua,mx_ua)
REAL :: theua(nz_ua,mx_ua)
!
!-----------------------------------------------------------------------
!
! Radar site variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=5) :: stnrad(mx_rad)
INTEGER :: isrcrad(0:mx_rad)
REAL :: latrad(mx_rad),lonrad(mx_rad)
REAL :: elvrad(mx_rad)
!
!-----------------------------------------------------------------------
!
! Radar observation variables
!
!-----------------------------------------------------------------------
!
INTEGER :: irad(mx_colrad)
INTEGER :: nlevrad(mx_colrad)
REAL :: latradc(mx_colrad),lonradc(mx_colrad)
REAL :: xradc(mx_colrad),yradc(mx_colrad)
REAL :: distrad(mx_colrad),uazmrad(mx_colrad),vazmrad(mx_colrad)
REAL :: hgtradc(nz_rdr,mx_colrad)
REAL :: theradc(nz_rdr,mx_colrad)
REAL :: obsrad(nvar_radin,nz_rdr,mx_colrad)
REAL :: odifrad(nvar_rad,nz_rdr,mx_colrad)
REAL :: corrad(nvar_rad,nz_rdr,mx_colrad)
REAL :: oanxrad(nvar_rad,nz_rdr,mx_colrad)
!
!-----------------------------------------------------------------------
!
! Retrieval radar variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=5) :: stnret(mx_ret)
INTEGER :: isrcret(0:mx_ret)
REAL :: latret(mx_ret),lonret(mx_ret)
REAL :: elvret(mx_ret)
!
!-----------------------------------------------------------------------
!
! Retrieval observation variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iret(mx_colret)
INTEGER :: nlevret(mx_colret)
REAL :: latretc(mx_colret),lonretc(mx_colret)
REAL :: xretc(mx_colret),yretc(mx_colret)
! real distret(mx_colret),uazmret(mx_colret),vazmret(mx_colret)
REAL :: hgtretc(nz_ret,mx_colret)
REAL :: theretc(nz_ret,mx_colret)
REAL :: obsret(nvar_anx,nz_ret,mx_colret)
REAL :: qrret(nz_ret,mx_colret)
REAL :: odifret(nvar_anx,nz_ret,mx_colret)
REAL :: corret(nvar_anx,nz_ret,mx_colret)
REAL :: oanxret(nvar_anx,nz_ret,mx_colret)
!
!-----------------------------------------------------------------------
!
! Namen de variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=6) :: var_sfc(nvar_sng)
DATA var_sfc/ 't ','td ','dd ','ff ', &
'urot ','vrot ','pstn ','pmsl ', &
'alt ','ceil ','low ','cover ', &
'solar ','vis'/
CHARACTER (LEN=6) :: var_ua(nvar_ua)
DATA var_ua / 'press ', &
'temp ','dewpt ','dd ','ff '/
CHARACTER (LEN=6) :: var_anx(nvar_anx)
DATA var_anx/ 'u ','v ', &
'press ','theta ','qv '/
!
!-----------------------------------------------------------------------
!
! Source-dependent parameters
! Qsrc is the expected square error it is used for
! setting the QC threshold (qclim) and for relative
! weighting of observations.
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=80) :: srcback
REAL :: hqback(nz_tab)
INTEGER :: nlvqback
REAL :: qback(nvar_anx,nz_tab)
CHARACTER (LEN=80) :: srcsng_full(nsrc_sng)
REAL :: qsrcsng(nvar_anx,nsrc_sng)
CHARACTER (LEN=80) :: srcua_full(nsrc_ua)
INTEGER :: nlvuatab(nsrc_ua)
REAL :: huaqsrc(nz_tab,nsrc_ua)
REAL :: qsrcua(nvar_anx,nz_tab,nsrc_ua)
CHARACTER (LEN=80) :: srcrad_full(nsrc_rad)
REAL :: qsrcrad(nvar_radin,nsrc_rad)
CHARACTER (LEN=80) :: srcret_full(nsrc_ret)
INTEGER :: nlvrttab(nsrc_ret)
REAL :: hrtqsrc(nz_tab,nsrc_ret)
REAL :: qsrcret(nvar_anx,nz_tab,nsrc_ret)
!
!-----------------------------------------------------------------------
!
! Quality Control Parameters and Variables
!
!-----------------------------------------------------------------------
!
REAL :: rmiss
PARAMETER (rmiss=-99.)
!
INTEGER :: qualrdsng(mx_sng,nvar_sng,ntime)
INTEGER :: qualsng(nvar_anx,mx_sng)
REAL :: qobsng(nvar_anx,mx_sng)
REAL :: qcmulsng(nsrc_sng)
REAL :: qcthrsng(nvar_anx,nsrc_sng)
!
INTEGER :: qualua(nvar_anx,nz_ua,mx_ua)
REAL :: qobsua(nvar_anx,nz_ua,mx_ua)
REAL :: qcmulua(nsrc_ua)
REAL :: qcthrua(nvar_anx,nsrc_ua)
!
INTEGER :: qualrad(nvar_rad,nz_rdr,mx_colrad)
REAL :: qobsrad(nvar_rad,nz_rdr,mx_colrad)
REAL :: qcmulrad(nsrc_rad)
REAL :: qcthrrad(nvar_radin,nsrc_rad)
!
INTEGER :: qualret(nvar_anx,nz_ret,mx_colret)
REAL :: qobsret(nvar_anx,nz_ret,mx_colret)
REAL :: qcmulret(nsrc_ret)
REAL :: qcthrret(nvar_anx,nsrc_ret)
!
!-----------------------------------------------------------------------
!
! Climin is the minimum possible value of each observed variable.
! Climax is the maximum possible value of each observed variable.
! Used for surface data only.
!
!-----------------------------------------------------------------------
!
REAL :: climin(nvar_sng),climax(nvar_sng)
DATA climin / -50., -50., 0., 0., &
-100., -100., 700., 880., &
880., 0., 0., 0., &
0., 0./
DATA climax / 125., 125., 360., 100., &
100., 100., 1100., 1090., &
1090., 30000.,30000., 1., &
1500., 150./
!
!-----------------------------------------------------------------------
!
! Filenames and such
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=12) :: suffix
CHARACTER (LEN=128) :: froot
CHARACTER (LEN=132) :: qclist
CHARACTER (LEN=132) :: qcsave
CHARACTER (LEN=132) :: stats
INTEGER :: istatus,jstatus
INTEGER :: nobsng,nobsrd(ntime)
INTEGER :: i4timanx
INTEGER :: lenfnm,maxsuf,lensuf,dotloc,mxroot,lenroot
INTEGER :: iqclist,iqcsave,iwstat
CHARACTER (LEN=132) :: status_file
INTEGER n_used_sng(nsrc_sng), n_used_ua(nsrc_sng)
INTEGER jsta, klev, ngoodlev
!
!-----------------------------------------------------------------------
!
! Physical constants
!
!-----------------------------------------------------------------------
!
REAL :: kts2ms,mb2pa,f2crat
PARAMETER (kts2ms=0.514444, &
mb2pa=100., &
f2crat=(5./9.))
!
!-----------------------------------------------------------------------
!
! ptlapse is minimum potential temperature lapse rate to use in
! superadiabatic check. Used to ensure initial static stability.
!
!-----------------------------------------------------------------------
REAL :: ptlapse
PARAMETER (ptlapse = (1.0/4000.) ) ! deg K / m
!
!-----------------------------------------------------------------------
!
! Function f_qvsat and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
REAL :: f_qvsat
!fpp$ expand (f_qvsat)
!dir$ inline always f_qvsat
!
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=8) :: rdsource
CHARACTER (LEN=80) :: basdmpfn
CHARACTER (LEN=80) :: header
INTEGER :: i,j,k,ios,ktab,isrc,jsrc,ivar,ifile
INTEGER :: grdbas,lbasdmpf,rbufsz
INTEGER :: nobsua,ncolrad,ncolret
REAL :: tgrid,qvmin,qvsat,qsrcmax,rngsq,sprkm
REAL :: temp,tvbar,qvprt,qtot,pconst,pr1,pr2
!
INTEGER :: ixrad(mx_rad),jyrad(mx_rad)
INTEGER :: retqropt
PARAMETER (retqropt=1)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
nt = 1 ! Number of time levels of data
ncat = 4 ! Number of correlation categories
!-----------------------------------------------------------------------
! Set grid mode to non-nesting and grid index, mgrid, to 1.
!-----------------------------------------------------------------------
nestgrd=0
mgrid=1
!-----------------------------------------------------------------------
! read in ARPS namelists
!-----------------------------------------------------------------------
CALL initpara
(nx,ny,nz,nstyps)
IF (lbcopt == 2) THEN
WRITE (*,*) "INITADAS: resetting lbcopt to 1 ", &
"& lateral bc's to 4"
lbcopt = 1
ebc = 4
wbc = 4
nbc = 4
sbc = 4
END IF
!-----------------------------------------------------------------------
! Read in adas namelist parameters
!-----------------------------------------------------------------------
!
CALL initadas
IF(incrdmp > 0 ) THEN
nxndg=nx
nyndg=ny
nzndg=nz
ELSE
nxndg=1
nyndg=1
nzndg=1
END IF
!
!-----------------------------------------------------------------------
!
! Allocate adas arrays
!
!-----------------------------------------------------------------------
!
ALLOCATE( CALL check_alloc_status
(istatus, "adas:x")
ALLOCATE( CALL check_alloc_status
(istatus, "adas:y")
ALLOCATE( CALL check_alloc_status
(istatus, "adas:z")
ALLOCATE(hterain(nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:hterain")
ALLOCATE(mapfct (nx,ny,8),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:mapfct")
ALLOCATE(zp (nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:zp")
ALLOCATE(j1 (nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:j1")
ALLOCATE(j2 (nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:j2")
ALLOCATE(j3 (nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:j3")
ALLOCATE(aj3x(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:aj3x")
ALLOCATE(aj3y(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:aj3y")
ALLOCATE(aj3z(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:aj3z")
ALLOCATE(j3inv(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:j3inv")
ALLOCATE(u (nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:u")
ALLOCATE(v (nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:v")
ALLOCATE(w (nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:w")
ALLOCATE(wcont(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:wcont")
ALLOCATE(ptprt(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:ptprt")
ALLOCATE(pprt (nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:pprt")
ALLOCATE(qv (nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:qv")
ALLOCATE(qc (nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:qc")
ALLOCATE(qr (nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:qr")
ALLOCATE(qi (nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:qi")
ALLOCATE(qs (nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:qs")
ALLOCATE(qh (nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:qh")
ALLOCATE(tke (nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:tke")
ALLOCATE(ubar (nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:ubar")
ALLOCATE(vbar (nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:vbar")
ALLOCATE(ptbar(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:ptbar")
ALLOCATE(pbar (nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:pbar")
ALLOCATE(rhostr(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:rhostr")
ALLOCATE(qvbar(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:qvbar")
ALLOCATE(udteb(ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:udteb")
ALLOCATE(udtwb(ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:udtbw")
ALLOCATE(vdtnb(nx,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:vdtnb")
ALLOCATE(vdtsb(nx,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:vdtsb")
ALLOCATE(pdteb(ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:pdteb")
ALLOCATE(pdtwb(ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:pdtwb")
ALLOCATE(pdtnb(nx,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:pdtnb")
ALLOCATE(pdtsb(nx,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:pdtsb")
ALLOCATE(trigs1(3*(nx-1)/2+1),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:trigs1")
ALLOCATE(trigs2(3*(ny-1)/2+1),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:trigs2")
ALLOCATE(ifax1(13),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:ifax1")
ALLOCATE(ifax2(13),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:ifax2")
ALLOCATE(vwork1(nx+1,ny+1),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:vwork1")
ALLOCATE(vwork2(ny,nx+1),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:vwork2")
ALLOCATE(wsave1(3*(ny-1)+15),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:wsave1")
ALLOCATE(wsave2(3*(nx-1)+15),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:wsave2")
ALLOCATE(ppi (nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:ppi")
ALLOCATE(csndsq(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:csndsq")
ALLOCATE(radfrc(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:radfrc")
ALLOCATE(radsw (nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:radsw")
ALLOCATE(rnflx (nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:rnflx")
ALLOCATE(soiltyp(nx,ny,nstyps),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:soiltyp")
ALLOCATE(stypfrct(nx,ny,nstyps),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:stypfrct")
ALLOCATE(vegtyp (nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:vegtyp")
ALLOCATE(lai (nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:lai")
ALLOCATE(roufns (nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:roufns")
ALLOCATE(veg (nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:veg")
ALLOCATE(tsfc (nx,ny,0:nstyps),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:tsfc")
ALLOCATE(qvsfc (nx,ny,0:nstyps),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:qvsfc")
ALLOCATE(tsoil (nx,ny,0:nstyps),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:tsoil")
ALLOCATE(wetsfc (nx,ny,0:nstyps),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:wetsfc")
ALLOCATE(wetdp (nx,ny,0:nstyps),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:wetdp")
ALLOCATE(wetcanp(nx,ny,0:nstyps),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:wetcanp")
ALLOCATE(snowdpth(nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:snowdpth")
ALLOCATE(ptcumsrc(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:ptcumsrc")
ALLOCATE(qcumsrc (nx,ny,nz,5),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:qcumsrc")
ALLOCATE(w0avg (nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:w0avg")
ALLOCATE(nca (nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:nca")
ALLOCATE(kfraincv (nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:kfraincv")
ALLOCATE(cldefi (nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:cldefi")
ALLOCATE(xland (nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:xland")
ALLOCATE(bmjraincv (nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:bmjraincv")
ALLOCATE(raing (nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:raing")
ALLOCATE(rainc (nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:rainc")
ALLOCATE(prcrate(nx,ny,4),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:prcrate")
ALLOCATE(usflx (nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:usflx")
ALLOCATE(vsflx (nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:vsflx")
ALLOCATE(ptsflx(nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:ptsflx")
ALLOCATE(qvsflx(nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:qvsflx")
ALLOCATE(uincr(nxndg,nyndg,nzndg),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:uincr")
ALLOCATE(vincr(nxndg,nyndg,nzndg),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:vincr")
ALLOCATE(wincr(nxndg,nyndg,nzndg),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:wincr")
ALLOCATE(pincr(nxndg,nyndg,nzndg),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:pincr")
ALLOCATE(ptincr(nxndg,nyndg,nzndg),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:ptincr")
ALLOCATE(qvincr(nxndg,nyndg,nzndg),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:qvincr")
ALLOCATE(qcincr(nxndg,nyndg,nzndg),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:qcincr")
ALLOCATE(qrincr(nxndg,nyndg,nzndg),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:qrincr")
ALLOCATE(qiincr(nxndg,nyndg,nzndg),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:qiincr")
ALLOCATE(qsincr(nxndg,nyndg,nzndg),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:qsincr")
ALLOCATE(qhincr(nxndg,nyndg,nzndg),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:qhincr")
ALLOCATE(knt(nvar_anx,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:knt")
ALLOCATE(rngsqi(nvar_anx),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:rngsqi")
ALLOCATE(wgtsum(nvar_anx,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:wgtsum")
ALLOCATE(zsum(nvar_anx,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:zsum")
ALLOCATE(sqsum(nvar_anx,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:sqsum")
ALLOCATE(xs(nx),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:xs")
ALLOCATE(ys(ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:ys")
ALLOCATE(zs(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:zs")
ALLOCATE(latgr(nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:latgr")
ALLOCATE(longr(nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:longr")
ALLOCATE(tem1(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:tem1")
ALLOCATE(tem2(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:tem2")
ALLOCATE(tem3(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:tem3")
ALLOCATE(tem4(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:tem4")
ALLOCATE(tem5(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:tem5")
ALLOCATE(tem6(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:tem6")
ALLOCATE(tem7(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:tem7")
ALLOCATE(tem8(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:tem8")
ALLOCATE(tem9(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:tem9")
ALLOCATE(tem10(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:tem10")
ALLOCATE(ibegin(ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:ibegin")
ALLOCATE(iend(ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:iend")
ALLOCATE(anx(nx,ny,nz,nvar_anx),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:anx")
ALLOCATE(xcor(ncat,ncat),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:xcor")
ALLOCATE(icatg(nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:icatg")
WRITE (*,*) "Finished allocating arrays"
!-----------------------------------------------------------------------
! Initialize the allocated arrays to zero.
!-----------------------------------------------------------------------
x=0.0
y=0.0
z=0.0
hterain=0.0
mapfct =0.0
zp =0.0
j1 =0.0
j2 =0.0
j3 =0.0
aj3x=0.0
aj3y=0.0
aj3z=0.0
j3inv=0.0
u =0.0
v =0.0
w =0.0
wcont=0.0
ptprt=0.0
pprt =0.0
qv =0.0
qc =0.0
qr =0.0
qi =0.0
qs =0.0
qh =0.0
tke =0.0
ubar =0.0
vbar =0.0
ptbar=0.0
pbar =0.0
rhostr=0.0
qvbar=0.0
udteb=0.0
udtwb=0.0
vdtnb=0.0
vdtsb=0.0
pdteb=0.0
pdtwb=0.0
pdtnb=0.0
pdtsb=0.0
trigs1=0.0
trigs2=0.0
ifax1=0.0
ifax2=0.0
vwork1=0.0
vwork2=0.0
wsave1=0.0
wsave2=0.0
ppi =0.0
csndsq=0.0
radfrc=0.0
radsw =0.0
rnflx =0.0
soiltyp=0.0
stypfrct=0.0
vegtyp =0.0
lai =0.0
roufns =0.0
veg =0.0
tsfc =0.0
qvsfc =0.0
tsoil =0.0
wetsfc =0.0
wetdp =0.0
wetcanp=0.0
snowdpth=0.0
ptcumsrc=0.0
qcumsrc =0.0
w0avg =0.0
nca =0.0
kfraincv =0.0
cldefi = 0.0
xland = 0.0
bmjraincv = 0.0
raing =0.0
rainc =0.0
prcrate=0.0
usflx =0.0
vsflx =0.0
ptsflx=0.0
qvsflx=0.0
knt=0.0
rngsqi=0.0
wgtsum=0.0
zsum=0.0
sqsum=0.0
xs=0.0
ys=0.0
zs=0.0
latgr=0.0
longr=0.0
tem1=0.0
tem2=0.0
tem3=0.0
tem4=0.0
tem5=0.0
tem6=0.0
tem7=0.0
tem8=0.0
tem9=0.0
tem10=0.0
!wdt-williams update
ibegin=0
iend=0
anx=0.0
xcor=0.0
icatg=0.0
!
!-----------------------------------------------------------------------
!
! Set expected _squared_ background error for each variable
! This will depend on the particular background field used,
! season, age of forecast, etc.
!
!-----------------------------------------------------------------------
!
PRINT *, 'Reading ', TRIM(backerrfil)
OPEN(4,FILE=trim(backerrfil),STATUS='old')
READ(4,'(a80)') srcback
READ(4,'(a80)') header
WRITE(6,'(//,a,/a)') 'Background Std Error for',srcback
WRITE(6,'(1x,a)') &
' k hgt(m) u(m/s) v(m/s) pres(mb) temp(K) RH(%)'
DO ktab=1,nz_tab
READ(4,*,END=2) hqback(ktab), &
qback(3,ktab), &
qback(4,ktab), &
qback(5,ktab), &
qback(1,ktab), &
qback(2,ktab)
WRITE(6,'(1x,i4,f10.2,5f8.2)') ktab,hqback(ktab), &
(qback(ivar,ktab),ivar=1,5)
qback(3,ktab)=100.*qback(3,ktab)
qback(5,ktab)=0.01*qback(5,ktab)
END DO
2 nlvqback=ktab-1
CLOSE(4)
!
isrc=1
DO jsrc=1,nsrc_sng
IF(srcsng(jsrc) /= 'NULL') THEN
PRINT *, 'Reading ', TRIM(sngerrfil(jsrc))
OPEN(4,ERR=3,FILE=trim(sngerrfil(jsrc)),STATUS='old')
READ(4,'(a8)',ERR=3,END=3) rdsource
IF(rdsource /= srcsng(jsrc)) THEN
WRITE(6,'(a,i4,a,a,a,a,a)') &
' Mismatch of source names',jsrc, &
' read ',rdsource,' expected ',srcsng(jsrc)
STOP
END IF
READ(4,'(a80)',ERR=3,END=3) srcsng_full(isrc)
READ(4,*,ERR=3,END=3) qcmulsng(isrc)
READ(4,'(a80)',ERR=3,END=3) header
READ(4,*,ERR=3,END=3) &
qsrcsng(3,isrc), &
qsrcsng(4,isrc), &
qsrcsng(5,isrc), &
qsrcsng(1,isrc), &
qsrcsng(2,isrc)
WRITE(6,'(//,a,a,/a)') 'Single-level std error for ', &
srcsng(isrc),srcsng_full(isrc)
WRITE(6,'(a,f5.1)') ' QC multiplier: ',qcmulsng(isrc)
WRITE(6,'(1x,a)') &
' u (m/s) v (m/s) pres(mb) temp(K) RH(%)'
WRITE(6,'(1x,5f8.2)') (qsrcsng(ivar,isrc),ivar=1,5)
qsrcsng(3,isrc)=100.*qsrcsng(3,isrc)
qsrcsng(5,isrc)=0.01*qsrcsng(5,isrc)
CLOSE(4)
isrc=isrc+1
END IF
END DO
3 CONTINUE
!
isrc=1
DO jsrc=1,nsrc_ua
IF(srcua(jsrc) /= 'NULL') THEN
PRINT *, 'Reading ', TRIM(uaerrfil(jsrc))
OPEN(4,ERR=6,FILE=trim(uaerrfil(jsrc)),STATUS='old')
READ(4,'(a8)',ERR=6,END=6) rdsource
IF(rdsource /= srcua(jsrc)) THEN
WRITE(6,'(a,i4,a,a,a,a,a)') &
' Mismatch of source names',jsrc, &
' read ',rdsource,' expected ',srcua(jsrc)
STOP
END IF
READ(4,'(a80)',ERR=6,END=6) srcua_full(isrc)
READ(4,*,ERR=6,END=6) qcmulua(isrc)
READ(4,'(a80)',ERR=6,END=6) header
WRITE(6,'(//,a,a,/a)') 'UA data std error for ', &
srcua(isrc),srcua_full(isrc)
WRITE(6,'(a,f5.1)') ' QC multiplier: ',qcmulua(isrc)
WRITE(6,'(1x,a)') &
' k hgt(m) u(m/s) v(m/s) pres(mb) temp(K) RH(%)'
DO ktab=1,nz_tab
READ(4,*,ERR=5,END=5) huaqsrc(ktab,isrc), &
qsrcua(3,ktab,isrc), &
qsrcua(4,ktab,isrc), &
qsrcua(5,ktab,isrc), &
qsrcua(1,ktab,isrc), &
qsrcua(2,ktab,isrc)
WRITE(6,'(1x,i4,f10.2,5f8.2)') ktab,huaqsrc(ktab,isrc), &
(qsrcua(ivar,ktab,isrc),ivar=1,5)
qsrcua(3,ktab,isrc)=100.*qsrcua(3,ktab,isrc)
qsrcua(5,ktab,isrc)=0.01*qsrcua(5,ktab,isrc)
END DO
5 nlvuatab(isrc)=ktab-1
CLOSE(4)
isrc=isrc+1
END IF
END DO
6 CONTINUE
!
isrc=1
DO jsrc=1,nsrc_rad
IF(srcrad(jsrc) /= 'NULL') THEN
PRINT *, 'Reading ', TRIM(raderrfil(jsrc))
OPEN(4,ERR=7,FILE=trim(raderrfil(jsrc)),STATUS='old')
READ(4,'(a8)',ERR=7,END=7) rdsource
IF(rdsource /= srcrad(jsrc)) THEN
WRITE(6,'(a,i4,a,a,a,a,a)') &
' Mismatch of source names',jsrc, &
' read ',rdsource,' expected ',srcrad(jsrc)
STOP
END IF
READ(4,'(a80)',ERR=7,END=7) srcrad_full(isrc)
READ(4,*,ERR=7,END=7) qcmulrad(isrc)
READ(4,'(a80)',ERR=7,END=7) header
READ(4,*,ERR=7,END=7) &
qsrcrad(1,isrc), &
qsrcrad(2,isrc), &
qsrcrad(3,isrc)
WRITE(6,'(//,a,a,/a)') 'Radar data std error for ', &
srcrad(isrc),srcrad_full(isrc)
WRITE(6,'(a,f5.1)') ' QC multiplier: ',qcmulrad(isrc)
WRITE(6,'(1x,a)') &
' ref(dBz) Vrad(m/s) SpWid(m/s)'
WRITE(6,'(1x,4f8.2)') &
(qsrcrad(ivar,isrc),ivar=1,nvar_radin)
CLOSE(4)
isrc=isrc+1
END IF
END DO
7 CONTINUE
!
isrc=1
DO jsrc=1,nsrc_ret
IF(srcret(jsrc) /= 'NULL') THEN
PRINT *, 'Reading ', TRIM(reterrfil(jsrc))
OPEN(4,FILE=trim(reterrfil(jsrc)),ERR=10,STATUS='old')
READ(4,'(a8)',ERR=10,END=10) rdsource
IF(rdsource /= srcret(jsrc)) THEN
WRITE(6,'(a,i4,a,a,a,a,a)') &
' Mismatch of source names',jsrc, &
' read ',rdsource,' expected ',srcret(jsrc)
STOP
END IF
READ(4,'(a80)',ERR=10,END=10) srcret_full(isrc)
READ(4,*,ERR=10,END=10) qcmulret(isrc)
READ(4,'(a80)',ERR=10,END=10) header
WRITE(6,'(//,a,a,/a)') 'Retrieval std error for ', &
srcret(isrc),srcret_full(isrc)
WRITE(6,'(a,f5.1)') ' QC multiplier: ',qcmulrad(isrc)
WRITE(6,'(1x,a)') &
' k hgt(m) u(m/s) v(m/s) pres(mb) temp(K) RH(%)'
DO ktab=1,nz_tab
READ(4,*,END=9,ERR=9) hrtqsrc(ktab,isrc), &
qsrcret(3,ktab,isrc), &
qsrcret(4,ktab,isrc), &
qsrcret(5,ktab,isrc), &
qsrcret(1,ktab,isrc), &
qsrcret(2,ktab,isrc)
WRITE(6,'(1x,i4,f10.2,5f8.2)') ktab,hrtqsrc(ktab,isrc), &
(qsrcret(ivar,ktab,isrc),ivar=1,5)
qsrcret(3,ktab,isrc)=100.*qsrcret(3,ktab,isrc)
qsrcret(5,ktab,isrc)=0.01*qsrcret(5,ktab,isrc)
END DO
9 nlvrttab(isrc)=ktab-1
CLOSE(4)
isrc=isrc+1
END IF
END DO
10 CONTINUE
!
!-----------------------------------------------------------------------
!
! Change standard error to standard error variance by squaring
! Calculate quality control thresholds
!
!-----------------------------------------------------------------------
!
DO ktab=1,nlvqback
DO ivar=1,nvar_anx
qback(ivar,ktab)=qback(ivar,ktab)*qback(ivar,ktab)
END DO
END DO
!
DO isrc=1,nsrc_sng
DO ivar=1,nvar_anx
qsrcsng(ivar,isrc)=qsrcsng(ivar,isrc)*qsrcsng(ivar,isrc)
qcthrsng(ivar,isrc)=qcmulsng(isrc)* &
SQRT(qback(ivar,1)+qsrcsng(ivar,isrc))
END DO
END DO
!
DO isrc=1,nsrc_ua
DO ivar=1,nvar_anx
qsrcmax=0.
DO ktab=1,nlvuatab(isrc)
qsrcua(ivar,ktab,isrc) = &
qsrcua(ivar,ktab,isrc)*qsrcua(ivar,ktab,isrc)
qsrcmax=AMAX1(qsrcmax,qsrcua(ivar,ktab,isrc))
END DO
qcthrua(ivar,isrc)=qcmulua(isrc)* &
SQRT(qback(ivar,1)+qsrcmax)
END DO
END DO
!
DO isrc=1,nsrc_rad
DO ivar=1,nvar_radin
qsrcrad(ivar,isrc) = &
qsrcrad(ivar,isrc)*qsrcrad(ivar,isrc)
qcthrrad(ivar,isrc)=qcmulrad(isrc)* &
SQRT(qback(ivar,1)+qsrcrad(ivar,isrc))
END DO
END DO
!
DO isrc=1,nsrc_ret
DO ivar=1,nvar_anx
qsrcmax=0.
DO ktab=1,nlvrttab(isrc)
qsrcret(ivar,ktab,isrc) = &
qsrcret(ivar,ktab,isrc)*qsrcret(ivar,ktab,isrc)
qsrcmax=AMAX1(qsrcmax,qsrcret(ivar,ktab,isrc))
END DO
qcthrret(ivar,isrc)=qcmulret(isrc)* &
SQRT(qback(ivar,1)+qsrcmax)
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Set cross-correlations between numerical categories.
! Roughly 1=clear,2=some evidence of outflow,3=precip,4=conv precip
! This could be read in as a table.
!
!-----------------------------------------------------------------------
!
DO j=1,ncat
DO i=1,ncat
xcor(i,j)=1.0-(IABS(i-j))*0.25
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Initialize grids, forming first guess fields based on
! the model initial option specified in the ARPS init file.
!
!-----------------------------------------------------------------------
!
CALL initgrdvar
(nx,ny,nz,1,nstyps,1, &
x,y,z,zp,hterain,mapfct, &
j1,j2,j3,aj3x,aj3y,aj3z,j3inv, &
u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke, &
udteb, udtwb, vdtnb, vdtsb, &
pdteb,pdtwb ,pdtnb ,pdtsb, &
trigs1,trigs2,ifax1,ifax2, &
wsave1,wsave2,vwork1,vwork2, &
ubar,vbar,ptbar,pbar,tem10,tem10, &
rhostr,tem10,qvbar,ppi,csndsq, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,qvsfc, &
ptcumsrc,qcumsrc,w0avg,nca,kfraincv, &
cldefi,xland,bmjraincv, &
raing,rainc,prcrate,exbcbuf, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
tem1,tem2,tem3,tem4,tem5, &
tem6,tem7,tem8,tem9)
!
!-----------------------------------------------------------------------
!
! Deallocate some arrays needed only for initgrdvar
!
!-----------------------------------------------------------------------
!
DEALLOCATE(ptcumsrc)
DEALLOCATE(qcumsrc)
DEALLOCATE(w0avg)
DEALLOCATE(kfraincv)
DEALLOCATE(nca)
DEALLOCATE(cldefi)
DEALLOCATE(xland)
!
!-----------------------------------------------------------------------
!
! Set location of scalar fields.
!
!-----------------------------------------------------------------------
!
CALL xytoll
(nx,ny,x,y,latgr,longr)
DO i=1,nx-1
xs(i)=0.5*(x(i)+x(i+1))
END DO
xs(nx)=2.0*xs(nx-1)-xs(nx-2)
DO j=1,ny-1
ys(j)=0.5*(y(j)+y(j+1))
END DO
ys(ny)=2.0*ys(ny-1)-ys(ny-2)
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
zs(i,j,k)=0.5*(zp(i,j,k)+zp(i,j,k+1))
END DO
END DO
END DO
DO i=1,nx-1
DO j=1,ny-1
zs(i,j,nz)=2.*zs(i,j,nz-1)-zs(i,j,nz-2)
END DO
END DO
CALL ctim2abss
( year,month,day,hour,minute,second, i4timanx)
!
!-----------------------------------------------------------------------
!
! Identify the background field correlation category for
! each 2-d point. This is based on precip rate, cumulus
! parameterization switch, and surface relative humidity.
!
!-----------------------------------------------------------------------
!
CALL setcat
(nx,ny,nz,ccatopt,zs, &
ptprt,pprt,qv,qc,qr,qi,qs,qh, &
ptbar,pbar, &
prcrate,icatg)
!
!-----------------------------------------------------------------------
!
! Load "background fields" for analysis
! Note, analysis is done at scalar points
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
anx(i,j,k,1)=0.5*(u(i,j,k)+u(i+1,j,k))
anx(i,j,k,2)=0.5*(v(i,j,k)+v(i,j+1,k))
anx(i,j,k,3)=pbar(i,j,k)+pprt(i,j,k)
anx(i,j,k,4)=ptbar(i,j,k)+ptprt(i,j,k)
anx(i,j,k,5)=qv(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Save background fields on staggered grid in the increment
! arrays.
!
!-----------------------------------------------------------------------
!
IF(incrdmp > 0) THEN
CALL incrsave
(nx,ny,nz,nxndg,nyndg,nzndg, &
u,v,w,pprt,ptprt,qv, &
qc,qr,qi,qs,qh, &
uincr,vincr,wincr,pincr,ptincr,qvincr, &
qcincr,qrincr,qiincr,qsincr,qhincr, &
istatus)
END IF
!
!-----------------------------------------------------------------------
!
! Filename creation
!
!-----------------------------------------------------------------------
!
WRITE(stats,'(a,a)') runname(1:lfnkey),'.lst'
CALL getunit
(iwstat)
PRINT *, 'Writing ', TRIM(stats)
OPEN(iwstat,IOSTAT=ios,FILE=trim(stats),STATUS='unknown')
IF(ios /= 0) iwstat=0
maxsuf=LEN(suffix)
mxroot=LEN(froot)
nobsng=0
DO ifile=1,nsngfil
PRINT *, 'Processing file ', ifile, ' ', TRIM(sngfname(ifile))
lenfnm=LEN(sngfname(ifile))
CALL strlnth
(sngfname(ifile),lenfnm)
CALL exsufx
(sngfname(ifile),lenfnm,suffix,maxsuf,dotloc,lensuf)
IF(lensuf == 3. .AND. suffix(1:3) == 'lso') THEN
CALL exfroot
(sngfname(ifile),lenfnm,froot,mxroot,lenroot)
WRITE(qclist,'(a,a)') froot(1:lenroot),'.sqc'
WRITE(qcsave,'(a,a)') froot(1:lenroot),'.lsq'
!
!-----------------------------------------------------------------------
!
! Open the files for listing QC info
! To suppress listing set the unit numbers to zero
!
!-----------------------------------------------------------------------
!
CALL getunit
(iqclist)
PRINT *, 'Opening qclist: ', TRIM(qclist)
OPEN(iqclist,IOSTAT=ios,FILE=trim(qclist),STATUS='unknown')
IF(ios /= 0) iqclist=0
CALL getunit
(iqcsave)
PRINT *, 'Opening qclist: ', TRIM(qcsave)
OPEN(iqcsave,IOSTAT=ios,FILE=trim(qcsave),STATUS='unknown')
IF(ios /= 0) iqcsave=0
!
!-----------------------------------------------------------------------
!
! Read surface data, QC and convert units.
!
!-----------------------------------------------------------------------
!
sprkm=0.001*sprdist
rngsq=sfcqcrng*sfcqcrng
PRINT *, 'Calling prepsfcobs'
CALL prepsfcobs(ntime,mx_sng, &
nvar_sng,nvar_anx,nsrc_sng,nobsng, &
sngfname(ifile),sngtmchk(ifile),blackfil, &
var_sfc,var_anx,srcsng,qsrcsng, &
rmiss,iqspr,sprkm,nobsrd,timesng, &
stnsng,csrcsng,isrcsng,latsng,lonsng,hgtsng,xsng,ysng, &
wx,kloud,idp3,store_emv,store_hgt,store_amt, &
obrdsng,obsng,qualrdsng,qobsng,qualsng, &
ival,climin,climax, &
rngsq,klim,wlim,qcthrsng, &
knt,wgtsum,zsum,sqsum,iqclist,iqcsave,jstatus)
IF(iqclist /= 0) THEN
CLOSE(iqclist)
CALL retunit(iqclist)
END IF
IF(iqcsave /= 0) THEN
CLOSE(iqcsave)
CALL retunit(iqcsave)
END IF
ELSE
!
!-----------------------------------------------------------------------
!
! Read other single-level data and convert units.
!
!-----------------------------------------------------------------------
!
PRINT *, 'Calling prepsglobs'
CALL prepsglobs
(mx_sng,ntime,srcsng,nsrc_sng, &
sngfname(ifile),stnsng,latsng,lonsng,xsng,ysng, &
hgtsng,obsng,qobsng,qualsng,isrcsng,qsrcsng, &
csrcsng,nx,ny,nz,nvar_anx,anx,xs,ys,zp, &
tem1,tem2,tem3,tem4,tem5,tem6, &
rmiss,nobsng,istatus)
END IF
END DO
!
!-----------------------------------------------------------------------
!
! Calculate initial obs differences for single-level data
!
!-----------------------------------------------------------------------
!
PRINT *, 'Calling grdtosng'
CALL grdtosng
(nx,ny,nz,nz_tab,mx_sng,nvar_anx,nobsng, &
xs,ys,zp,icatg,anx,qback,hqback,nlvqback, &
tem1,tem2,tem3,tem4,tem5,tem6, &
stnsng,isrcsng,icatsng,hgtsng,xsng,ysng, &
obsng,qobsng,qualsng, &
odifsng,oanxsng,thesng)
!
!-----------------------------------------------------------------------
!
! Read upper-air data, QC and convert units
!
!-----------------------------------------------------------------------
!
PRINT *, 'Calling prepuaobs'
CALL prepuaobs
(nvar_anx,nz_ua,mx_ua,nz_tab,nsrc_ua, &
mx_ua_file,nuafil,uafname,srcua, &
stnua,elevua,xua,yua,hgtua,obsua, &
qsrcua,huaqsrc,nlvuatab, &
qobsua,qualua,isrcua,nlevsua, &
rmiss,nobsua,istatus)
!
!-----------------------------------------------------------------------
!
! Calculate initial obs differences for upper-air data
!
!-----------------------------------------------------------------------
!
PRINT *, 'Calling grdtoua'
CALL grdtoua
(nx,ny,nz,nz_tab,nz_ua,mx_ua,nvar_anx,nobsua, &
xs,ys,zp,anx,qback,hqback,nlvqback, &
tem1,tem2,tem3,tem4,tem5,tem6, &
stnua,isrcua,elevua,xua,yua,hgtua, &
obsua,qobsua,qualua,nlevsua, &
odifua,oanxua,theua)
!
!-----------------------------------------------------------------------
!
! Read radar data, unfold and convert units
!
!-----------------------------------------------------------------------
!
PRINT *, 'Calling prepradar'
CALL prepradar
(nx,ny,nz,nz_tab,nvar_anx,nvar_radin, &
nvar_rad,mx_rad,nsrc_rad,nz_rdr,mx_colrad,mx_pass, &
raduvobs,radrhobs,radistride,radkstride, &
iuserad,npass,refrh,rhradobs, &
xs,ys,zs,hterain,anx,qback,hqback,nlvqback, &
nradfil,radfname,srcrad,isrcrad,qsrcrad,qcthrrad, &
stnrad,latrad,lonrad,elvrad, &
latradc,lonradc,xradc,yradc,irad,nlevrad, &
distrad,uazmrad,vazmrad,hgtradc,theradc, &
obsrad,oanxrad,odifrad,qobsrad,qualrad, &
ncolrad,istatus,tem1,tem2,tem3,tem4,tem5,tem6)
!
!-----------------------------------------------------------------------
!
! Read retrieval data.
!
!-----------------------------------------------------------------------
!
PRINT *, 'Calling prepretr'
CALL prepretr
(nx,ny,nz,nvar_anx, &
nz_ret,mx_ret,mx_colret,nz_tab,nsrc_ret, &
nretfil,retfname, &
isrcret,srcret,nlvrttab,qsrcret,hrtqsrc, &
stnret,latret,lonret,elvret, &
latretc,lonretc,xretc,yretc,iret,nlevret, &
hgtretc,obsret,qrret,qobsret,qualret, &
rmiss,ncolret,tem1,istatus)
PRINT *, 'Calling grdtoret'
CALL grdtoret
(nx,ny,nz,nz_tab, &
nz_ret,mx_ret,mx_colret,nvar_anx,ncolret, &
xs,ys,zp,anx,qback,hqback,nlvqback, &
tem1,tem2,tem3,tem4,tem5,tem6, &
stnret,iret,xretc,yretc,hgtretc, &
obsret,qobsret,qualret,nlevret, &
odifret,oanxret,theretc)
!
!-----------------------------------------------------------------------
!
! Quality-control observation differences
!
!-----------------------------------------------------------------------
!
PRINT *, 'Calling qcdiff'
CALL qcdiff
(nvar_anx,nvar_rad,nvar_radin,mx_sng,nsrc_sng, &
nz_ua,mx_ua,nsrc_ua, &
nz_rdr,mx_rad,mx_colrad,nsrc_rad, &
nz_ret,mx_ret,mx_colret,nsrc_ret, &
nobsng,nobsua,ncolrad,ncolret,var_anx, &
stnsng,isrcsng,hgtsng,obsng,oanxsng,odifsng, &
qcthrsng,qualsng, &
stnua,isrcua,hgtua,obsua,oanxua,odifua, &
qcthrua,qualua,nlevsua, &
stnrad,irad,isrcrad,hgtradc,obsrad,odifrad, &
qcthrrad,qualrad,nlevrad, &
stnret,iret,isrcret,hgtretc, &
obsret,oanxret,odifret, &
qcthrret,qualret,nlevret, &
wgtsum,zsum)
!
!-----------------------------------------------------------------------
!
! Analyze data
!
!-----------------------------------------------------------------------
!
!wdt
CALL FLUSH
(06)
PRINT *, 'Calling anxiter'
CALL anxiter
(nx,ny,nz, &
nvar_anx,nvar_radin,nvar_rad,nz_ua,nz_rdr,nz_ret, &
mx_sng,mx_ua,mx_rad,mx_colrad,mx_ret,mx_colret, &
nsrc_sng,nsrc_ua,nsrc_rad,nsrc_ret,ncat, &
mx_pass,npass,iwstat,xs,ys,zs,icatg,xcor,var_anx, &
xsng,ysng,hgtsng,thesng, &
obsng,odifsng,qobsng,qualsng,isrcsng,icatsng,nobsng, &
xua,yua,hgtua,theua, &
obsua,odifua,qobsua,qualua,isrcua,nlevsua,nobsua, &
elvrad,xradc,yradc, &
distrad,uazmrad,vazmrad,hgtradc,theradc, &
obsrad,odifrad,qobsrad,qualrad, &
irad,isrcrad,nlevrad,ncolrad, &
xretc,yretc,hgtretc,theretc, &
obsret,odifret,qobsret,qualret, &
iret,isrcret,nlevret,ncolret, &
srcsng,srcua,srcrad,srcret, &
ianxtyp,iusesng,iuseua,iuserad,iuseret, &
xyrange,kpvrsq,wlim,zrange,zwlim, &
thrng,rngsqi,knt,wgtsum,zsum, &
corsng,corua,corrad,corret, &
oanxsng,oanxua,oanxrad,oanxret, &
anx,tem1,tem2,tem3,ibegin,iend,istatus)
!
!wdt
!-----------------------------------------------------------------------
!
! Write to ARPSControl status file
!
!-----------------------------------------------------------------------
! Only u (index 1) is checked. Eventually should check all
! Only 1 level of UA data is checked.
!
! 10 Good
! 100 Good
! 200 Good - superob
! -9 Outside domain
! -99 Initializing value
! -199 Rejected in quality test
!
! > 0 : used++ : in file, in domain, good point
! = -199 : rej++ : in file, in domain, bad point
! = -9 : : in file, outside domain
! = 0 : : not in file
!
! Good data has isrc* > 0 and qual* > 0.
! single-point: SA (sao), MESO (okmeso)
WRITE (status_file, '(2A)') runname(1:lfnkey),'.adasstat'
PRINT *, 'Creating ', TRIM(status_file)
OPEN (UNIT=09, FILE=status_file)
n_used_sng = 0
DO jsta=1,nobsng
IF (isrcsng(jsta) .GT. 0 .AND. qualsng(1,jsta) .GT. 0) THEN
n_used_sng(isrcsng(jsta)) = n_used_sng(isrcsng(jsta)) + 1
ELSE
END IF
END DO
DO jsrc=1,nsrc_sng
IF (srcsng(jsrc)(:2) .EQ. 'SA') THEN
WRITE (09,9901) '$sao_n_used = ', n_used_sng(jsrc)
ELSE IF (srcsng(jsrc)(:4) .EQ. 'MESO') THEN
WRITE (09,9901) '$meso_n_used = ', n_used_sng(jsrc)
ELSE IF (srcsng(jsrc)(:5) .EQ. 'MDCRS') THEN
WRITE (09,9901) '$mdcrs_n_used = ', n_used_sng(jsrc)
!wdt update 2002-03-05 from kwthomas -- added buoy
ELSE IF (srcsng(jsrc)(:4) .EQ. 'BUOY') THEN
WRITE (09,9901) '$buoy_n_used = ', n_used_sng(jsrc)
ELSE IF (srcsng(jsrc)(:4) .NE. 'NULL') THEN
WRITE (09,'(3A)') '# Unable to process srcsng: ', &
srcsng(jsrc), stnsng(jsta,1)
END IF
END DO
! upper-air: NWS RAOB (raob), WPDN PRO (pro)
n_used_ua = 0
DO jsta=1,nobsua
IF ((isrcua(jsta) <= 0) .or. (isrcua(jsta) > nsrc_ua)) THEN
WRITE (09,*) '# Unable to process: irsrcua=', &
isrcua(jsta), ' ', TRIM(stnua(jsta))
ELSE IF (srcua(isrcua(jsta)) .EQ. 'NWS RAOB') THEN
WRITE (09,9902) 'raob', TRIM(stnua(jsta)), nlevsua(jsta)
n_used_ua(isrcua(jsta)) = n_used_ua(isrcua(jsta)) + 1
ELSE IF (srcua(isrcua(jsta)) .EQ. 'WPDN PRO') THEN
WRITE (09,9902) 'pro', TRIM(stnua(jsta)), nlevsua(jsta)
n_used_ua(isrcua(jsta)) = n_used_ua(isrcua(jsta)) + 1
ELSE IF (srcua(isrcua(jsta))(:4) .NE. 'NULL') THEN
WRITE (09,'(3A)') '# Unable to process: ', &
srcua(isrcua(jsta)), TRIM(stnua(jsta))
ENDIF
END DO
!depricated:
! DO jsta=1,nobsua
!! criterion: 3 good u obs out of the first 5 levels
! IF (isrcua(jsta) .GT. 0) THEN
! ngoodlev = 0
!
! DO klev=1,5 !nlevsua(jsta)
! IF (qualua(1,klev,jsta) .GT. 0) ngoodlev = ngoodlev + 1
! END DO
!
! IF (ngoodlev .GE. 3) THEN
! IF (srcua(isrcua(jsta)) .EQ. 'NWS RAOB') THEN
! WRITE (09,9902) 'raob', TRIM(stnua(jsta)), nlevsua(jsta)
! ELSE IF (srcua(isrcua(jsta)) .EQ. 'WPDN PRO') THEN
! WRITE (09,9902) 'pro', TRIM(stnua(jsta)), nlevsua(jsta)
! ELSE IF (srcua(jsrc)(:4) .NE. 'NULL') THEN
! WRITE (09,'(3A)') '# Unable to process: ', &
! srcua(jsrc), TRIM(stnua(jsta))
! END IF
! n_used_ua(isrcua(jsta)) = n_used_ua(isrcua(jsta)) + 1
! END IF
!
! END IF ! good isrcua
! END DO
DO jsrc=1,nsrc_ua
IF (srcua(jsrc)(:8) .EQ. 'NWS RAOB') THEN
WRITE (09,9901) '$raob_n_used = ', n_used_ua(jsrc)
ELSE IF (srcua(jsrc)(:8) .EQ. 'WPDN PRO') THEN
WRITE (09,9901) '$pro_n_used = ', n_used_ua(jsrc)
END IF
END DO
9902 FORMAT ('${', A, '}{stn_levels}{"', A, '"} = ', I6, ';')
9901 FORMAT (A,I6,';')
! wdt
CALL FLUSH
(06)
CLOSE (09)
!-----------------------------------------------------------------------
!
! Use reflectivity to establish cloud water
!
!-----------------------------------------------------------------------
!
IF( radcldopt > 0 ) CALL radmcro
(nx,ny,nz, &
mx_rad,nsrc_rad,nvar_radin,nz_rdr,mx_colrad,mx_pass, &
radistride,radkstride,iuserad,npass, &
nradfil,radfname,srcrad,isrcrad,stnrad, &
latrad,lonrad,elvrad,latradc,lonradc,irad, &
nlevrad,xradc,yradc,hgtradc,obsrad,ncolrad, &
radqvopt,radqcopt,radqropt,radptopt, &
refsat,rhrad, &
refcld,cldrad,ceilopt,ceilmin,dzfill, &
refrain,radsetrat,radreflim,radptgain, &
xs,ys,zs,zp,j3, &
anx(1,1,1,3),anx(1,1,1,4),anx(1,1,1,5), &
ptbar,qvbar,rhostr, &
qc,qr,qi,qs,qh, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7)
!
!-----------------------------------------------------------------------
!
! Transfer analysed fields to ARPS variables for writing.
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
pprt(i,j,k)=anx(i,j,k,3)-pbar(i,j,k)
ptprt(i,j,k)=anx(i,j,k,4)-ptbar(i,j,k)
tgrid =anx(i,j,k,4) * (anx(i,j,k,3)/p0) ** rddcp
qvsat=f_qvsat
( anx(i,j,k,3),tgrid)
qvmin=rhmin*qvsat
qv(i,j,k)=MAX(qvmin,MIN(anx(i,j,k,5),qvsat))
END DO
END DO
END DO
IF(spradopt == 1) THEN
!-----------------------------------------------------------------------
!
! Check for superadiabatic layers in each column.
! For spradopt=1,
! Where superadiabatic levels are found, make the next level down
! cooler, applying the minimum potential temperature lapse rate,
! ptlapse.
!
!-----------------------------------------------------------------------
!
DO j=1,ny-1
DO i=1,nx-1
ptprt(i,j,nz-1)=anx(i,j,nz-1,4)-ptbar(i,j,nz-1)
DO k=nz-2,2,-1
anx(i,j,k,4)=MIN(anx(i,j,k,4), &
(anx(i,j,k+1,4)-ptlapse*(zs(i,j,k+1)-zs(i,j,k))))
ptprt(i,j,k)=anx(i,j,k,4)-ptbar(i,j,k)
END DO
END DO
END DO
ELSE IF(spradopt == 2) THEN
!
!-----------------------------------------------------------------------
!
! For spradopt=2,
! Where superadiabatic levels are found, make the next level up
! warmer, applying the minimum potential temperature lapse rate,
! ptlapse.
!
!-----------------------------------------------------------------------
!
DO j=1,ny-1
DO i=1,nx-1
ptprt(i,j,2)=anx(i,j,2,4)-ptbar(i,j,2)
DO k=3,nz-1
anx(i,j,k,4)=MAX(anx(i,j,k,4), &
(anx(i,j,k-1,4)+ptlapse*(zs(i,j,k)-zs(i,j,k-1))))
ptprt(i,j,k)=anx(i,j,k,4)-ptbar(i,j,k)
END DO
END DO
END DO
END IF
!
DO k=1,nz-1
DO j=1,ny-1
DO i=2,nx-1
u(i,j,k)=0.5*(anx(i,j,k,1)+anx(i-1,j,k,1))
END DO
END DO
END DO
DO k=1,nz-1
DO j=1,ny-1
u(1,j,k)=u(2,j,k)
u(nx,j,k)=u(nx-1,j,k)
END DO
END DO
!
DO k=1,nz-1
DO j=2,ny-1
DO i=1,nx-1
v(i,j,k)=0.5*(anx(i,j,k,2)+anx(i,j-1,k,2))
END DO
END DO
END DO
DO k=1,nz-1
DO i=1,nx-1
v(i,1,k)=v(i,2,k)
v(i,ny,k)=v(i,ny-1,k)
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Adjust wind fields to minimize 3-d divergence
! This is the first run-through for "w" in order to obtain
! a good estimate of w for the cloud analysis. Note here wndadj
! is set to "2" for integration of 2-d divergence for "w" and
! the cloud_w opt is set to 0.
!
!-----------------------------------------------------------------------
!
CALL adjuvw
( nx,ny,nz,u,v,w, &
wcont,ptprt,ptbar, &
zp,j1,j2,j3,aj3z,mapfct,rhostr,w_cld, &
2,obropt,obrzero,0, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8)
!
!-----------------------------------------------------------------------
!
! Enforce vertical w boundary conditions
!
!-----------------------------------------------------------------------
!
CALL rhouvw
(nx,ny,nz,rhostr,tem1,tem2,tem3)
CALL vbcw
(nx,ny,nz,w,wcont,tbc,bbc,u,v, &
rhostr,tem1,tem2,tem3, &
j1,j2,j3)
!
!-----------------------------------------------------------------------
!
! Free-up some memory.
!
!-----------------------------------------------------------------------
!
DEALLOCATE(anx)
!
!-----------------------------------------------------------------------
!
! Call cloud analysis.
!
!-----------------------------------------------------------------------
!
IF( cloudopt > 0 ) THEN
ALLOCATE(ref_mos_3d(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:ref_mos_3d")
ref_mos_3d=0.
ALLOCATE(w_cld(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:w_cld")
w_cld=0.
ALLOCATE(cldpcp_type_3d(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:cldpcp_type_3d")
cldpcp_type_3d=0
ALLOCATE(icing_index_3d(nx,ny,nz),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:icing_index_3d")
icing_index_3d=0
ALLOCATE(l_mask_pcptype(nx,ny),STAT=istatus)
CALL check_alloc_status
(istatus, "adas:l_mask_pcptype")
l_mask_pcptype=.false.
CALL cmpclddrv
(nx,ny,nz,i4timanx, &
xs,ys,zs,j3,hterain,latgr,longr, &
nobsng,stnsng,csrcsng,xsng,ysng, &
timesng,latsng,lonsng,hgtsng, &
kloud,store_amt,store_hgt, &
stnrad,isrcrad,latrad,lonrad,elvrad,ixrad,jyrad, &
pprt,ptprt,qv,qc,qr,qi,qs,qh,w, &
pbar,ptbar,qvbar,rhostr, &
ref_mos_3d,w_cld,cldpcp_type_3d, &
icing_index_3d,l_mask_pcptype, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,tem9,tem10)
END IF
!
IF(retqropt > 0 .AND. ncolret > 0) CALL retmcro
(nx,ny,nz, &
mx_ret,nsrc_ret,nvar_anx,nz_ret,mx_colret, &
srcret,isrcret,iret,nlevret, &
xretc,yretc,hgtretc,qrret,ncolret, &
dzfill, &
xs,ys,zs,qr)
!
!-----------------------------------------------------------------------
!
! Create rhobar from rhostr
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem6(i,j,k)=rhostr(i,j,k)/j3(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Calculate and store the sound wave speed squared in csndsq.
! Use original definition of sound speed.
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
csndsq(i,j,k)= cpdcv*pbar(i,j,k)/tem6(i,j,k)
END DO
END DO
END DO
!
IF(hydradj == 1 .OR. hydradj == 2) THEN
pconst=0.5*g*dz
!
!-----------------------------------------------------------------------
!
! Create thermal buoyancy at each scalar point,
! which is stored in tem2
!
!-----------------------------------------------------------------------
!
DO k=1,nz
DO j=1,ny
DO i=1,nx
qvprt=qv(i,j,k)-qvbar(i,j,k)
qtot=qc(i,j,k)+qr(i,j,k)+ &
qi(i,j,k)+qs(i,j,k)+qh(i,j,k)
tem2(i,j,k)=j3(i,j,k)*tem6(i,j,k)* &
g*( (ptprt(i,j,k)/ptbar(i,j,k)) + &
(qvprt/(rddrv+qvbar(i,j,k))) - &
((qvprt+qtot)/(1.+qvbar(i,j,k))) )
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Average thermal buoyancy to w points
! which is stored in tem3
!
!-----------------------------------------------------------------------
!
CALL avgsw
(tem2,nx,ny,nz,1,nx,1,ny, tem3)
IF(hydradj == 1) THEN
DO i=1,nx
DO j=1,ny
tem1(i,j,1)=pprt(i,j,1)
DO k=2,nz-2
tem1(i,j,k)= &
( (1. - (pconst*j3(i,j,k-1)/csndsq(i,j,k-1)) )* &
tem1(i,j,k-1) + &
dz*tem3(i,j,k) ) / &
(1. + (pconst*j3(i,j,k)/csndsq(i,j,k)) )
IF(j == 26 .AND. MOD(i,10) == 0) THEN
IF(k == nz-2) PRINT *, ' Point i= ',i, ' j=26'
PRINT *, ' k,pprt,tem1,diff =', &
k,pprt(i,j,k),tem1(i,j,k), &
(tem1(i,j,k)-pprt(i,j,k))
END IF
pprt(i,j,k)=tem1(i,j,k)
END DO
pprt(i,j,nz-1)=pprt(i,j,nz-2)
END DO
END DO
ELSE IF(hydradj == 2) THEN
DO i=1,nx
DO j=1,ny
tem1(i,j,nz-1)=pprt(i,j,nz-1)
DO k=nz-2,2,-1
tem1(i,j,k)= &
( (1.+ (pconst*j3(i,j,k+1)/csndsq(i,j,k+1)) )* &
tem1(i,j,k+1) - &
dz*tem3(i,j,k+1) ) / &
(1.- (pconst*j3(i,j,k )/csndsq(i,j,k )) )
IF(j == 26 .AND. MOD(i,10) == 0) THEN
IF(k == nz-2) PRINT *, ' Point i= ',i, ' j=26'
PRINT *, ' k,pprt,tem1,diff =', &
k,pprt(i,j,k),tem1(i,j,k), &
(tem1(i,j,k)-pprt(i,j,k))
END IF
pprt(i,j,k)=tem1(i,j,k)
END DO
pprt(i,j,1)=pprt(i,j,2)
END DO
END DO
END IF
ELSE IF (hydradj == 3) THEN
!
!-----------------------------------------------------------------------
!
! Calculate total pressure, store in tem1.
! Calculate virtual temperature, store in tem2.
!
!-----------------------------------------------------------------------
!
DO k=1,nz
DO j=1,ny
DO i=1,nx
tem1(i,j,k) = pbar(i,j,k)+pprt(i,j,k)
temp = (ptbar(i,j,k)+ptprt(i,j,k))* &
((tem1(i,j,k)/p0)**rddcp)
tem2(i,j,k) = temp*(1.0+rvdrd*qv(i,j,k))/ &
(1.0+qv(i,j,k))
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Integrate hydrostatic equation, begining at k=2
!
!-----------------------------------------------------------------------
!
pconst=-g/rd
DO k=3,nz-1
DO j=1,ny-1
DO i=1,nx-1
tvbar=0.5*(tem2(i,j,k)+tem2(i,j,k-1))
tem1(i,j,k)=tem1(i,j,k-1)* &
EXP(pconst*(zs(i,j,k)-zs(i,j,k-1))/tvbar)
pprt(i,j,k)=tem1(i,j,k)-pbar(i,j,k)
END DO
END DO
END DO
DO j=1,ny-1
DO i=1,nx-1
tvbar=tem2(i,j,2)
tem1(i,j,1)=tem1(i,j,2)* &
EXP(pconst*(zs(i,j,1)-zs(i,j,2))/tvbar)
pprt(i,j,1)=tem1(i,j,1)-pbar(i,j,1)
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Adjust wind fields to minimize 3-d divergence
!
!-----------------------------------------------------------------------
!
CALL adjuvw
( nx,ny,nz,u,v,w, &
wcont,ptprt,ptbar, &
zp,j1,j2,j3,aj3z,mapfct,rhostr,w_cld, &
wndadj,obropt,obrzero,cldwopt, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8)
!
!-----------------------------------------------------------------------
!
! Enforce vertical w boundary conditions
!
!-----------------------------------------------------------------------
!
CALL rhouvw
(nx,ny,nz,rhostr,tem1,tem2,tem3)
CALL vbcw
(nx,ny,nz,w,wcont,tbc,bbc,u,v, &
rhostr,tem1,tem2,tem3, &
j1,j2,j3)
!
!-----------------------------------------------------------------------
!
! Enforce boundary conditions
! Use zero gradient for lateral bc, and rigid for top,bottom.
!
!-----------------------------------------------------------------------
!
ebc = 3
wbc = 3
nbc = 3
sbc = 3
tbc = 1
bbc = 1
!
DO k=1,nz
DO j=1,ny
pdteb(j,k)=0.
pdtwb(j,k)=0.
END DO
DO i=1,nx
pdtnb(i,k)=0.
pdtsb(i,k)=0.
END DO
END DO
!
CALL bcu
(nx,ny,nz,0., u, &
pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc)
CALL bcv
(nx,ny,nz,0., v, &
pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc)
CALL lbcw
(nx,ny,nz,0.0,w,wcont, &
pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc)
CALL bcsclr
(nx,ny,nz,0.,ptprt,ptprt,ptprt, &
pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc)
CALL bcsclr
(nx,ny,nz,0., qv,qv,qv, &
pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc)
CALL bcsclr
(nx,ny,nz,0., qc,qc,qc, &
pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc)
CALL bcsclr
(nx,ny,nz,0., qr,qr,qr, &
pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc)
CALL bcsclr
(nx,ny,nz,0., qi,qi,qi, &
pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc)
CALL bcsclr
(nx,ny,nz,0., qs,qs,qs, &
pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc)
CALL bcsclr
(nx,ny,nz,0., qh,qh,qh, &
pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc)
!
!-----------------------------------------------------------------------
!
! Apply extrapolated gradient bc to pressure for geostrophic
! type balance to winds along lateral boundaries. Consistent
! with horizontal winds constant.
!
!-----------------------------------------------------------------------
!
DO k=2,nz-2
DO j=1,ny-1
pprt(1,j,k)=2.*pprt(2,j,k)-pprt(3,j,k)
END DO
END DO
DO k=2,nz-2
DO j=1,ny-1
pprt(nx-1,j,k)=2.*pprt(nx-2,j,k)-pprt(nx-3,j,k)
END DO
END DO
DO k=2,nz-2
DO i=1,nx-1
pprt(i,ny-1,k)=2.*pprt(i,ny-2,k)-pprt(i,ny-3,k)
END DO
END DO
DO k=2,nz-2
DO i=1,nx-1
pprt(i,1,k)=2.*pprt(i,2,k)-pprt(i,3,k)
END DO
END DO
!
!-----------------------------------------------------------------------
!
! For top and bottom bc, apply hydrostatic pressure equation to
! total pressure, then subtract pbar.
!
!-----------------------------------------------------------------------
!
pconst=-g/rd
DO j=1,ny-1
DO i=1,nx-1
pr2 = pbar(i,j,2)+pprt(i,j,2)
temp = (ptbar(i,j,2)+ptprt(i,j,2))* &
((pr2/p0)**rddcp)
tvbar = temp*(1.0+rvdrd*qv(i,j,2))/ &
(1.0+qv(i,j,2))
pr1=pr2*EXP(pconst*(zs(i,j,1)-zs(i,j,2))/tvbar)
pprt(i,j,1)=pr1-pbar(i,j,1)
pr2 = pbar(i,j,nz-2)+pprt(i,j,nz-2)
temp = (ptbar(i,j,nz-2)+ptprt(i,j,nz-2))* &
((pr2/p0)**rddcp)
tvbar= temp*(1.0+rvdrd*qv(i,j,nz-2))/ &
(1.0+qv(i,j,nz-2))
pr1=pr2*EXP(pconst*(zs(i,j,nz-1)-zs(i,j,nz-2))/tvbar)
pprt(i,j,nz-1)=pr1-pbar(i,j,nz-1)
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Adjust surface (skin) temperature based on radiation and
! current air temperature.
!
!-----------------------------------------------------------------------
!
IF(tsfcadj > 0 ) THEN
IF (radopt == 2) THEN
rbufsz = n2d_radiat*nx*ny + n3d_radiat*nx*ny*nz
ELSE
rbufsz = 1
END IF
ALLOCATE(rsirbm(nx,ny))
rsirbm = 0.
ALLOCATE(rsirdf(nx,ny))
rsirdf = 0.
ALLOCATE(rsuvbm(nx,ny))
rsuvbm = 0.
ALLOCATE(rsuvdf(nx,ny))
rsuvdf = 0.
ALLOCATE(cosz(nx,ny))
cosz = 0.
ALLOCATE(cosss(nx,ny))
cosss = 0.
ALLOCATE(fdirir(nx,ny))
fdirir = 0.
ALLOCATE(fdifir(nx,ny))
fdifir = 0.
ALLOCATE(fdirpar(nx,ny))
fdirpar= 0.
ALLOCATE(fdifpar(nx,ny))
fdifpar= 0.
ALLOCATE(st4(nx,ny))
st4 = 0.
ALLOCATE(tem11(nx,ny,nz))
tem11 = 0.
ALLOCATE(tem12(nx,ny,nz))
tem12 = 0.
ALLOCATE(tem13(nx,ny,nz))
tem13 = 0.
ALLOCATE(tem14(nx,ny,nz))
tem14 = 0.
ALLOCATE(tem15(nx,ny,nz))
tem15 = 0.
ALLOCATE(tem16(nx,ny,nz))
tem16 = 0.
ALLOCATE(radbuf(rbufsz))
radbuf = 0.
CALL adjtsfc
(nx,ny,nz,rbufsz, &
ptprt,pprt,qv,qc,qr,qi,qs,qh, &
ptbar,pbar,ppi,rhostr, &
x,y,z,zp,j3inv, &
soiltyp, tsfc, wetsfc,snowdpth, &
radfrc, radsw, rnflx, &
rsirbm,rsirdf,rsuvbm,rsuvdf, &
cosz, cosss, &
fdirir,fdifir,fdirpar,fdifpar,st4, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8, &
tem9,tem10,tem11,tem12,tem13,tem14,tem15,tem16, &
radbuf)
DEALLOCATE(rsirbm)
DEALLOCATE(rsirdf)
DEALLOCATE(rsuvbm)
DEALLOCATE(rsuvdf)
DEALLOCATE(cosz)
DEALLOCATE(cosss)
DEALLOCATE(fdirir)
DEALLOCATE(fdifir)
DEALLOCATE(fdirpar)
DEALLOCATE(fdifpar)
DEALLOCATE(st4)
DEALLOCATE(tem11)
DEALLOCATE(tem12)
DEALLOCATE(tem13)
DEALLOCATE(tem14)
DEALLOCATE(tem15)
DEALLOCATE(tem16)
DEALLOCATE(radbuf)
END IF
!
!-----------------------------------------------------------------------
!
! Calculate analysis increments and write to file, if desired.
!
!-----------------------------------------------------------------------
!
IF(incrdmp > 0) THEN
CALL incrcalc
(nx,ny,nz,nxndg,nyndg,nzndg, &
u,v,w,pprt,ptprt,qv, &
qc,qr,qi,qs,qh, &
uincr,vincr,wincr,pincr,ptincr,qvincr, &
qcincr,qrincr,qiincr,qsincr,qhincr, &
istatus)
CALL incrdump
(nxndg,nyndg,nzndg,incrdmp,incdmpf, &
uincr,vincr,wincr,pincr,ptincr,qvincr, &
qcincr,qrincr,qiincr,qsincr,qhincr, &
uincdmp,vincdmp,wincdmp, &
pincdmp,ptincdmp,qvincdmp, &
qcincdmp,qrincdmp,qiincdmp,qsincdmp,qhincdmp, &
istatus)
END IF
!
!-----------------------------------------------------------------------
!
! Data dump of the model grid and base state arrays:
!
!-----------------------------------------------------------------------
!
DO k=1,nz
DO j=1,ny
DO i=1,nx
tem1(i,j,k)=0.
END DO
END DO
END DO
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem2(i,j,k)=rhostr(i,j,k)/j3(i,j,k)
END DO
END DO
END DO
!
IF( hdmpfmt == 5 ) GO TO 700
IF( hdmpfmt == 9 ) GO TO 700
!
!-----------------------------------------------------------------------
!
! Find a unique name basdmpfn(1:lbasdmpf) for the grid and
! base state array dump file
!
!-----------------------------------------------------------------------
!
CALL gtbasfn
(runname(1:lfnkey),dirname,ldirnam,hdmpfmt, &
mgrid,nestgrd,basdmpfn,lbasdmpf)
WRITE(6,'(1x,a,a)') &
'Data dump of grid and base state arrays into file ', &
basdmpfn(1:lbasdmpf)
CALL setcornerll
( nx,ny,x,y ) ! set the corner lat/lon
grdbas = 1 ! Dump out grid and base state arrays only
CALL dtadump
(nx,ny,nz,nstyps, &
hdmpfmt,nchdmp,basdmpfn(1:lbasdmpf), &
grdbas,filcmprs, &
u,v,w,ptprt,pprt,qv, &
qc,qr,qi,qs,qh, &
tem1,tem1,tem1, &
ubar,vbar,tem1,ptbar,pbar,tem2,qvbar, &
x,y,z,zp,hterain,j1,j2,j3, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
tem3,tem4,tem5)
!
!-----------------------------------------------------------------------
!
! Find a unique name hdmpfn(1:ldmpf) for history dump data set
! at time 'curtim'.
!
!-----------------------------------------------------------------------
!
700 CONTINUE
CALL gtdmpfn
(runname(1:lfnkey),dirname, &
ldirnam,curtim,hdmpfmt, &
mgrid,nestgrd, hdmpfn, ldmpf)
WRITE(6,'(1x,a,a)') 'History data dump in file ',hdmpfn(1:ldmpf)
grdbas = 0 ! No base state or grid array is dumped.
CALL dtadump
(nx,ny,nz,nstyps, &
hdmpfmt,nchdmp,hdmpfn(1:ldmpf), &
grdbas,filcmprs, &
u,v,w,ptprt,pprt,qv, &
qc,qr,qi,qs,qh, &
tem1,tem1,tem1, &
ubar,vbar,tem1,ptbar,pbar,tem2,qvbar, &
x,y,z,zp,hterain, j1,j2,j3, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
tem3,tem4,tem5)
END PROGRAM arpsdas
!
SUBROUTINE exsufx(fname,lenfnm,suffix,maxsuf,dotloc,lensuf) 2
IMPLICIT NONE
INTEGER :: lenfnm
INTEGER :: maxsuf
CHARACTER (LEN=1) :: fname(lenfnm)
CHARACTER (LEN=1) :: suffix(maxsuf)
INTEGER :: lensuf
INTEGER :: dotloc
!
INTEGER :: i
!
DO i=lenfnm,1,-1
IF(fname(i) == '.') GO TO 200
END DO
dotloc=0
lensuf=0
DO i=1,maxsuf
suffix(i)=' '
END DO
RETURN
200 CONTINUE
dotloc=i
lensuf=MIN(maxsuf,lenfnm-i)
DO i=1,lensuf
suffix(i)=fname(dotloc+i)
END DO
RETURN
END SUBROUTINE exsufx
!
SUBROUTINE exfroot(fname,lenfnm,froot,mxroot,lenroot) 1
IMPLICIT NONE
INTEGER :: lenfnm
INTEGER :: mxroot
CHARACTER (LEN=1) :: fname(lenfnm)
CHARACTER (LEN=1) :: froot(mxroot)
INTEGER :: lenroot
!
INTEGER :: i
INTEGER :: slashloc
INTEGER :: dotloc
!
DO i=lenfnm,1,-1
IF(fname(i) == '/') EXIT
END DO
! 101 CONTINUE
slashloc=i
DO i=slashloc,lenfnm
IF(fname(i) == '.') EXIT
END DO
! 201 CONTINUE
dotloc=i
lenroot=(dotloc-slashloc)-1
DO i=1,lenroot
froot(i)=fname(slashloc+i)
END DO
RETURN
END SUBROUTINE exfroot
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SETCAT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE setcat(nx,ny,nz,ccatopt,zs, & 1,2
ptprt,pprt,qv,qc,qr,qi,qs,qh, &
ptbar,pbar, &
prcrate,icatg)
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Assign categories to grid locations according to the background
! state to account for decorrelation across areas with active
! parameterized convection and those without.
!
!
!-----------------------------------------------------------------------
!
!
! AUTHOR: Keith Brewster
! 8/7/98
!
! MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz
INTEGER :: ccatopt
REAL :: zs (nx,ny,nz) ! The physical height coordinate defined at
! w-point of the staggered grid.
REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature (K)
REAL :: pprt (nx,ny,nz) ! Perturbation pressure (Pascal)
REAL :: qv (nx,ny,nz) ! Water vapor specific humidity (kg/kg)
REAL :: qc (nx,ny,nz) ! Cloud water mixing ratio (kg/kg)
REAL :: qr (nx,ny,nz) ! Rain water mixing ratio (kg/kg)
REAL :: qi (nx,ny,nz) ! Cloud ice mixing ratio (kg/kg)
REAL :: qs (nx,ny,nz) ! Snow mixing ratio (kg/kg)
REAL :: qh (nx,ny,nz) ! Hail mixing ratio (kg/kg)
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal).
REAL :: prcrate(nx,ny,4) ! precipitation rate (kg/(m**2*s))
! prcrate(1,1,1) = total precipitation rate
! prcrate(1,1,2) = grid scale precip. rate
! prcrate(1,1,3) = cumulus precip. rate
! prcrate(1,1,4) = microphysics precip. rate
INTEGER :: icatg(nx,ny)
!
REAL :: prctrw,bltop,rhrkul,thevrkul
PARAMETER (prctrw=1.0E-04, & ! (kg/(m**2*s))
bltop=1000., & ! m
rhrkul=0.8, &
thevrkul=9.0) ! K**2
!
INTEGER :: i,j
! integer knt,k
INTEGER :: knt1,knt2a,knt2b,knt3,knt4
REAL :: pr,temp,qvsat,rh
! real rhsum,rhmean
! real the,thesum,thesq,themean,thevar
!
!-----------------------------------------------------------------------
!
! Function f_qvsat and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
REAL :: f_qvsat
!fpp$ expand (f_qvsat)
!dir$ inline always f_qvsat
!
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
!
INCLUDE 'phycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
knt1=0
knt2a=0
knt2b=0
knt3=0
knt4=0
!
!-----------------------------------------------------------------------
!
! Initialize with default of no precip or precip influence
!
!-----------------------------------------------------------------------
!
DO j=1,ny-1
DO i=1,nx-1
icatg(i,j)=1
END DO
END DO
!
IF(ccatopt == 1) THEN
DO j=1,ny-1
DO i=1,nx-1
!
!-----------------------------------------------------------------------
!
! Is convective precip turned on or heavy rain occuring?
!
!-----------------------------------------------------------------------
!
IF (prcrate(i,j,3) > 0. .OR. prcrate(i,j,1) > prctrw) THEN
knt4=knt4+1
icatg(i,j)=4
! print *, ' set i,j:',i,j,' cat 4',' prcrate(1):',
! : prcrate(i,j,1),' prcrate(3):',prcrate(i,j,3)
!
!-----------------------------------------------------------------------
!
! Is it currently raining?
!
!-----------------------------------------------------------------------
!
ELSE IF ( prcrate(i,j,1) > 0. ) THEN
knt3=knt3+1
icatg(i,j)=3
! print *, ' set i,j:',i,j,' cat 3',' prcrate(1):',
! : prcrate(i,j,1),' prcrate(3):',prcrate(i,j,3)
!
!-----------------------------------------------------------------------
!
! Evidence of rain-cooled air?
! Lapse rate close to a moist adiabatic lapse rate in the lowest 2-km,
! and a high relative humidity in that layer. Or a high relative
! humidity at the first level above ground.
!
!-----------------------------------------------------------------------
!
ELSE
!
!-----------------------------------------------------------------------
!
! For the layer 0-2km above first grid level (k=2).
! Calculate the mean relative humidity.
! Calculate the standard deviation of saturated equivalent
! potential temperature.
!
!-----------------------------------------------------------------------
!
! rhsum=0.
! thesum=0.
! thesq=0.
! knt=0
! DO 40 k=2,nz-1
! IF(k.gt.2 .and. (zs(i,j,k)-zs(i,j,2)).gt.bltop) GO TO 41
! knt=knt+1
! pr=pprt(i,j,k)+pbar(i,j,k)
! temp=(ptprt(i,j,k)+ptbar(i,j,k)) * (pr/p0)**rddcp
! qvsat=f_qvsat( pr,temp )
! rh=min(1.,(qv(i,j,k)/qvsat))
! the=F_PT2PTE( pr,(ptprt(i,j,k)+ptbar(i,j,k)),qvsat)
! rhsum=rhsum+rh
! thesum=thesum+the
! thesq=thesq+(the*the)
! 40 CONTINUE
! 41 CONTINUE
! IF(knt.gt.0) THEN
! rhmean=rhsum/float(knt)
! themean=thesum/float(knt)
! ELSE
! rhmean=0.
! themean=0.
! END IF
! IF(knt.gt.1) THEN
! thevar=(thesq-((thesum*thesum)/float(knt)))/float(knt-1)
! ELSE
! thevar=9999.
! END IF
!
! IF (rhmean.gt.rhrkul .and. thevar.lt.thevrkul) THEN
! knt2a=knt2a+1
! icatg(i,j)=2
! print *, ' set i,j:',i,j,' cat 2',' prcrate(1):',
! : prcrate(i,j,1),' prcrate(3):',prcrate(i,j,3)
! ELSE
!
pr=pprt(i,j,2)+pbar(i,j,2)
temp=(ptprt(i,j,2)+ptbar(i,j,2)) * (pr/p0)**rddcp
qvsat=f_qvsat
( pr,temp )
rh=qv(i,j,2)/qvsat
IF(rh > 0.8) THEN
knt2b=knt2b+1
icatg(i,j)=2
! print *, ' set i,j:',i,j,' cat 2*',' RH(2) = ',rh
END IF
! END IF
! IF (mod(i,20).eq.0 .and. mod(j,10).eq.0) THEN
! print *, i,j,' rhmean=',rhmean,
! : ' thevar=',thevar,' icat=',icatg(i,j)
! print *, ' '
! END IF
END IF
END DO
END DO
knt1=((nx-1)*(ny-1))-knt2a-knt2b-knt3-knt4
WRITE(6,'(a)') ' Precip correlation category assignment counts'
WRITE(6,'(a,i6)') ' cat 1: ',knt1
WRITE(6,'(a,i6)') ' cat 2a: ',knt2a
WRITE(6,'(a,i6)') ' cat 2b: ',knt2b
WRITE(6,'(a,i6)') ' cat 3: ',knt3
WRITE(6,'(a,i6)') ' cat 4: ',knt4
END IF
!
CALL iedgfill
(icatg,1,nx,1,nx-1, 1,ny,1,ny-1,1,1,1,1)
!
RETURN
END SUBROUTINE setcat