PROGRAM arpsdas,315
!
!##################################################################
!##################################################################
!###### ######
!###### 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.
!
! 06/10/2002 (Eric Kemp)
! Addition of new soil model variables.
!
! 04/09/2003 (K. Brewster)
! Modified I/O of source errors to stop if I/O problem.
! Correction of problem in super-obbing.
!
! 03/01/2005 (S. Leyton and D. Weber)
! Add MPI capabilities.
!
! 12/01/2005 (K. W. Thomas)
! Update MPI capabilities.
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INCLUDE 'alloc.inc'
INCLUDE 'adas.inc'
INCLUDE 'nudging.inc'
INCLUDE 'exbc.inc'
INCLUDE 'radcst.inc'
INCLUDE 'mp.inc' ! Message passing parameters
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 :: zpsoil(:,:,:) ! Soil level depth.
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 :: j3soil(:,:,:) ! Coordinate transformation Jacobian defined
! as d( zpsoil )/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
REAL, ALLOCATABLE :: j3soilinv (:,:,:)! Inverse of j3soil
!
!-----------------------------------------------------------------------
!
! 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
REAL, ALLOCATABLE :: radswnet(:,:) ! Net shortwave radiation
REAL, ALLOCATABLE :: radlwin(:,:) ! Incominging longwave radiation
!
!-----------------------------------------------------------------------
!
! 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 :: tsoil(:,:,:,:) ! Soil temperature (K)
REAL, ALLOCATABLE :: qsoil(:,:,:,:) ! Soil moisture
REAL, ALLOCATABLE :: qvsfc(:,:,:) ! Effective qv at sfc.
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(:), xslg(:)
REAL, ALLOCATABLE :: ys(:), yslg(:)
REAL, ALLOCATABLE :: zs(:,:,:)
REAL, ALLOCATABLE :: latgr(:,:)
REAL, ALLOCATABLE :: longr(:,:)
!
!-----------------------------------------------------------------------
!
! Temporary arrays. Some are only used for MPI
!
!-----------------------------------------------------------------------
!
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.
REAL, ALLOCATABLE :: tem11 (:,:,:) ! Temporary work array.
INTEGER, ALLOCATABLE :: item1(:) ! Temporary work array.
INTEGER, ALLOCATABLE :: ibegin(:)
INTEGER, ALLOCATABLE :: iend(:)
INTEGER, ALLOCATABLE :: tems1di(:)
INTEGER, ALLOCATABLE :: temr1di(:)
REAL, ALLOCATABLE :: tems1dr(:)
REAL, ALLOCATABLE :: temr1dr(:)
REAL, ALLOCATABLE :: tems2dr(:,:)
REAL, ALLOCATABLE :: temr2dr(:,:)
!
! Multi-level needs to handle "trnua", which is 2D. The above two temporary
! arrays need to be preserved for "anxiter".
!
REAL, ALLOCATABLE :: tems2drua(:,:)
REAL, ALLOCATABLE :: temr2drua(:,:)
REAL, ALLOCATABLE :: tems3dr(:,:,:)
REAL, ALLOCATABLE :: temr3dr(:,:,:)
!
!-----------------------------------------------------------------------
!
! 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 :: radbuf(:)
REAL, ALLOCATABLE :: sh(:,:) !added by DTD
REAL, ALLOCATABLE :: temrad1 (:,:,:) ! 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.
REAL, ALLOCATABLE :: tem17 (:,:,:) ! Temporary work array. added by DTD
REAL, ALLOCATABLE :: tem1soil (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem2soil (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem3soil (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem4soil (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem5soil (:,:,:) ! 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, 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 :: nzsoil ! Number of soil levels
INTEGER :: nxlg ! Number of grid points in the x-direction large grid
INTEGER :: nylg ! Number of grid points in the y-direction large grid
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
REAL :: rlim ! Largest possible radius of influence
!-----------------------------------------------------------------------
!
! Analysis parameters
!
!-----------------------------------------------------------------------
!
REAL, PARAMETER :: rhmin=0.05 ! rh safety net value to prevent neg qv
INTEGER, PARAMETER :: iqspr=3 ! Use qobs of pstn to combine x,y,elev
INTEGER, PARAMETER :: klim= 1 ! Min of one other sfc station for QC
!
!-----------------------------------------------------------------------
!
! Indices of specific observation variables
!
!-----------------------------------------------------------------------
!
INTEGER, PARAMETER :: ipres=3 ! Pressure
INTEGER, PARAMETER :: iptmp=4 ! Potential Temperature
INTEGER, PARAMETER :: iqv=5 ! Specific Humidity
!
!-----------------------------------------------------------------------
!
! 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
!
!-----------------------------------------------------------------------
!
! 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)
REAL :: lonsng(mx_sng,ntime)
REAL :: hgtsng(mx_sng,ntime)
REAL :: xsng(mx_sng)
REAL :: ysng(mx_sng)
REAL :: trnsng(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)
REAL :: yua(mx_ua)
REAL :: trnua(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)
REAL :: lonradc(mx_colrad)
REAL :: xradc(mx_colrad)
REAL :: yradc(mx_colrad)
REAL :: trnradc(mx_colrad)
REAL :: distrad(mx_colrad)
REAL :: uazmrad(mx_colrad)
REAL :: 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)
REAL :: lonret(mx_ret)
REAL :: elvret(mx_ret)
!
!-----------------------------------------------------------------------
!
! Retrieval observation variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iret(mx_colret)
INTEGER :: nlevret(mx_colret)
REAL :: latretc(mx_colret)
REAL :: lonretc(mx_colret)
REAL :: xretc(mx_colret)
REAL :: yretc(mx_colret)
REAL :: trnretc(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)
REAL :: barqclim(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=256) :: froot
CHARACTER (LEN=256) :: qclist
CHARACTER (LEN=256) :: qcsave
CHARACTER (LEN=256) :: stats ! runname + '.lst'
INTEGER :: istatus,jstatus,istat
INTEGER :: nobsng,nobsrd(ntime)
INTEGER :: i4timanx
INTEGER :: lenfnm,maxsuf,lensuf,dotloc,mxroot,lenroot
INTEGER :: iqclist,iqcsave,iwstat
CHARACTER (LEN=256) :: status_file ! runname + '.adasstat'
INTEGER n_used_sng(nsrc_sng), n_used_ua(nsrc_sng)
INTEGER jsta, klev, ngoodlev
!
!-----------------------------------------------------------------------
!
! MPI bookkeeping
!
!-----------------------------------------------------------------------
!
INTEGER :: indexsng(mx_sng) ! indexes each data point according to
! which local proc domain it resides
INTEGER, ALLOCATABLE :: ksng(:) ! Number of obs owned by each processor
INTEGER :: ksngmax ! Max "ksng" value.
INTEGER :: indexua(mx_ua) ! indexes each data point according to
! which local proc domain it resides
INTEGER, ALLOCATABLE :: kua(:) ! Number of obs owned by each processor
INTEGER :: kuamax ! Max "kua" value.
INTEGER :: indexrad(mx_rad) ! indexes each data point according to
! which local proc domain it resides
INTEGER :: indexret(mx_ret) ! indexes each data point according to
! which local proc domain it resides
INTEGER :: iproc ! Radii of influence "i" direction
INTEGER :: jproc ! Radii of influence "j" direction
INTEGER, ALLOCATABLE :: mpi_map(:,:) ! Map of communication entries
INTEGER :: nmap ! Number of entries in "mpi_map".
!
!-----------------------------------------------------------------------
!
! 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=256) :: basdmpfn
CHARACTER (LEN=80 ) :: header
INTEGER :: i,j,k,ios,ktab,isrc,jsrc,ivar,ifile
INTEGER :: grdbas,lbasdmpf,rbufsz
INTEGER :: nobsua,ncolrad,ncolret,nmeso,nsao
REAL :: tgrid,qvmin,qvsat,qsrcmax,rngsq
REAL :: temp,tvbar,qvprt,qtot,pconst,pr1,pr2
!
INTEGER :: ixrad(mx_rad),jyrad(mx_rad)
INTEGER :: retqropt
PARAMETER (retqropt=1)
INTEGER :: imstat
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! 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,nzsoil,nstyps)
IF (mp_opt > 0) THEN
nxlg = (nx - 3) * nproc_x + 3
nylg = (ny - 3) * nproc_y + 3
ELSE
nxlg = nx
nylg = ny
END IF
IF (lbcopt == 2) THEN
IF (myproc == 0) THEN
WRITE (*,*) "INITADAS: resetting lbcopt to 1 & lateral bc's to 4"
END IF
lbcopt = 1
ebc = 4
wbc = 4
nbc = 4
sbc = 4
END IF
!-----------------------------------------------------------------------
! Read in adas namelist parameters
!-----------------------------------------------------------------------
!
CALL initadas
IF (mp_opt > 0) THEN
CALL radinf(rlim)
!
! See what our processor radius of influence is in each direction.
!
iproc = INT( rlim / ( (nx-3) * dx )) + 1
jproc = INT( rlim / ( (ny-3) * dy )) + 1
IF (myproc == 0) THEN
WRITE(6,*) 'Process radius of influence x & y ',iproc,jproc
END IF
ELSE
!
! We'll allocate a trivial array.
!
iproc = 1
jproc = 1
ENDIF
IF(incrdmp > 0 ) THEN
nxndg=nx
nyndg=ny
nzndg=nz
ELSE
nxndg=1
nyndg=1
nzndg=1
END IF
!
!-----------------------------------------------------------------------
!
! Allocate adas arrays
!
!-----------------------------------------------------------------------
!
ALLOCATE(x(nx),STAT=istatus)
CALL check_alloc_status(istatus, "adas:x")
ALLOCATE(y(ny),STAT=istatus)
CALL check_alloc_status(istatus, "adas:y")
ALLOCATE(z(nz),STAT=istatus)
CALL check_alloc_status(istatus, "adas:z")
CALL check_alloc_status(istatus, "adas:y")
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(zpsoil(nx,ny,nzsoil),STAT=istatus)
CALL check_alloc_status(istatus, "adas:zpsoil")
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(j3soil(nx,ny,nzsoil),STAT=istatus)
CALL check_alloc_status(istatus, "adas:j3soil")
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(j3soilinv(nx,ny,nzsoil),STAT=istatus)
CALL check_alloc_status(istatus, "adas:j3soilinv")
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(radswnet (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "adas:radswnet")
ALLOCATE(radlwin (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "adas:radlwin")
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(qvsfc (nx,ny,0:nstyps),STAT=istatus)
CALL check_alloc_status(istatus, "adas:qvsfc")
ALLOCATE(tsoil (nx,ny,nzsoil,0:nstyps),STAT=istatus)
CALL check_alloc_status(istatus, "adas:tsoil")
ALLOCATE(qsoil (nx,ny,nzsoil,0:nstyps),STAT=istatus)
CALL check_alloc_status(istatus, "adas:qsoil")
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(xslg(nxlg),STAT=istatus)
CALL check_alloc_status(istatus, "adas:xslg")
ALLOCATE(yslg(nylg),STAT=istatus)
CALL check_alloc_status(istatus, "adas:yslg")
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(tem1soil(nx,ny,nzsoil),STAT=istatus)
CALL check_alloc_status(istatus, "adas:tem1soil")
ALLOCATE(tem2soil(nx,ny,nzsoil),STAT=istatus)
CALL check_alloc_status(istatus, "adas:tem2soil")
ALLOCATE(tem3soil(nx,ny,nzsoil),STAT=istatus)
CALL check_alloc_status(istatus, "adas:tem3soil")
ALLOCATE(tem4soil(nx,ny,nzsoil),STAT=istatus)
CALL check_alloc_status(istatus, "adas:tem4soil")
ALLOCATE(tem5soil(nx,ny,nzsoil),STAT=istatus)
CALL check_alloc_status(istatus, "adas:tem5soil")
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(tem11(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "adas:tem11")
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")
nmap = ( 2 * iproc + 1 ) * ( 2 * jproc + 1 )
ALLOCATE(mpi_map(nmap,2),STAT=istatus)
CALL check_alloc_status(istatus, "adas:mpi:map")
IF (myproc == 0) 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
zpsoil =0.0
j1 =0.0
j2 =0.0
j3 =0.0
j3soil =0.0
aj3x=0.0
aj3y=0.0
aj3z=0.0
j3inv=0.0
j3soilinv=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
radswnet =0.0
radlwin =0.0
soiltyp=0.0
stypfrct=0.0
vegtyp =0.0
lai =0.0
roufns =0.0
veg =0.0
qvsfc =0.0
tsoil =0.0
qsoil =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
xslg = 0.0
yslg = 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
qsrcsng = 0.0
ibegin=0
iend=0
anx=0.0
xcor=0.0
icatg=0.0
!
!-----------------------------------------------------------------------
!
! Read in quality control info. Most of the information will not be
! broadcast, as processor 0 will handle pre-processing of obs before
! the data gets sent to the other processors.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Set expected _squared_ background error for each variable
! This will depend on the particular background field used,
! season, age of forecast, etc.
!
!-----------------------------------------------------------------------
!
IF (myproc == 0) THEN
WRITE(6,'(/,a,a)') ' 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,*,iostat=istat) hqback(ktab), &
qback(3,ktab), &
qback(4,ktab), &
qback(5,ktab), &
qback(1,ktab), &
qback(2,ktab)
IF( istat /= 0 ) EXIT
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
nlvqback=ktab-1
CLOSE(4)
DO jsrc=1,nsrc_sng
IF(srcsng(jsrc) /= 'NULL') THEN
WRITE(6,'(/,a,a)') ' Reading ', TRIM(sngerrfil(jsrc))
OPEN(4,FILE=trim(sngerrfil(jsrc)),STATUS='old')
READ(4,'(a8)') 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)
CALL arpsstop("mismatch",1)
END IF
READ(4,'(a80)') srcsng_full(jsrc)
READ(4,*) qcmulsng(jsrc)
READ(4,'(a80)') header
READ(4,*) qsrcsng(3,jsrc), &
qsrcsng(4,jsrc), &
qsrcsng(5,jsrc), &
qsrcsng(1,jsrc), &
qsrcsng(2,jsrc)
WRITE(6,'(//,a,a,/a)') 'Single-level std error for ', &
srcsng(jsrc),srcsng_full(jsrc)
WRITE(6,'(a,f5.1)') ' QC multiplier: ',qcmulsng(jsrc)
WRITE(6,'(1x,a)') &
' u (m/s) v (m/s) pres(mb) temp(K) RH(%)'
WRITE(6,'(1x,5f8.2)') (qsrcsng(ivar,jsrc),ivar=1,5)
qsrcsng(3,jsrc)=100.*qsrcsng(3,jsrc)
qsrcsng(5,jsrc)=0.01*qsrcsng(5,jsrc)
CLOSE(4)
ELSE
DO k=1,npass
iusesng(jsrc,k) = 0
END DO
END IF
END DO
DO jsrc=1,nsrc_ua
IF(srcua(jsrc) /= 'NULL') THEN
WRITE(6,'(/,a,a)') ' Reading ', TRIM(uaerrfil(jsrc))
OPEN(4,FILE=trim(uaerrfil(jsrc)),STATUS='old')
READ(4,'(a8)') 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)
CALL arpsstop("mismatch",1)
END IF
READ(4,'(a80)') srcua_full(jsrc)
READ(4,*) qcmulua(jsrc)
READ(4,'(a80)') header
WRITE(6,'(//,a,a,/a)') 'UA data std error for ', &
srcua(jsrc),srcua_full(jsrc)
WRITE(6,'(a,f5.1)') ' QC multiplier: ',qcmulua(jsrc)
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,*,iostat=istat) huaqsrc(ktab,jsrc), &
qsrcua(3,ktab,jsrc), &
qsrcua(4,ktab,jsrc), &
qsrcua(5,ktab,jsrc), &
qsrcua(1,ktab,jsrc), &
qsrcua(2,ktab,jsrc)
IF( istat /= 0 ) EXIT
WRITE(6,'(1x,i4,f10.2,5f8.2)') ktab,huaqsrc(ktab,jsrc), &
(qsrcua(ivar,ktab,jsrc),ivar=1,5)
qsrcua(3,ktab,jsrc)=100.*qsrcua(3,ktab,jsrc)
qsrcua(5,ktab,jsrc)=0.01*qsrcua(5,ktab,jsrc)
END DO
nlvuatab(jsrc)=ktab-1
CLOSE(4)
ELSE
DO k=1,npass
iuseua(jsrc,k) = 0
END DO
END IF
END DO
DO jsrc=1,nsrc_rad
IF(srcrad(jsrc) /= 'NULL') THEN
WRITE(6,'(/,a,a)') ' Reading ', TRIM(raderrfil(jsrc))
OPEN(4,FILE=trim(raderrfil(jsrc)),STATUS='old')
READ(4,'(a8)') 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)
CALL arpsstop("mismatch",1)
END IF
READ(4,'(a80)') srcrad_full(jsrc)
READ(4,*) qcmulrad(jsrc)
READ(4,'(a80)') header
READ(4,*) qsrcrad(1,jsrc), &
qsrcrad(2,jsrc), &
qsrcrad(3,jsrc)
WRITE(6,'(//,a,a,/a)') 'Radar data std error for ', &
srcrad(jsrc),srcrad_full(jsrc)
WRITE(6,'(a,f5.1)') ' QC multiplier: ',qcmulrad(jsrc)
WRITE(6,'(1x,a)') &
' ref(dBz) Vrad(m/s) SpWid(m/s)'
WRITE(6,'(1x,4f8.2)') &
(qsrcrad(ivar,jsrc),ivar=1,nvar_radin)
CLOSE(4)
qsrcsng(5,jsrc)=0.01*qsrcsng(5,jsrc)
CLOSE(4)
ELSE
DO k=1,npass
iuserad(jsrc,k) = 0
END DO
END IF
END DO
DO jsrc=1,nsrc_ret
IF(srcret(jsrc) /= 'NULL') THEN
WRITE(6,'(/,a,a)') ' Reading ', TRIM(reterrfil(jsrc))
OPEN(4,FILE=trim(reterrfil(jsrc)),STATUS='old')
READ(4,'(a8)') 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)
CALL arpsstop("mismatch",1)
END IF
READ(4,'(a80)') srcret_full(jsrc)
READ(4,*) qcmulret(jsrc)
READ(4,'(a80)') header
WRITE(6,'(//,a,a,/a)') 'Retrieval std error for ', &
srcret(jsrc),srcret_full(jsrc)
WRITE(6,'(a,f5.1)') ' QC multiplier: ',qcmulrad(jsrc)
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,*,iostat=istat) hrtqsrc(ktab,jsrc), &
qsrcret(3,ktab,jsrc), &
qsrcret(4,ktab,jsrc), &
qsrcret(5,ktab,jsrc), &
qsrcret(1,ktab,jsrc), &
qsrcret(2,ktab,jsrc)
IF( istat /= 0 ) EXIT
WRITE(6,'(1x,i4,f10.2,5f8.2)') ktab,hrtqsrc(ktab,jsrc), &
(qsrcret(ivar,ktab,jsrc),ivar=1,5)
qsrcret(3,ktab,jsrc)=100.*qsrcret(3,ktab,jsrc)
qsrcret(5,ktab,jsrc)=0.01*qsrcret(5,ktab,jsrc)
END DO
nlvrttab(jsrc)=ktab-1
CLOSE(4)
ELSE
DO k=1,npass
iuseret(jsrc,k) = 0
END DO
END IF
END DO
!
!-----------------------------------------------------------------------
!
! 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))
barqclim(ivar,isrc)=qcmulsng(isrc)*SQRT(qsrcsng(ivar,isrc))
END DO
qcthrsng(iqv,isrc)=min(qcthrsng(iqv,isrc),0.75)
barqclim(iqv,isrc)=min(barqclim(iqv,isrc),0.75)
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
qcthrua(iqv,isrc)=min(qcthrua(iqv,isrc),0.75)
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
qcthrret(iqv,isrc)=min(qcthrret(iqv,isrc),0.75)
END DO
END IF
!
! Broadcast read variables to all procs. Some are dead, so we keep them
! that way. Others are only used by processor zero initially.
!
! Also broadcast variables that may have been updated.
!
CALL mpupdater(hqback, nz_tab)
CALL mpupdater(qback, (nvar_anx*nz_tab))
CALL mpupdatei(nlvqback, 1)
CALL mpupdatec(srcsng, (8*nsrc_sng))
CALL mpupdater(qcthrsng, (nvar_anx*nsrc_sng))
CALL mpupdater(qcthrua, (nvar_anx*nsrc_ua))
CALL mpupdater(qcthrrad, (nvar_radin*nsrc_rad))
CALL mpupdater(qcthrret, (nvar_anx*nsrc_ret))
CALL mpupdatei(iusesng,(nsrc_sng+1)*mx_pass)
CALL mpupdatei(iuseua,(nsrc_ua+1)*mx_pass)
CALL mpupdatei(iuserad,(nsrc_rad+1)*mx_pass)
CALL mpupdatei(iuseret,(nsrc_ret+1)*mx_pass)
CALL make_mpi_map(mpi_map,nmap,iproc,jproc,nx,ny)
!
!-----------------------------------------------------------------------
!
! 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,nzsoil,1,nstyps,1, &
x,y,z,zp,zpsoil,hterain,mapfct, &
j1,j2,j3,j3soil,aj3x,aj3y,aj3z,j3inv,j3soilinv, &
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, &
tsoil,qsoil,wetcanp,snowdpth,qvsfc, &
ptcumsrc,qcumsrc,w0avg,nca,kfraincv, &
cldefi,xland,bmjraincv, &
raing,rainc,prcrate,exbcbuf, &
radfrc,radsw,rnflx,radswnet,radlwin, &
usflx,vsflx,ptsflx,qvsflx, &
tem1soil,tem2soil,tem3soil,tem4soil,tem5soil, &
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
!
!-----------------------------------------------------------------------
!
! In the MPI world, some subroutines will need to know characteristics
! of the large grid. To simply later subroutine calls (one call not two
! surrounded by IF/ELSE/ENDIF, everybody that always needs the large
! grid will use the large grid variables. Large grid info is used by
! early decision making when determining what obs are needed.
!
!-----------------------------------------------------------------------
!
IF (mp_opt > 0 ) THEN
CALL mpimerge1dx(xs,nx,xslg)
CALL mpimerge1dy(ys,ny,yslg)
CALL mpupdater(xslg,nxlg)
CALL mpupdater(yslg,nylg)
ELSE
xslg = xs
yslg = ys
END IF
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
!
! Look at how "u" and "v" are defined as to whether we need to pass data
! to obtain the correct average of "u" and "v" at nx-1 and ny-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
IF (mp_opt > 0) THEN
DO k=1,nvar_anx
IF (k == 1) THEN ! u
j = 1
ELSE IF (k == 2) THEN ! v
j = 2
ELSE ! not u or v
j = 0
END IF
CALL mpsendrecv2dew(anx(1,1,1,k),nx,ny,nz,ebc,wbc,j,tem1)
CALL mpsendrecv2dns(anx(1,1,1,k),nx,ny,nz,nbc,sbc,j,tem1)
END DO
END IF
!-----------------------------------------------------------------------
!
! Save background fields on staggered grid in the increment
! arrays.
!
! PROBABLY NEEDS MPI WORK!
!
!-----------------------------------------------------------------------
!
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
!
!-----------------------------------------------------------------------
!
iwstat = 0
IF (myproc == 0) THEN
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
END IF
IF(ios /= 0) iwstat=0
maxsuf=LEN(suffix)
mxroot=LEN(froot)
nobsng=0
iqclist = 0
iqcsave = 0
!
! Read and preprocess the single level data, all done on processor 0.
! Broadcast the data to the rest of the ADAS world.
!
IF (myproc == 0) THEN
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.
!
!-----------------------------------------------------------------------
!
rngsq=sfcqcrng*sfcqcrng
IF (myproc == 0) PRINT *, 'Calling prepsfcobs'
CALL prepsfcobs(ntime,mx_sng, &
nvar_sng,nvar_anx,nsrc_sng,nobsng,ipres,iptmp,iqv, &
sngfname(ifile),sngtmchk(ifile),blackfil, &
var_sfc,var_anx,srcsng,qsrcsng, &
rmiss,iqspr,sprdist,nobsrd,timesng, &
stnsng,csrcsng,isrcsng,latsng,lonsng,hgtsng,xsng,ysng, &
nxlg,nylg,xslg,yslg, &
wx,kloud,idp3,store_emv,store_hgt,store_amt, &
obrdsng,obsng,qualrdsng,qobsng,qualsng, &
ival,climin,climax, &
rngsq,klim,wlim,qcthrsng,barqclim, &
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. Note that we are
! using the global domain here, even if we're in MPI mode, as only
! processor 0 is doing this.
!
!-----------------------------------------------------------------------
!
IF (myproc == 0) 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,nxlg,nylg,nz,nvar_anx,anx,xslg,yslg,zp, &
tem1,tem2,tem3,tem4,tem5,tem6, &
rmiss,nobsng,istatus)
END IF
END DO
END IF ! myproc == 0
!
!-----------------------------------------------------------------------
!
! MPI issues. Broadcast the number of data points, then the data. If
! there are no data values, we don't waste anyone's time.
!
!-----------------------------------------------------------------------
!
IF (mp_opt > 0) THEN
CALL mpupdatei(nobsng,1)
IF (nobsng > 0) THEN
CALL mpupdatei(timesng,mx_sng*ntime)
CALL mpupdatec(stnsng,mx_sng*ntime*5)
CALL mpupdatec(csrcsng,mx_sng*ntime*8)
CALL mpupdatei(isrcsng,mx_sng)
CALL mpupdater(latsng,mx_sng*ntime)
CALL mpupdater(lonsng,mx_sng*ntime)
CALL mpupdater(hgtsng,mx_sng*ntime)
CALL mpupdater(xsng,mx_sng)
CALL mpupdater(ysng,mx_sng)
CALL mpupdatec(wx,mx_sng*ntime*8)
CALL mpupdatei(kloud,mx_sng*ntime)
CALL mpupdatei(idp3,mx_sng*ntime)
CALL mpupdatec(store_emv,mx_sng*5*ntime)
CALL mpupdater(store_hgt,mx_sng*5*ntime)
CALL mpupdatec(store_amt,mx_sng*5*ntime*4)
CALL mpupdater(obrdsng,mx_sng*nvar_sng*ntime)
CALL mpupdater(obsng,nvar_anx*mx_sng)
CALL mpupdater(qobsng,nvar_anx*mx_sng)
CALL mpupdatei(qualsng,nvar_anx*mx_sng)
CALL mpupdater(knt,nvar_anx*nz)
CALL mpupdater(wgtsum,nvar_anx*nz)
CALL mpupdater(zsum,nvar_anx*nz)
!
!-----------------------------------------------------------------------
! We need to know which processor "owns" which obs, so they can be
! consulted on an "as needed" basis.
!-----------------------------------------------------------------------
!
ALLOCATE(item1(nobsng),STAT=istatus)
CALL check_alloc_status(istatus, "adas:item1:nobsng")
ALLOCATE(ksng(nprocs),STAT=istatus)
CALL check_alloc_status(istatus, "adas:ksng:nprocs")
CALL mpiprocess(nobsng,indexsng,nprocs,ksng,ksngmax, &
isrcsng,item1,nx,ny,xsng,ysng,xs,ys)
DEALLOCATE (item1)
!
!-----------------------------------------------------------------------
! Mark obs that we don't "own" so we don't do unnecessary
! computations.
!-----------------------------------------------------------------------
!
ALLOCATE(item1(0:nprocs-1),STAT=istatus)
CALL check_alloc_status(istatus, "adas:item1:nprocs")
item1 = 0
DO k=1,nmap
IF (mpi_map(k,1) .NE. -1 ) item1(mpi_map(k,1)) = 1
END DO
item1(myproc) = 1
DO k=1,nobsng
IF (indexsng(k) == -1) CYCLE
IF (isrcsng(k) == -1) CYCLE
IF (item1(indexsng(k)) == 0) isrcsng(k) = 0
END DO
DEALLOCATE (item1)
!
!-----------------------------------------------------------------------
!
! Although we're similar to PREPSFCGLOBS, and need to do one final task
! (compute heights), we now use the local grid variables, not the
! global.
!
!-----------------------------------------------------------------------
!
CALL prepsglmdcrs_sm(mx_sng,ntime,latsng,lonsng, &
xsng,ysng,hgtsng,csrcsng, &
nx,ny,nz,nvar_anx,anx,xs,ys,zp, &
tem1,tem2,tem3,tem4,tem5,tem6, &
rmiss,nobsng,indexsng,istatus)
END IF
!
! Collect and distribute the data.
!
IF (mp_opt > 0 .AND. nobsng > 0) THEN
ALLOCATE(tems1dr(ksngmax),STAT=istat)
CALL check_alloc_status(istat, "mpi_1dr_collect:tmps")
ALLOCATE(temr1dr(ksngmax),STAT=istat)
CALL check_alloc_status(istat, "mpi_1dr_collect:tmpr")
CALL mpi_1dr_collect(hgtsng,nobsng,indexsng, &
nprocs, ksng, ksngmax, mpi_map, nmap, tems1dr, temr1dr)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Calculate initial obs differences for single-level data
!
!-----------------------------------------------------------------------
!
IF (myproc == 0 ) PRINT *, 'Calling grdtosng'
CALL grdtosng(nx,ny,nxlg,nylg,nz,nz_tab,mx_sng,nvar_anx,nobsng, &
xs,ys,xslg,yslg,zp,icatg,anx,qback,hqback,nlvqback, &
tem1,tem2,tem3,tem4,tem5,tem6, &
stnsng,isrcsng,icatsng,hgtsng,xsng,ysng, &
obsng,qobsng,qualsng, &
odifsng,oanxsng,thesng,trnsng,indexsng)
!
!-----------------------------------------------------------------------
! Some data needs to be collected and distributed.
!-----------------------------------------------------------------------
!
IF ( mp_opt > 0 .AND. nobsng > 0 ) THEN
ALLOCATE(tems2dr(nvar_anx,ksngmax),STAT=istat)
CALL check_alloc_status(istat, "mpi_2dr_collect:tmps")
ALLOCATE(temr2dr(nvar_anx,ksngmax),STAT=istat)
CALL check_alloc_status(istat, "mpi_2dr_collect:tmpr")
!
! 1D real temporary arrays already allocated, but 1D ints aren't.
!
ALLOCATE(tems1di(ksngmax),STAT=istat)
CALL check_alloc_status(istat, "mpi_1di_collect:tmps")
ALLOCATE(temr1di(ksngmax),STAT=istat)
CALL check_alloc_status(istat, "mpi_1di_collect:tmpr")
CALL mpi_2dr_collect( qobsng, nvar_anx, mx_sng, nobsng, indexsng, &
nprocs, ksng, ksngmax, mpi_map, nmap, tems2dr, temr2dr)
CALL mpi_2dr_collect( odifsng, nvar_anx, mx_sng, nobsng, indexsng, &
nprocs, ksng, ksngmax, mpi_map, nmap, tems2dr, temr2dr)
CALL mpi_2dr_collect( oanxsng, nvar_anx, mx_sng, nobsng, indexsng, &
nprocs, ksng, ksngmax, mpi_map, nmap, tems2dr, temr2dr)
CALL mpi_1di_collect( icatsng, nobsng, indexsng, &
nprocs, ksng, ksngmax, mpi_map, nmap, tems1di, temr1di)
CALL mpi_1dr_collect( thesng, nobsng, indexsng, &
nprocs, ksng, ksngmax, mpi_map, nmap, tems1dr, temr1dr)
CALL mpi_1dr_collect( trnsng, nobsng, indexsng, &
nprocs, ksng, ksngmax, mpi_map, nmap, tems1dr, temr1dr)
!
! 2D space can't be deallocated, as the arrays get passed to "anxiter".
!
DEALLOCATE(tems1di)
DEALLOCATE(temr1di)
DEALLOCATE(tems1dr)
DEALLOCATE(temr1dr)
END IF
!
!-----------------------------------------------------------------------
!
! Read upper-air data, QC and convert units
!
!-----------------------------------------------------------------------
!
IF (myproc == 0) THEN
PRINT *, 'Calling prepuaobs'
CALL prepuaobs(nx,ny,nz,nvar_anx, &
nz_ua,mx_ua,nz_tab,nsrc_ua,mx_ua_file, &
anx,xs,ys,zp,tem1,tem2,tem3,tem4,tem5,tem6,tem7, &
nuafil,uafname,srcua, &
stnua,elevua,xua,yua,hgtua,obsua, &
qsrcua,huaqsrc,nlvuatab, &
qobsua,qualua,isrcua,nlevsua, &
rmiss,nobsua,istatus)
END IF
!
!-----------------------------------------------------------------------
!
! MPI issues. Broadcast the number of data points, then the data. If
! there are no data values, we don't waste anyone's time.
!
!-----------------------------------------------------------------------
!
IF (mp_opt > 0) THEN
CALL mpupdatei(nobsua,1)
IF (nobsua > 0) THEN
CALL mpupdatec(stnua,mx_ua*5)
CALL mpupdater(elevua,mx_ua)
CALL mpupdater(xua,mx_ua)
CALL mpupdater(yua,mx_ua)
CALL mpupdater(hgtua,nz_ua*mx_ua)
CALL mpupdater(obsua,nvar_anx*nz_ua*mx_ua)
CALL mpupdater(qobsua,nvar_anx*nz_ua*mx_ua)
CALL mpupdatei(nlvuatab,nsrc_ua)
CALL mpupdatei(qualua,nvar_anx*nz_ua*mx_ua)
CALL mpupdatei(isrcua,mx_ua)
CALL mpupdatei(nlevsua,mx_ua)
!
!-----------------------------------------------------------------------
! We need to know which processor "owns" which obs, so they can be
! consulted on an "as needed" basis.
!-----------------------------------------------------------------------
!
ALLOCATE(item1(nobsua),STAT=istatus)
CALL check_alloc_status(istatus, "adas:item1:nobsua")
ALLOCATE(kua(nprocs),STAT=istatus)
CALL check_alloc_status(istatus, "adas:kua:nprocs")
CALL mpiprocess(nobsua,indexua,nprocs,kua,kuamax, &
isrcua,item1,nx,ny,xua,yua,xs,ys)
DEALLOCATE (item1)
!
!-----------------------------------------------------------------------
! Mark obs that we don't "own" so we don't do unnecessary
! computations.
!-----------------------------------------------------------------------
!
ALLOCATE(item1(0:nprocs-1),STAT=istatus)
CALL check_alloc_status(istatus, "adas:item1:nprocs")
item1 = 0
DO k=1,nmap
IF (mpi_map(k,1) .NE. -1 ) item1(mpi_map(k,1)) = 1
END DO
item1(myproc) = 1
!
!-----------------------------------------------------------------------
! Multi-level obs are never combined, so we don't have to worry.
!-----------------------------------------------------------------------
!
DO k=1,nobsua
IF (indexua(k) == -1) CYCLE
IF (item1(indexua(k)) == 0) isrcua(k) = 0
END DO
DEALLOCATE (item1)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Calculate initial obs differences for upper-air data
!
!-----------------------------------------------------------------------
!
IF (myproc == 0) PRINT *, 'Calling grdtoua'
CALL grdtoua(nx,ny,nxlg,nylg,nz,nz_tab,nz_ua,mx_ua,nvar_anx,nobsua, &
xs,ys,xslg,yslg,zp,anx,qback,hqback,nlvqback, &
tem1,tem2,tem3,tem4,tem5,tem6, &
stnua,isrcua,elevua,xua,yua,hgtua, &
obsua,qobsua,qualua,nlevsua, &
odifua,oanxua,theua,trnua,indexua)
!
!-----------------------------------------------------------------------
! Some data needs to be collected and distributed.
!-----------------------------------------------------------------------
!
IF ( mp_opt > 0 .AND. nobsua > 0 ) THEN
ALLOCATE(tems3dr(nvar_anx,nz_ua,kuamax),STAT=istat)
CALL check_alloc_status(istat, "mpi_2dr_collect:tmps")
ALLOCATE(temr3dr(nvar_anx,nz_ua,kuamax),STAT=istat)
CALL check_alloc_status(istat, "mpi_2dr_collect:tmpr")
ALLOCATE(tems1dr(kuamax),STAT=istat)
CALL check_alloc_status(istat, "mpi_1dr_collect:tmps")
ALLOCATE(temr1dr(kuamax),STAT=istat)
CALL check_alloc_status(istat, "mpi_1dr_collect:tmpr")
ALLOCATE(tems2drua(nz_ua,kuamax),STAT=istat)
CALL check_alloc_status(istat, "mpi_2dr_collect:tmpsua")
ALLOCATE(temr2drua(nz_ua,kuamax),STAT=istat)
CALL check_alloc_status(istat, "mpi_2dr_collect:tmprua")
CALL mpi_3dr_collect( qobsua, nvar_anx, nz_ua, mx_ua, nobsua, indexua, &
nprocs, kua, kuamax, mpi_map, nmap, tems3dr, temr3dr)
CALL mpi_3dr_collect( qualua, nvar_anx, nz_ua, mx_ua, nobsua, indexua, &
nprocs, kua, kuamax, mpi_map, nmap, tems3dr, temr3dr)
CALL mpi_3dr_collect( odifua, nvar_anx, nz_ua, mx_ua, nobsua, indexua, &
nprocs, kua, kuamax, mpi_map, nmap, tems3dr, temr3dr)
CALL mpi_3dr_collect( oanxua, nvar_anx, nz_ua, mx_ua, nobsua, indexua, &
nprocs, kua, kuamax, mpi_map, nmap, tems3dr, temr3dr)
CALL mpi_1dr_collect( trnua, nobsua, indexua, &
nprocs, kua, kuamax, mpi_map, nmap, tems1dr, temr1dr)
CALL mpi_2dr_collect( theua, nz_ua, mx_ua, nobsua, indexua, &
nprocs, kua, kuamax, mpi_map, nmap, tems2drua, temr2drua)
DEALLOCATE(tems1dr)
DEALLOCATE(temr1dr)
DEALLOCATE(tems2drua)
DEALLOCATE(temr2drua)
END IF
!
!-----------------------------------------------------------------------
!
! Read radar data, unfold and convert units
!
!-----------------------------------------------------------------------
!
IF (myproc == 0) 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,trnradc, &
obsrad,oanxrad,odifrad,qobsrad,qualrad, &
ncolrad,istatus,tem1,tem2,tem3,tem4,tem5,tem6)
!
!-----------------------------------------------------------------------
!
! Read retrieval data.
!
! Retrieval code has *NOT* been tested, so it may not be MPI ready!
!
!-----------------------------------------------------------------------
!
IF (myproc == 0) 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)
!
!-----------------------------------------------------------------------
!
! Now get the retrieval data so that each processor only
! contains the data it needs for its local domain plus the radius of
! influence, where applicable.
!
! Not guaranteed MPI ready.
!
!-----------------------------------------------------------------------
!
IF (myproc == 0) 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,trnretc)
!
!-----------------------------------------------------------------------
!
! Quality-control observation differences
!
!-----------------------------------------------------------------------
!
IF (myproc == 0) PRINT *, 'Calling qcdiff'
CALL qcdiff(nvar_anx,nvar_rad,nvar_radin,mx_sng,nsrc_sng, &
indexsng,indexua, &
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)
IF ( mp_opt > 0 .AND. nobsng > 0) THEN
CALL mpi_2dr_collect( qualsng, nvar_anx, mx_sng, nobsng, indexsng, &
nprocs, ksng, ksngmax, mpi_map, nmap, tems2dr, temr2dr)
END IF
IF ( mp_opt > 0 .AND. nobsua > 0 ) THEN
CALL mpi_3dr_collect( qualua, nvar_anx, nz_ua, mx_ua, nobsua, indexua, &
nprocs, kua, kuamax, mpi_map, nmap, tems3dr, temr3dr)
END IF
!
!-----------------------------------------------------------------------
!
! Analyze data
!
!-----------------------------------------------------------------------
!
IF (myproc == 0) THEN
CALL FLUSH (6)
PRINT *, 'Calling anxiter'
END IF
!
! Notice that we are passing "nxlg/nylg" instead of "nx/ny".
!
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, &
mpi_map,nmap, &
xsng,ysng,hgtsng,thesng,trnsng, &
obsng,odifsng,qobsng,qualsng,isrcsng,icatsng,nobsng, &
indexsng,nprocs,ksng,ksngmax, &
indexua,kua,kuamax, &
xua,yua,hgtua,theua,trnua, &
obsua,odifua,qobsua,qualua,isrcua,nlevsua,nobsua, &
elvrad,xradc,yradc, &
distrad,uazmrad,vazmrad,hgtradc,theradc,trnradc, &
obsrad,odifrad,qobsrad,qualrad, &
irad,isrcrad,nlevrad,ncolrad, &
xretc,yretc,hgtretc,theretc,trnretc, &
obsret,odifret,qobsret,qualret, &
iret,isrcret,nlevret,ncolret, &
srcsng,srcua,srcrad,srcret, &
ianxtyp,iusesng,iuseua,iuserad,iuseret, &
xyrange,kpvrsq,wlim,zrange,zwlim, &
thrng,trnropt,trnrcst,trnrng,rngsqi,knt,wgtsum,zsum, &
corsng,corua,corrad,corret, &
oanxsng,oanxua,oanxrad,oanxret, &
anx,tem1,tem2,tem3,ibegin,iend, &
tems2dr,temr2dr,tems3dr,temr3dr, &
istatus)
!
!-----------------------------------------------------------------------
!
! 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)
IF (myproc == 0) THEN
WRITE (status_file, '(2A)') runname(1:lfnkey),'.adasstat'
PRINT *, 'Creating ', TRIM(status_file)
OPEN (UNIT=9, FILE=status_file)
END IF
!
! MPI is slightly more complicated. Values for stations not owned by the
! processor may say valid when then aren't.
!
n_used_sng = 0
DO jsta=1,nobsng
IF (mp_opt > 0) THEN
IF (indexsng(jsta) .NE. myproc ) qualsng(1,jsta) = -9
END IF
IF (isrcsng(jsta) .GT. 0 .AND. qualsng(1,jsta) .GT. 0) THEN
n_used_sng(isrcsng(jsta)) = n_used_sng(isrcsng(jsta)) + 1
END IF
END DO
!
! Combine all mesonet and sao types into a single counter for
! each group. KB 5/14/03
!
nmeso=0
nsao=0
DO jsrc=1,nsrc_sng
IF (mp_opt > 0) THEN
CALL mpsumi(n_used_sng(jsrc),1)
END IF
IF (srcsng(jsrc)(:2) .EQ. 'SA') THEN
nsao = nsao + n_used_sng(jsrc)
ELSE IF (srcsng(jsrc)(:4) .EQ. 'ASOS') THEN
nsao = nsao + n_used_sng(jsrc)
ELSE IF (srcsng(jsrc)(:4) .EQ. 'AWOS') THEN
nsao = nsao + n_used_sng(jsrc)
ELSE IF (srcsng(jsrc)(:5) .EQ. 'SYNOP') THEN
nsao = nsao + n_used_sng(jsrc)
ELSE IF (srcsng(jsrc)(:4) .EQ. 'MESO') THEN
nmeso = nmeso + n_used_sng(jsrc)
ELSE IF (srcsng(jsrc)(:5) .EQ. 'WTXMN') THEN
nmeso = nmeso + n_used_sng(jsrc)
ELSE IF (srcsng(jsrc)(:5) .EQ. 'ARMMN') THEN
nmeso = nmeso + n_used_sng(jsrc)
ELSE IF (srcsng(jsrc)(:7) .EQ. 'COAGMET') THEN
nmeso = nmeso + n_used_sng(jsrc)
ELSE IF (srcsng(jsrc)(:7) .EQ. 'HPLAINS') THEN
nmeso = nmeso + n_used_sng(jsrc)
ELSE IF (srcsng(jsrc)(:4) .EQ. 'ISWS') THEN
nmeso = nmeso + n_used_sng(jsrc)
END IF
END DO
IF (myproc == 0) THEN
WRITE (9,9901) '$sao_n_used = ', nsao
WRITE (9,9901) '$meso_n_used = ', nmeso
END IF
DO jsrc=1,nsrc_sng
IF (srcsng(jsrc)(:5) .EQ. 'MDCRS') THEN
IF (myproc == 0) THEN
WRITE (9,9901) '$mdcrs_n_used = ', n_used_sng(jsrc)
END IF
ELSE IF (srcsng(jsrc)(:4) .EQ. 'BUOY') THEN
IF (myproc == 0) THEN
WRITE (9,9901) '$buoy_n_used = ', n_used_sng(jsrc)
END IF
ELSE IF (srcsng(jsrc)(:4) .NE. 'NULL') THEN
IF (myproc == 0) THEN
WRITE (9,'(3A)') '# Unable to process srcsng: ', &
! srcsng(jsrc), stnsng(jsta,1)
srcsng(jsrc)
END IF
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
IF (myproc == 0) THEN
WRITE (9,*) '# Unable to process: irsrcua=', &
isrcua(jsta), ' ', TRIM(stnua(jsta))
END IF
ELSE IF (srcua(isrcua(jsta)) .EQ. 'NWS RAOB') THEN
IF (myproc == 0) THEN
WRITE (9,9902) 'raob', TRIM(stnua(jsta)), nlevsua(jsta)
END IF
n_used_ua(isrcua(jsta)) = n_used_ua(isrcua(jsta)) + 1
ELSE IF (srcua(isrcua(jsta)) .EQ. 'WPDN PRO') THEN
IF (myproc == 0) THEN
WRITE (9,9902) 'pro', TRIM(stnua(jsta)), nlevsua(jsta)
END IF
n_used_ua(isrcua(jsta)) = n_used_ua(isrcua(jsta)) + 1
ELSE IF (srcua(isrcua(jsta))(:4) .NE. 'NULL') THEN
IF (myproc == 0) THEN
WRITE (9,'(3A)') '# Unable to process: ', &
srcua(isrcua(jsta)), TRIM(stnua(jsta))
END IF
ENDIF
END DO
!deprecated:
! 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 (9,9902) 'raob', TRIM(stnua(jsta)), nlevsua(jsta)
! ELSE IF (srcua(isrcua(jsta)) .EQ. 'WPDN PRO') THEN
! WRITE (9,9902) 'pro', TRIM(stnua(jsta)), nlevsua(jsta)
! ELSE IF (srcua(jsrc)(:4) .NE. 'NULL') THEN
! WRITE (9,'(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
IF (myproc == 0) THEN
WRITE (9,9901) '$raob_n_used = ', n_used_ua(jsrc)
ENDIF
ELSE IF (srcua(jsrc)(:8) .EQ. 'WPDN PRO') THEN
IF (myproc == 0) THEN
WRITE (9,9901) '$pro_n_used = ', n_used_ua(jsrc)
ENDIF
END IF
END DO
9902 FORMAT ('${', A, '}{stn_levels}{"', A, '"} = ', I6, ';')
9901 FORMAT (A,I6,';')
CALL FLUSH (6)
CLOSE (9)
!-----------------------------------------------------------------------
!
! 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 (mp_opt > 0) THEN
CALL mpsendrecv2dew(pprt,nx,ny,nz,ebc,wbc,0,tem1)
CALL mpsendrecv2dns(pprt,nx,ny,nz,nbc,sbc,0,tem1)
CALL mpsendrecv2dew(ptprt,nx,ny,nz,ebc,wbc,0,tem1)
CALL mpsendrecv2dns(ptprt,nx,ny,nz,nbc,sbc,0,tem1)
CALL mpsendrecv2dew(qv,nx,ny,nz,ebc,wbc,0,tem1)
CALL mpsendrecv2dns(qv,nx,ny,nz,nbc,sbc,0,tem1)
END IF
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=-99.
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,indexsng,stnsng,isrcsng,csrcsng,xsng,ysng, &
nxlg,nylg,xslg,yslg, &
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,tem11)
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, &
ebc,wbc,nbc,sbc) ! should be global value. Use local value instead,
! since ADAS is still not parallized.
CALL bcv(nx,ny,nz,0., v, &
pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc, &
ebc,wbc,nbc,sbc) ! should be global value.
CALL lbcw(nx,ny,nz,0.0,w,wcont, &
pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc, &
ebc,wbc,nbc,sbc) ! should be global value.
CALL bcsclr(nx,ny,nz,0.,ptprt,ptprt,ptprt, &
pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc, &
ebc,wbc,nbc,sbc) ! should be global value.
CALL bcsclr(nx,ny,nz,0., qv,qv,qv, &
pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc, &
ebc,wbc,nbc,sbc) ! should be global value.
CALL bcsclr(nx,ny,nz,0., qc,qc,qc, &
pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc, &
ebc,wbc,nbc,sbc) ! should be global value.
CALL bcsclr(nx,ny,nz,0., qr,qr,qr, &
pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc, &
ebc,wbc,nbc,sbc) ! should be global value.
CALL bcsclr(nx,ny,nz,0., qi,qi,qi, &
pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc, &
ebc,wbc,nbc,sbc) ! should be global value.
CALL bcsclr(nx,ny,nz,0., qs,qs,qs, &
pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc, &
ebc,wbc,nbc,sbc) ! should be global value.
CALL bcsclr(nx,ny,nz,0., qh,qh,qh, &
pdteb,pdtwb,pdtnb,pdtsb,ebc,wbc,nbc,sbc,tbc,bbc, &
ebc,wbc,nbc,sbc) ! should be global value.
!
!-----------------------------------------------------------------------
!
! 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(tsfcopt > 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(temrad1(nx,ny,nz))
temrad1 = 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(tem17(nx,ny,nz))
tem17 = 0
ALLOCATE(radbuf(rbufsz))
radbuf = 0.
ALLOCATE(sh(nx,ny))
sh = 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, tsoil(1,1,1,0),qsoil(1,1,1,0),snowdpth, &
radfrc, radsw, rnflx, radswnet, radlwin, &
rsirbm,rsirdf,rsuvbm,rsuvdf, &
cosz, cosss, &
fdirir,fdifir,fdirpar,fdifpar, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8, &
tem9,temrad1,tem11,tem12,tem13,tem14,tem15,tem16, &
radbuf, sh, tem17)
! DTD: placed call to fix_soil_nstp here, to propagate adjusted soil temperatures above
! to the other soil types, if nstyp > 1
! The call to adtsfc above only adjusts the soil temperatures for one soil type
CALL fix_soil_nstyp(nx,ny,nzsoil,1,nstyp,tsoil, &
qsoil,wetcanp)
DEALLOCATE(rsirbm)
DEALLOCATE(rsirdf)
DEALLOCATE(rsuvbm)
DEALLOCATE(rsuvdf)
DEALLOCATE(cosz)
DEALLOCATE(cosss)
DEALLOCATE(fdirir)
DEALLOCATE(fdifir)
DEALLOCATE(fdirpar)
DEALLOCATE(fdifpar)
DEALLOCATE(tem12)
DEALLOCATE(tem13)
DEALLOCATE(tem14)
DEALLOCATE(tem15)
DEALLOCATE(tem16)
DEALLOCATE(radbuf)
DEALLOCATE(sh)
DEALLOCATE(tem17)
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,incrhdfcompr,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:
!
!-----------------------------------------------------------------------
!
IF (mp_opt > 0) THEN
CALL mpsendrecv2dew(u,nx,ny,nz,ebc,wbc,1,tem1)
CALL mpsendrecv2dns(u,nx,ny,nz,nbc,sbc,1,tem1)
CALL mpsendrecv2dew(v,nx,ny,nz,ebc,wbc,2,tem1)
CALL mpsendrecv2dns(v,nx,ny,nz,nbc,sbc,2,tem1)
CALL mpsendrecv2dew(w,nx,ny,nz,ebc,wbc,3,tem1)
CALL mpsendrecv2dns(w,nx,ny,nz,nbc,sbc,3,tem1)
CALL mpsendrecv2dew(ptprt,nx,ny,nz,ebc,wbc,0,tem1)
CALL mpsendrecv2dns(ptprt,nx,ny,nz,nbc,sbc,0,tem1)
CALL mpsendrecv2dew(pprt,nx,ny,nz,ebc,wbc,0,tem1)
CALL mpsendrecv2dns(pprt,nx,ny,nz,nbc,sbc,0,tem1)
CALL mpsendrecv2dew(qv,nx,ny,nz,ebc,wbc,0,tem1)
CALL mpsendrecv2dns(qv,nx,ny,nz,nbc,sbc,0,tem1)
CALL mpsendrecv2dew(qc,nx,ny,nz,ebc,wbc,0,tem1)
CALL mpsendrecv2dns(qc,nx,ny,nz,nbc,sbc,0,tem1)
CALL mpsendrecv2dew(qr,nx,ny,nz,ebc,wbc,0,tem1)
CALL mpsendrecv2dns(qr,nx,ny,nz,nbc,sbc,0,tem1)
CALL mpsendrecv2dew(qi,nx,ny,nz,ebc,wbc,0,tem1)
CALL mpsendrecv2dns(qi,nx,ny,nz,nbc,sbc,0,tem1)
CALL mpsendrecv2dew(qs,nx,ny,nz,ebc,wbc,0,tem1)
CALL mpsendrecv2dns(qs,nx,ny,nz,nbc,sbc,0,tem1)
CALL mpsendrecv2dew(qh,nx,ny,nz,ebc,wbc,0,tem1)
CALL mpsendrecv2dns(qh,nx,ny,nz,nbc,sbc,0,tem1)
END IF
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)
IF (myproc == 0) 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
! blocking inserted for ordering i/o for message passing
DO i=0,nprocs-1,max_fopen
! DO i=0,0,max_fopen
IF(myproc >= i.AND.myproc <= i+max_fopen-1)THEN
CALL dtadump(nx,ny,nz,nzsoil,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,zpsoil, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsoil,qsoil,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx,radswnet,radlwin, &
usflx,vsflx,ptsflx,qvsflx, &
tem3,tem4,tem5)
END IF
IF (mp_opt > 0) CALL mpbarrier
END DO
!
!-----------------------------------------------------------------------
!
! 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)
IF (myproc == 0) &
WRITE(6,'(1x,a,a)') 'History data dump in file ',hdmpfn(1:ldmpf)
grdbas = 0 ! No base state or grid array is dumped.
! blocking inserted for ordering i/o for message passing
DO i=0,nprocs-1,max_fopen
! DO i=0,0,max_fopen
IF(myproc >= i.AND.myproc <= i+max_fopen-1)THEN
CALL dtadump(nx,ny,nz,nzsoil,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,zpsoil, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsoil,qsoil,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx,radswnet,radlwin, &
usflx,vsflx,ptsflx,qvsflx, &
tem3,tem4,tem5)
END IF
IF (mp_opt > 0) CALL mpbarrier
END DO
IF (mp_opt > 0) THEN
IF (myproc == 0) WRITE(6,'(/,5x,a)') '==== ADAS_MPI terminated normally ===='
CALL mpexit(0)
ELSE
WRITE(6,'(/,5x,a)') '==== ADAS terminated normally ===='
STOP
END IF
END PROGRAM arpsdas
!
SUBROUTINE exsufx(fname,lenfnm,suffix,maxsuf,dotloc,lensuf) 4
IMPLICIT NONE
INTEGER :: lenfnm
INTEGER :: maxsuf
CHARACTER (LEN=lenfnm) :: fname
CHARACTER (LEN=maxsuf) :: suffix
INTEGER :: lensuf
INTEGER :: dotloc
!
INTEGER :: i,j
!
dotloc=0
lensuf=0
DO i=1,maxsuf
suffix(i:i)=' '
END DO
DO i=lenfnm,1,-1
IF(fname(i:i) == '.') EXIT
END DO
dotloc=i
IF(dotloc > 0) THEN
lensuf=MIN(maxsuf,lenfnm-i)
DO i=1,lensuf
j=dotloc+i
suffix(i:i)=fname(j:j)
END DO
END IF
RETURN
END SUBROUTINE exsufx
!
SUBROUTINE exfroot(fname,lenfnm,froot,mxroot,lenroot) 3
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
slashloc=i
DO i=slashloc,lenfnm
IF(fname(i) == '.') EXIT
END DO
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, & 3,10
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:
!
! 2/2/06 Kevin Thomas
! MPI the counts.
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'mp.inc'
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 :: nxlg,nylg
INTEGER :: knt1,knt2a,knt2b,knt3,knt4
INTEGER :: is_dup
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
!
! We jump thru some extra hoops to make sure the MPI and non-MPI counts
! are the same. If we don't, overlap points can be counted twice.
!
!-----------------------------------------------------------------------
!
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
is_dup = 0
IF ((mp_opt > 0) .AND. (j >= ny-2) .AND. (loc_y .NE. nproc_y)) is_dup=1
DO i=1,nx-1
IF ((mp_opt > 0) .AND. (i >= nx-2) .AND. (loc_x .NE. nproc_x)) is_dup=1
!
!-----------------------------------------------------------------------
!
! Is convective precip turned on or heavy rain occuring?
!
!-----------------------------------------------------------------------
!
IF (prcrate(i,j,3) > 0. .OR. prcrate(i,j,1) > prctrw) THEN
IF (is_dup == 0) 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
IF (is_dup == 0) 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
IF (is_dup == 0) 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
IF (mp_opt > 0 ) THEN
CALL mpsumi(knt2a,1)
CALL mpsumi(knt2b,1)
CALL mpsumi(knt3,1)
CALL mpsumi(knt4,1)
nxlg = (nx - 3) * nproc_x + 3
nylg = (ny - 3) * nproc_y + 3
knt1=((nxlg-1)*(nylg-1))-knt2a-knt2b-knt3-knt4
ELSE
knt1=((nx-1)*(ny-1))-knt2a-knt2b-knt3-knt4
END IF
IF (myproc == 0 ) THEN
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
END IF
!
CALL iedgfill(icatg,1,nx,1,nx-1, 1,ny,1,ny-1,1,1,1,1)
!
RETURN
END SUBROUTINE setcat