PROGRAM arpsextsnd,116
!
!##################################################################
!##################################################################
!###### ######
!###### PROGRAM EXTSND ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Extracts soundings from ARPS model history data.
! User specifies an ARPS history file and either a lat,lon or
! an x-y location.
!
! Creates a file with the sounding data in an ASCII format.
! The format matches one used by UNIDATA's GEMPAK software for ASCII
! output of sounding data (GEMPAK program SNLIST).
!
! Because of that compatibility, the output file can be
! read by certain sounding plotting and analysis programs on the
! computers at the University of Oklahoma (OU).
!
! The following instructions apply to ARPS users at OU:
!
! After running this program use
! ~rcarpent/bin/skewt -sfc -hodo data_file
!
! where data_file is the output of this program (a name specified
! in the input of this program). An NCARgraphics gmeta file
! is created by skewt.
!
! It is also possible to use this file as input to the standard
! skewt plotting programs on the metgem and Rossby computers.
! See the author for details.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster
! April 1992.
!
! MODIFICATION HISTORY:
!
! October, 1992 (K. Brewster)
! Conversion to ARPS 3.0.
!
! February, 1993 (K. Brewster)
! Additional documentation for ARPS 3.1 release
!
! 4/13/93 (M. Xue)
! Modified to conform to the new data dump format.
!
! 10/19/94 (KB)
! Modified to conform to yet another data dump format.
! Corrections made for change in x,y definition in dump file.
!
! 06/19/95 (KB)
! Modified documentation and output formats to update info.
!
! 03/07/96 (KB)
! Changed calling arguments in FINDLC, moved subroutines to a library.
! Corrected a bug in ymap(ny) calculation.
!
! 11/06/96 (KB)
! Unified interpolation calls with other programs using interpolation.
! Added capability to process multiple soundings in one run.
!
! 2000/05/19 (Gene Bassett)
! Converted to F90, creating allocation and arpsextsnd main
! subroutines.
!
! 2000/07/28 (Ming Xue)
! Converted to F90 free format. Use ALLOCATABLE instead of
! POINTER allocation to avoid double memory usage.
!
! Changed to namelist format input.
!
! 2001/12/05 (Ming Hu)
! Add an output file which is sounding used in
! ARPS as base field like May20.snd
!
! 05/28/2002 (J. Brotzge)
! Added tsoil/qsoil to accomodate new soil scheme.
!
! 1 June 2002 Eric Kemp
! Soil variable updates.
!
! 2005/03/30 (Kevin W. Thomas)
! MPI this program so that large domains get use splitfiles and still
! get soundings.
!
! 04/10/2005 (Keith Brewster)
! Added vertical velocity variable to GEMPAK sounding file to
! accomodate NWS NSHARP program.
!
!-----------------------------------------------------------------------
!
! DATA ARRAYS READ IN:
!
! x x coordinate of grid points in physical/comp. space (m)
! y y coordinate of grid points in physical/comp. space (m)
! z z coordinate of grid points in computational space (m)
! zp z coordinate of grid points in physical space (m)
! zpsoil z coordinate of grid points in the soil (m)
!
! uprt Perturbation x component of velocity (m/s)
! vprt Perturbation y component of velocity (m/s)
! wprt Perturbation z component of velocity (m/s)
!
! ptprt Perturbation potential temperature (K)
! pprt Perturbation pressure (Pascal)
!
! qvprt Perturbation water vapor mixing ratio (kg/kg)
! qc Cloud water mixing ratio (kg/kg)
! qr Rainwater mixing ratio (kg/kg)
! qi Cloud ice mixing ratio (kg/kg)
! qs Snow mixing ratio (kg/kg)
! qh Hail mixing ratio (kg/kg)
!
! tke Turbulent Kinetic Energy ((m/s)**2)
! kmh Horizontal turb. mixing coef. for momentum ( m**2/s )
! kmv Vertical turb. mixing coef. for momentum ( m**2/s )
!
! ubar Base state x velocity component (m/s)
! vbar Base state y velocity component (m/s)
! wbar Base state z velocity component (m/s)
! ptbar Base state potential temperature (K)
! pbar Base state pressure (Pascal)
! rhobar Base state density (kg/m**3)
! qvbar Base state water vapor mixing ratio (kg/kg)
!
! soiltyp Soil type
! vegtyp Vegetation type
! lai Leaf Area Index
! roufns Surface roughness
! veg Vegetation fraction
!
! tsoil Soil temperature (K)
! qsoil Soil moisture (m**3/m**3)
! wetcanp Canopy water amount
!
! raing Grid supersaturation rain
! rainc Cumulus convective rain
! prcrate Precipitation rates
!
! radfrc Radiation forcing (K/s)
! radsw Solar radiation reaching the surface
! rnflx Net radiation flux absorbed by surface
! radswnet Net shortwave radiation, SWin - SWout
! radlwin Incoming longwave radiation
!
! usflx Surface flux of u-momentum (kg/(m*s**2))
! vsflx Surface flux of v-momentum (kg/(m*s**2))
! ptsflx Surface heat flux (K*kg/(m**2 * s ))
! qvsflx Surface moisture flux of (kg/(m**2 * s))
!
! CALCULATED DATA ARRAYS:
!
! su Sounding x component of velocity (m/s)
! sv Sounding y component of velocity (m/s)
! sw Sounding w component of velocity (m/s)
! stheta Sounding potential temperature (K)
! spres Sounding pressure (mb)
! stemp Sounding temperature (C)
! sdewp Sounding dew-point (C)
! sdrct Sounding wind direction (degrees)
! ssped Sounding wind speed (m/s)
! somega Sounding omega vertical velocity (Pa/s)
! shght Sounding height (m)
!
! WORK ARRAYS:
!
! tem1 Temporary work array.
!
! Temporary arrays are defined and used differently by each
! subroutine.
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz,nzsoil
INTEGER :: nxlg,nylg
INTEGER :: nstyps
INTEGER, PARAMETER :: maxsnd = 1000
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
INCLUDE 'mp.inc'
!
!-----------------------------------------------------------------------
!
! Arrays to be read in:
!
!-----------------------------------------------------------------------
!
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(:,:,:) ! The physical height coordinate defined at
! w-point of the soil grid.
REAL, ALLOCATABLE :: uprt (:,:,:) ! Perturbation u-velocity (m/s)
REAL, ALLOCATABLE :: vprt (:,:,:) ! Perturbation v-velocity (m/s)
REAL, ALLOCATABLE :: wprt (:,:,:) ! Perturbation w-velocity (m/s)
REAL, ALLOCATABLE :: ptprt (:,:,:) ! Perturbation potential temperature (K)
REAL, ALLOCATABLE :: pprt (:,:,:) ! Perturbation pressure (Pascal)
REAL, ALLOCATABLE :: qvprt (:,:,:) ! Perturbation 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 :: kmh (:,:,:) ! Horizontal turb. mixing coef. for
! momentum. ( m**2/s )
REAL, ALLOCATABLE :: kmv (:,:,:) ! Vertical turb. mixing coef. for
! momentum. ( m**2/s )
REAL, ALLOCATABLE :: ubar (:,:,:) ! Base state u-velocity (m/s)
REAL, ALLOCATABLE :: vbar (:,:,:) ! Base state v-velocity (m/s)
REAL, ALLOCATABLE :: wbar (:,:,:) ! Base state w-velocity (m/s)
REAL, ALLOCATABLE :: ptbar (:,:,:) ! Base state potential temperature (K)
REAL, ALLOCATABLE :: pbar (:,:,:) ! Base state pressure (Pascal)
REAL, ALLOCATABLE :: rhobar (:,:,:) ! Base state air density (kg/m**3)
REAL, ALLOCATABLE :: qvbar (:,:,:) ! Base state water vapor specific
! humidity (kg/kg)
INTEGER, ALLOCATABLE :: soiltyp (:,:,:) ! Soil type
REAL, ALLOCATABLE :: stypfrct(:,:,:) ! Soil type
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 (m**3/m**3)
REAL, ALLOCATABLE :: wetcanp(:,:,:) ! Canopy water amount
REAL, ALLOCATABLE :: snowdpth(:,:) ! Snow depth (m)
REAL, ALLOCATABLE :: raing(:,:) ! Grid supersaturation rain
REAL, ALLOCATABLE :: rainc(:,:) ! Cumulus convective rain
REAL, ALLOCATABLE :: prcrate(:,:,:) ! precipitation rate (kg/(m**2*s))
! prcrate(1,1,1) = total precip. rate
! prcrate(1,1,2) = grid scale precip. rate
! prcrate(1,1,3) = cumulus precip. rate
! prcrate(1,1,4) = microphysics precip. rate
REAL, ALLOCATABLE :: radfrc(:,:,:) ! Radiation forcing (K/s)
REAL, ALLOCATABLE :: radsw (:,:) ! Solar radiation reaching the surface
REAL, ALLOCATABLE :: rnflx (:,:) ! Net radiation flux absorbed by surface
REAL, ALLOCATABLE :: radswnet(:,:) ! Net shortwave radiation
REAL, ALLOCATABLE :: radlwin(:,:) ! Incoming longwave radiation
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)
!
!-----------------------------------------------------------------------
!
! Computed variables
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: xs(:) ! x location of scalar points
REAL, ALLOCATABLE :: ys(:) ! y location of scalar points
!
!-----------------------------------------------------------------------
!
! Extracted sounding variables
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: su(:,:)
REAL, ALLOCATABLE :: sv(:,:)
REAL, ALLOCATABLE :: sw(:,:)
REAL, ALLOCATABLE :: stheta(:,:)
REAL, ALLOCATABLE :: sqv(:,:)
REAL, ALLOCATABLE :: spres(:,:)
REAL, ALLOCATABLE :: stemp(:,:)
REAL, ALLOCATABLE :: sdewp(:,:)
REAL, ALLOCATABLE :: sdrct(:,:)
REAL, ALLOCATABLE :: ssped(:,:)
REAL, ALLOCATABLE :: somega(:,:)
REAL, ALLOCATABLE :: shght(:,:)
!
!-----------------------------------------------------------------------
!
! Work Arrays
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: dxfld(:)
REAL, ALLOCATABLE :: dyfld(:)
REAL, ALLOCATABLE :: rdxfld(:)
REAL, ALLOCATABLE :: rdyfld(:)
REAL, ALLOCATABLE :: slopey(:,:,:)
REAL, ALLOCATABLE :: alphay(:,:,:)
REAL, ALLOCATABLE :: betay(:,:,:)
REAL, ALLOCATABLE :: tem1(:,:,:)
REAL, ALLOCATABLE :: tem2(:,:,:)
REAL, ALLOCATABLE :: tem3(:,:,:)
REAL :: xpt(maxsnd),ypt(maxsnd)
CHARACTER (LEN=256) :: ofile(maxsnd)
INTEGER :: ipt(maxsnd),jpt(maxsnd)
INTEGER :: is_good(maxsnd)
!
!-----------------------------------------------------------------------
!
! Sounding indentification variables
! These are required for the proper operation of the
! plotting programs on the metgem computer at the
! University of Oklahoma.
!
!-----------------------------------------------------------------------
!
!wdt update
CHARACTER (LEN=8) :: stid(maxsnd)
REAL :: slat(maxsnd),slon(maxsnd),selev(maxsnd)
INTEGER :: istnm(maxsnd)
!
!-----------------------------------------------------------------------
!
! Map projection variables
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: xmap(:)
REAL, ALLOCATABLE :: ymap(:)
REAL, ALLOCATABLE :: latgr(:,:)
REAL, ALLOCATABLE :: longr(:,:)
INTEGER :: i4time,iyr,imo,iday,ihr,imin,isec
REAL :: latnot(2)
!
!-----------------------------------------------------------------------
!
! Misc. internal variables
!
!-----------------------------------------------------------------------
!
REAL :: time,xctr,yctr,xll,yll,xsndmap,ysndmap
REAL :: ustorm,vstorm
REAL :: xmin, xmax, ymin, ymax
INTEGER :: i,j,k,locopt,isnd,nsnd
INTEGER :: ireturn,hinfmt,lengbf,lenfil,nchin
CHARACTER (LEN=256) :: filename,grdbasfn
CHARACTER (LEN=80) :: sndtime,snddate,sndlocation
NAMELIST /message_passing/ nproc_x,nproc_y,readsplit
NAMELIST /input/ hinfmt,grdbasfn,filename,ustorm,vstorm, &
locopt,nsnd,slat,slon,xpt,ypt,ofile,stid,istnm
INTEGER :: nlevs, istatus
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
CALL mpinit_proc
IF(myproc == 0) WRITE(6,'(/6(/1x,a)/)') &
'###############################################################',&
'# #',&
'# Welcome to EXTSND, a program that reads in a history file #',&
'# generated by ARPS and extracts a sounding at an x-y point. #',&
'# #',&
'###############################################################'
!
!-----------------------------------------------------------------------
!
! Read in message passing options.
!
!-----------------------------------------------------------------------
!
IF (myproc == 0) THEN
READ (5,message_passing)
WRITE(6,'(a)')' Namelist block message_passing sucessfully read.'
END IF
CALL mpupdatei(nproc_x,1)
CALL mpupdatei(nproc_y,1)
CALL mpupdatei(readsplit,1)
IF (mp_opt == 0 ) THEN
nproc_x = 1
nproc_y = 1
readsplit = 0
END IF
!
!-----------------------------------------------------------------------
!
! Initialize message passing variables.
!
!-----------------------------------------------------------------------
!
CALL mpinit_var
xpt = -999.0
ypt = -999.0
IF (myproc == 0 ) THEN
READ(5,input)
lengbf=LEN_trim(grdbasfn)
WRITE(6,'(/a,a)')' The grid/base name is ', grdbasfn(1:lengbf)
lenfil = LEN_trim(filename)
WRITE(6,'(a,a)')' The data set name is ', filename(1:lenfil)
IF( nsnd > maxsnd ) then
WRITE(6,'(a,/a,i5)') &
'The number of sounding locations to be extracted exceeded maximum ',&
'allowed. nsnd is reset to ', maxsnd
nsnd = maxsnd
ENDIF
ENDIF
CALL mpupdatec(grdbasfn,256)
CALL mpupdatec(filename,256)
CALL mpupdatei(lengbf,1)
CALL mpupdatei(lenfil,1)
CALL mpupdatei(hinfmt,1)
CALL mpupdater(ustorm,1)
CALL mpupdater(vstorm,1)
CALL mpupdatei(locopt,1)
CALL mpupdatei(nsnd,1)
CALL mpupdater(slat,maxsnd)
CALL mpupdater(slon,maxsnd)
CALL mpupdatec(ofile,256*maxsnd)
CALL mpupdatec(stid,8*maxsnd)
CALL mpupdatei(istnm,maxsnd)
!
!-----------------------------------------------------------------------
!
! Obtain the grid dimensions from input data.
!
!-----------------------------------------------------------------------
!
IF (mp_opt > 0 .AND. readsplit == 0) THEN
WRITE(filename,'(a,a,2i2.2)') filename(1:lenfil),'_',loc_x,loc_y
lenfil = lenfil + 5
WRITE(grdbasfn,'(a,a,2i2.2)') grdbasfn(1:lengbf),'_',loc_x,loc_y
lengbf = lengbf + 5
ENDIF
IF (myproc == 0) THEN
CALL get_dims_from_data(hinfmt,grdbasfn(1:lengbf), &
nx,ny,nz,nzsoil,nstyps, ireturn)
IF (mp_opt > 0 .AND. readsplit > 0) THEN
!
! Fiddle with nx/ny, which apparently are wrong.
!
nx = (nx - 3) / nproc_x + 3
ny = (ny - 3) / nproc_y + 3
END IF
IF (nstyps <= 0) nstyps = 1
nstyp = nstyps ! Copy to global variable
IF( ireturn /= 0 ) THEN
PRINT*,'Problem occured when trying to get dimensions from data.'
PRINT*,'Program stopped.'
CALL arpsstop('Problem with data.',1)
END IF
END IF
CALL mpupdatei(nx,1)
CALL mpupdatei(ny,1)
CALL mpupdatei(nz,1)
CALL mpupdatei(nzsoil,1)
CALL mpupdatei(nstyps,1)
CALL mpupdatei(nstyp,1)
IF (myproc == 0 ) &
WRITE(6,'(4(a,i5))') 'nx =',nx,', ny=',ny,', nz=',nz,', nzsoil=',nzsoil
!
!-----------------------------------------------------------------------
!
! Allocate arrays
!
!-----------------------------------------------------------------------
!
ALLOCATE(x (nx),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:x')
ALLOCATE(y (ny),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:y')
ALLOCATE(z (nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:z')
ALLOCATE(zp (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:zp')
ALLOCATE(zpsoil (nx,ny,nzsoil),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:zpsoil')
ALLOCATE(uprt (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:uprt')
ALLOCATE(vprt (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:vprt')
ALLOCATE(wprt (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:wprt')
ALLOCATE(ptprt (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:ptprt')
ALLOCATE(pprt (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:pprt')
ALLOCATE(qvprt (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:qvprt')
ALLOCATE(qc (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:qc')
ALLOCATE(qr (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:qr')
ALLOCATE(qi (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:qi')
ALLOCATE(qs (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:qs')
ALLOCATE(qh (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:qh')
ALLOCATE(tke (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:tke')
ALLOCATE(kmh (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:kmh')
ALLOCATE(kmv (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:kmv')
ALLOCATE(ubar (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:ubar')
ALLOCATE(vbar (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:vbar')
ALLOCATE(wbar (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:wbar')
ALLOCATE(ptbar (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:ptbar')
ALLOCATE(pbar (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:pbar')
ALLOCATE(rhobar (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:rhobar')
ALLOCATE(qvbar (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:qvbar')
ALLOCATE(soiltyp (nx,ny,nstyps),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:soiltyp')
ALLOCATE(stypfrct(nx,ny,nstyps),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:stypfrct')
ALLOCATE(vegtyp(nx,ny),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:vegtyp')
ALLOCATE(lai (nx,ny),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:lai')
ALLOCATE(roufns (nx,ny),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:roufns')
ALLOCATE(veg (nx,ny),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:veg')
ALLOCATE(tsoil (nx,ny,nzsoil,0:nstyps),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:tsoil')
ALLOCATE(qsoil (nx,ny,nzsoil,0:nstyps),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:qsoil')
ALLOCATE(wetcanp(nx,ny,0:nstyps),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:wetcanp')
ALLOCATE(snowdpth(nx,ny),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:snowdpth')
ALLOCATE(raing(nx,ny),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:raing')
ALLOCATE(rainc(nx,ny),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:rainc')
ALLOCATE(prcrate(nx,ny,4),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:prcrate')
ALLOCATE(radfrc(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:radfrc')
ALLOCATE(radsw (nx,ny),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:radsw')
ALLOCATE(rnflx (nx,ny),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:rnflx')
ALLOCATE(radswnet(nx,ny),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:radswnet')
ALLOCATE(radlwin(nx,ny),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:radlwin')
ALLOCATE(usflx (nx,ny),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:usflx')
ALLOCATE(vsflx (nx,ny),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:vsflx')
ALLOCATE(ptsflx(nx,ny),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:ptsflx')
ALLOCATE(qvsflx(nx,ny),stat=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:qvsflx')
ALLOCATE(xs(nx),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:xs')
ALLOCATE(ys(ny),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:ys')
ALLOCATE(su(nz,maxsnd),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:su')
ALLOCATE(sv(nz,maxsnd),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:sv')
ALLOCATE(sw(nz,maxsnd),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:sw')
ALLOCATE(stheta(nz,maxsnd),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:stheta')
ALLOCATE(sqv(nz,maxsnd),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:sqv')
ALLOCATE(spres(nz,maxsnd),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:spres')
ALLOCATE(stemp(nz,maxsnd),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:stemp')
ALLOCATE(sdewp(nz,maxsnd),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:sdewp')
ALLOCATE(sdrct(nz,maxsnd),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:sdrct')
ALLOCATE(ssped(nz,maxsnd),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:ssped')
ALLOCATE(somega(nz,maxsnd),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:somega')
ALLOCATE(shght(nz,maxsnd),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:shght')
ALLOCATE(dxfld(nx),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:dxfld')
ALLOCATE(dyfld(ny),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:dyfld')
ALLOCATE(rdxfld(nx),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:rdxfld')
ALLOCATE(rdyfld(ny),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:rdyfld')
ALLOCATE(slopey(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:slopey')
ALLOCATE(alphay(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:alphay')
ALLOCATE(betay(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:betay')
ALLOCATE(tem1(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:tem1')
ALLOCATE(tem2(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:tem2')
ALLOCATE(tem3(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:tem3')
ALLOCATE(xmap(nx),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:xmap')
ALLOCATE(ymap(ny),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:ymap')
ALLOCATE(latgr(nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:latgr')
ALLOCATE(longr(nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, 'arpsextsnd:longr')
x =0.0
y =0.0
z =0.0
zp =0.0
zpsoil =0.0
uprt =0.0
vprt =0.0
wprt =0.0
ptprt =0.0
pprt =0.0
qvprt =0.0
qc =0.0
qr =0.0
qi =0.0
qs =0.0
qh =0.0
tke =0.0
kmh =0.0
kmv =0.0
ubar =0.0
vbar =0.0
wbar =0.0
ptbar =0.0
pbar =0.0
radsw =0.0
qvbar =0.0
soiltyp =0
stypfrct=0.0
vegtyp=0.0
lai =0.0
roufns =0.0
veg =0.0
tsoil =0.0
qsoil =0.0
wetcanp=0.0
snowdpth=0.0
raing=0.0
rainc=0.0
prcrate=0.0
radfrc=0.0
radsw =0.0
rnflx =0.0
radswnet =0.0
radlwin =0.0
usflx =0.0
vsflx =0.0
ptsflx=0.0
qvsflx=0.0
xs=0.0
ys=0.0
su=0.0
sv=0.0
sw=0.0
stheta=0.0
sqv=0.0
spres=0.0
stemp=0.0
sdewp=0.0
sdrct=0.0
ssped=0.0
somega=0.0
shght=0.0
dxfld=0.0
dyfld=0.0
rdxfld=0.0
rdyfld=0.0
slopey=0.0
alphay=0.0
betay=0.0
tem1=0.0
tem2=0.0
tem3=0.0
!
!-----------------------------------------------------------------------
!
! Read all input data arrays
!
!-----------------------------------------------------------------------
!
CALL dtaread(nx,ny,nz,nzsoil,nstyps, &
hinfmt,nchin,grdbasfn(1:lengbf),lengbf, &
filename(1:lenfil),lenfil,time, &
x,y,z,zp,zpsoil, uprt ,vprt ,wprt ,ptprt, pprt , &
qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv, &
ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsoil,qsoil,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx,radswnet,radlwin, &
usflx,vsflx,ptsflx,qvsflx, &
ireturn, tem1,tem2,tem3)
!
!-----------------------------------------------------------------------
!
! ireturn = 0 for a successful read
!
!-----------------------------------------------------------------------
!
IF( ireturn == 0 ) THEN ! successful read
!
!-----------------------------------------------------------------------
!
! Calculate scalar locations
!
!-----------------------------------------------------------------------
!
DO i=1,nx-1
xs(i)=0.5*(x(i)+x(i+1))
END DO
DO j=1,ny-1
ys(j)=0.5*(y(j)+y(j+1))
END DO
dx=x(2)-x(1)
dy=y(2)-y(1)
!
!-----------------------------------------------------------------------
!
! Compute sounding time and write it
!
!-----------------------------------------------------------------------
!
CALL ctim2abss( year,month,day,hour,minute,second,i4time)
i4time=i4time + nint(time)
CALL abss2ctim( i4time, iyr, imo, iday, ihr, imin, isec )
WRITE(6,'(a,i4.4,a1,i2.2,a1,i2.2,a1,i2.2,a1,i2.2,a1,i2.2)') &
' Time of history data: ', &
iyr,'/',imo,'/',iday,':',ihr,':',imin,':',isec
!
!-----------------------------------------------------------------------
!
! Find and write the lot, lon of the extraction point in
! SNLIST file.
!
!-----------------------------------------------------------------------
!
latnot(1)=trulat1
latnot(2)=trulat2
CALL setmapr(mapproj,sclfct,latnot,trulon)
CALL lltoxy(1,1,ctrlat,ctrlon,xctr,yctr)
PRINT *, ' dx= ',dx,' dy= ',dy
nxlg = (nx-3)*nproc_x+3
nylg = (ny-3)*nproc_y+3
xll=xctr-(0.5*(nxlg-3)*dx)
yll=yctr-(0.5*(nylg-3)*dy)
DO i=1,nx-1
xmap(i)=xll+xs(i)
END DO
xmap(nx)=2.*xmap(nx-1)-xmap(nx-2)
DO j=1,ny-1
ymap(j)=yll+ys(j)
END DO
ymap(ny)=2.*ymap(ny-1)-ymap(ny-2)
CALL xytoll(nx,ny,xmap,ymap,latgr,longr)
!
! Set the range for checking the soundings. Make sure that more than one
! processor doesn't check a point, as each processor writes out its own
! files!
!
xmin = xs(2)
xmax = xs(nx-1)
ymin = ys(2)
ymax = ys(ny-1)
DO isnd=1,nsnd
IF(locopt == 1) THEN
IF(slat(isnd) < -90.) EXIT
CALL lltoxy(1,1,slat(isnd),slon(isnd),xpt(isnd),ypt(isnd))
xpt(isnd)=xpt(isnd)-xll
ypt(isnd)=ypt(isnd)-yll
ELSE
xpt(isnd)=xpt(isnd)*1000.
ypt(isnd)=ypt(isnd)*1000.
xsndmap=xpt(isnd)+xll
ysndmap=ypt(isnd)+yll
CALL xytoll(1,1,xsndmap,ysndmap,slat(isnd),slon(isnd))
END IF
print*,' isnd, xpt, ypt =', isnd, xpt(isnd), ypt(isnd)
!
! Check to see if the sounding is out of the grid box. If so, computed
! data will be garbage!
!
is_good(isnd)=1
IF( xpt(isnd) < xmin .or. xpt(isnd) > xmax .or. &
ypt(isnd) < ymin .or. ypt(isnd) > ymax ) THEN
write(6,*)
write(6,*) 'Station ',stid(isnd),' is outside of the grid.'
write(6,*)
is_good(isnd)=0
cycle
END IF
!
!-----------------------------------------------------------------------
!
! Find location in ARPS grid.
!
!-----------------------------------------------------------------------
!
CALL setijloc(1,1,nx,ny,xpt(isnd),ypt(isnd), &
xs,ys,ipt(isnd),jpt(isnd))
!
WRITE (6,'(2x,a,f12.2,a,f12.2,/2X,a,f12.2,a,f12.2,/2X,a,i6,a,i6)')&
' location x= ',(0.001*xpt(isnd)),' y= ',(0.001*ypt(isnd)), &
' lat= ',slat(isnd),' lon= ',slon(isnd), &
' found near pt i= ',ipt(isnd),' j= ',jpt(isnd)
WRITE (6,'(12x,a,f12.2,a,f12.2,/10X,a,f12.2,a,f12.2)') &
' x= ',(0.001*xs(ipt(isnd))), &
' y= ',(0.001*ys(jpt(isnd))), &
' lat= ',latgr(ipt(isnd),jpt(isnd)), &
' lon= ',longr(ipt(isnd),jpt(isnd))
IF(ofile(isnd)(1:4) == ' ') ofile(isnd)='SNLIST.FIL'
END DO
!
!-----------------------------------------------------------------------
!
! Interpolate (in the horizontal) for the whole vertical column.
!
!-----------------------------------------------------------------------
!
CALL setdxdy(nx,ny, &
1,nx-1,1,ny-1, &
xs,ys,dxfld,dyfld,rdxfld,rdyfld)
CALL colint(nx,ny,nz,maxsnd,nsnd, &
xs,ys,zp,xpt,ypt,ipt,jpt, &
uprt, vprt, wprt, ptprt, pprt, qvprt, &
ubar, vbar, wbar, ptbar, pbar, qvbar, is_good, &
su,sv,somega,stheta,spres,shght,sqv,selev, &
dxfld,dyfld,rdxfld,rdyfld, &
slopey,alphay,betay, &
tem1,nlevs)
!
!-----------------------------------------------------------------------
!
! Add back storm motion that ARPS subtracts from the sounding winds.
!
!-----------------------------------------------------------------------
!
DO isnd=1,nsnd
IF (is_good(isnd) == 0 ) cycle
DO k=1,nlevs
su(k,isnd)=su(k,isnd)+ustorm
sv(k,isnd)=sv(k,isnd)+vstorm
END DO
!
!-----------------------------------------------------------------------
!
! Convert sounding to desired units/quantities
!
!-----------------------------------------------------------------------
!
CALL cnvsnd(su(1,isnd),sv(1,isnd),sw(1,isnd),stheta(1,isnd), &
spres(1,isnd),sqv(1,isnd),slon(isnd), &
sdrct(1,isnd),ssped(1,isnd),somega(i,isnd), &
stemp(1,isnd),sdewp(1,isnd),nlevs)
!
!-----------------------------------------------------------------------
!
! Output sounding to look like GEMPAK SNLIST.FIL
! to use the Skew-T and Hodograph programs
! e.g., those by Bill McCaul and Rich Carpenter and NSHARP
!
! Example of what the header looks like:
!
!23456789012345678901234567890123456789012345678901234567890
!SNPARM = PRES;HGHT;TMPC;DWPC;DRCT;SPED;OMEG
!
!STID = SEP STNM = 72260 TIME = 920308/1500
!SLAT = 32.21 SLON = -98.18 SELEV = 399.0
!
! PRES HGHT TMPC DWPC DRCT SPED OMEG
!
!-----------------------------------------------------------------------
!
OPEN(3,FILE=trim(ofile(isnd)),STATUS='unknown')
WRITE(3,'(/a/)') ' SNPARM = PRES;HGHT;TMPC;DWPC;DRCT;SPED;OMEG'
iyr=MOD(iyr,100)
!wdt update
WRITE(3,'(a,a8,3x,a,i8,3x,a,i2.2,i2.2,i2.2,a1,i2.2,i2.2)') &
' STID = ',stid(isnd), &
'STNM = ',istnm(isnd), &
'TIME = ',iyr,imo,iday,'/',ihr,imin
WRITE(3,'(a,f8.2,3x,a,f8.2,3x,a,f7.1/)') &
' SLAT = ',slat(isnd),'SLON = ',slon(isnd), &
'SELV = ',selev(isnd)
WRITE(3,'(6x,a)') &
'PRES HGHT TMPC DWPC DRCT SPED OMEG'
!
!-----------------------------------------------------------------------
!
! Print out of sounding
!
!-----------------------------------------------------------------------
!
DO k=1,nlevs
WRITE(3,'(1x,7(F9.2))') &
spres(k,isnd),shght(k,isnd), &
stemp(k,isnd),sdewp(k,isnd), &
sdrct(k,isnd),ssped(k,isnd),somega(k,isnd)
END DO
CLOSE(3)
write(6,'(a,a,a)') &
'Sounding file ',trim(ofile(isnd)),' was produced.'
!-----------------------------------------------------------------------
!
! OUTPUT sounding according to follow format which can be used in
! ARPS as base field like May20.snd
!
! Sounding File Format for ARPS 3.0
!
! Record 1: a header line, such as "1-D Sounding Input for ARPS"
! (skipped)
! Record 2: miscellaneous description of sounding (skipped)
! Record 3: time of sounding (character*72)
! Record 4: date of sounding (character*72)
! Record 5: location of sounding (character*72)
! Record 6: three character strings designate the sounding
! data type, e.g.
! 'pressure' 'potential_temperature' 'relative_humidity'.
! Only the first character of the strings
! is decoded, and thus the first character should not be
! left blank. Note that either upper or lower case may be
! used. A more detailed explanation is provided in the
! portion of the code where these strings are declared.
! ------------------------------
! The following three strings designate the type of sounding data.
!
! if height(1:1)='h' or 'H', the sounding is given on height levels
! if height(1:1)='p' or 'P', the sounding is given on pressure levels
!
! if therm(1:1) ='t' or 'T', the sounding is specified in temperature
! if therm(1:1) ='p' or 'P', the sounding is specified in potential
! temperature.
!
! if humid(1:1) ='s' or 'S', the soundings uses specific humidity
! if humid(1:1) ='r' or 'R', the sounding uses relative humidity,
! if humid(1:1) ='d' or 'D', the soundings uses dewpoint temperature,
!
! if wind(1:1) ='x' or 'Y', the sounding is specified in x and y
! component of velocity (u and v)
! if wind(1:1) ='d' or 'D', the sounding is specified in direction
! and speed in m/s.
! if wind(1:1) ='k' or 'K', the sounding is specified in direction
! and speed in knots.
!--------------------------------
! Record 7: Ground-level height (m) and the correspoding pressure
! (Pascal) when sounding is specified at height levels.
! When it is given on pressure levels, this record is
! not used. The last sounding data is assumed to be at
! the ground level (z=0).
! Record 8: number of data levels in sounding (variable lvlsnd )
! Record 9: a line of data column labels (not read in)
! Record 10 to the end:
! sounding data in the order 'z/p, pt/t, qv/rh, u, v'
!
! For records 10 to the end, there is one line of data
! corresponding to each sounding level.
!
! Important convention:
! Line 10 corresponds to the level k = lvlsnd
! Line 11 corresponds to the level k = (lvlsnd -1)
! etc.
!
! Units of the data:
! pressure: Pascals, height: meters, temperature, dewpoint,
! and potential temperature: degrees Kelvin, specific
! humidity kg/kg, relative humidity: 0 to 1.
!
! CAUTION: The sounding data have to be listed in the order of
! decreasing height or increasing pressure.
!
!-----------------------------------------------------------------------
!
WRITE(*,*) year,month,day,hour,minute,second
WRITE(sndtime,'(i2,a,i2,a,a)') ihr,':',imin,':','00 CST'
IF(ihr <=9 ) sndtime(1:1)='0'
IF(imin <=9 ) sndtime(4:4)='0'
WRITE(snddate,'(i3,i3,a,i5,a,f6.1,a,f6.1,a)') iday,imo,',',year, &
' - Domain speed plused (umove=',ustorm,',vmove=',ustorm,')'
WRITE(sndlocation,'(a,f8.2,3x,a,f8.2,3x,a,f7.1)') &
'SLAT=',slat(isnd),'SLON=',slon(isnd),'SELV=',selev(isnd)
OPEN(3,FILE=trim(ofile(isnd))//'forARPS',STATUS='unknown')
WRITE(3,*) '1-D Sounding Input for ARPS'
WRITE(3,*) 'supercell storm'
WRITE(3,'(1x,a)') sndtime
WRITE(3,'(a)') snddate
WRITE(3,'(1x,a)') sndlocation
WRITE(3,*) '''height'' ''potential temperature'' ''specific humidity'' ''uv'' '
WRITE(3,*) shght(1,isnd), spres(1,isnd)*100
WRITE(3,*) nlevs
WRITE(3,*) ' height potential temperature SH U V'
!
DO k=nlevs,1,-1
WRITE(3,'(1x,6(F14.5))') &
shght(k,isnd), &
stheta(k,isnd),sqv(k,isnd), &
su(k,isnd),sv(k,isnd)
END DO
CLOSE(3)
write(6,'(a,a,a)') &
'Sounding file ',trim(ofile(isnd))//'forARPS',' was produced.'
!-----------------------------------------------------------------------
END DO ! the number of sounding
ELSE
WRITE(6,'(a)') ' Error reading data. EXTSND ends'
END IF
IF (mp_opt > 0) CALL mpexit(0)
STOP
END PROGRAM ARPSEXTSND
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE COLINT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE colint(nx,ny,nz,maxsnd,nsnd, & 8,9
xs,ys,zp,xpt,ypt,ipt,jpt, &
uprt, vprt, wprt, ptprt, pprt, qvprt, &
ubar, vbar, wbar, ptbar, pbar, qvbar, is_good, &
su,sv,sw,stheta,spres,shght,sqv,selev, &
dxfld,dyfld,rdxfld,rdyfld, &
slopey,alphay,betay, &
tem1,nlevs)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Interpolates ARPS history data in the horizontal to create
! a column of data located at point xpt, ypt.
!
! Bilinear interpolation is used.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster
! April 1992.
!
! MODIFICATION HISTORY:
!
! October, 1992 (K. Brewster)
! Conversion to ARPS 3.0.
!
! October, 1994 (K. Brewster)
! Conversion to ARPS 4.0.
!
! 04/10/2005 (K. Brewster)
! Addition of vertical velocity fields.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx,ny,nz Dimensions of ARPS grids.
!
! xs x coordinate of scalar points in physical/comp. space (m)
! ys y coordinate of scalar points in physical/comp. space (m)
! zp z coordinate of scalar grid points in physical space (m)
!
! xpt x coordinate of desired sounding (m)
! ypt y coordinate of desired sounding (m)
!
! ipt i index of grid point just west of xpt,ypt
! jpt j index of grid point just south of xpt,ypt
!
! uprt x component of perturbation velocity (m/s)
! vprt y component of perturbation velocity (m/s)
! vprt z component of perturbation velocity (m/s)
!
! ptprt Perturbation potential temperature (K)
! pprt Perturbation pressure (Pascal)
!
! qvprt Perturbation water vapor mixing ratio (kg/kg)
!
! ubar Base state x velocity component (m/s)
! vbar Base state y velocity component (m/s)
! wbar Base state z velocity component (m/s)
! ptbar Base state potential temperature (K)
! pbar Base state pressure (Pascal)
! qvbar Base state water vapor mixing ratio (kg/kg)
!
! is_good Make no computations if this sounding is outside the grid.
!
! OUTPUT:
!
! su Interpolated u wind component. (m/s)
! sv Interpolated v wind component. (m/s)
! sw Interpolated vertical velocity. (m/s)
! stheta Interpolated potential temperature (K).
! spres Interpolated pressure. (Pascals)
! shght Interpolated height (meters)
! sqv Interpolated water vapor mixing ratio (kg/kg).
! selev Interpolated surface elevation (m)
! nlevs Number of above-ground sounding levels.
!
! WORK ARRAYS:
!
! tem1 Temporary work array.
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
! Arguments -- location data
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Dimensions of ARPS grids.
INTEGER :: maxsnd
INTEGER :: nsnd
REAL :: xs(nx) ! x coordinate of grid points in
! physical/comp. space (m)
REAL :: ys(ny) ! y coordinate of grid points in
! physical/comp. space (m)
REAL :: zp(nx,ny,nz) ! z coordinate of grid points in
! physical space (m)
REAL :: xpt(maxsnd) ! location to find in x coordinate (m)
REAL :: ypt(maxsnd) ! location to find in y coordinate (m)
INTEGER :: ipt(maxsnd) ! i index to the west of desired
! location
INTEGER :: jpt(maxsnd) ! j index to the south of desired
! location
INTEGER :: is_good(maxsnd) ! Zero when sounding is outside the grid
!
!-----------------------------------------------------------------------
!
! Arguments -- model data
!
!-----------------------------------------------------------------------
!
REAL :: uprt (nx,ny,nz) ! Perturbation u-velocity (m/s)
REAL :: vprt (nx,ny,nz) ! Perturbation v-velocity (m/s)
REAL :: wprt (nx,ny,nz) ! Perturbation w-velocity (m/s)
REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature (K)
REAL :: pprt (nx,ny,nz) ! Perturbation pressure (Pascal)
REAL :: qvprt (nx,ny,nz) ! Perturbation water vapor specific
! humidity (kg/kg)
REAL :: ubar (nx,ny,nz) ! Base state u-velocity (m/s)
REAL :: vbar (nx,ny,nz) ! Base state v-velocity (m/s)
REAL :: wbar (nx,ny,nz) ! Base state w-velocity (m/s)
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal)
REAL :: qvbar (nx,ny,nz) ! Base state water vapor specific
! humidity (kg/kg)
REAL :: dxfld(nx)
REAL :: dyfld(ny)
REAL :: rdxfld(nx)
REAL :: rdyfld(ny)
REAL :: slopey(nx,ny,nz)
REAL :: alphay(nx,ny,nz)
REAL :: betay(nx,ny,nz)
!
!-----------------------------------------------------------------------
!
! Arguments -- Extracted sounding variables
!
!-----------------------------------------------------------------------
!
REAL :: su(nz,maxsnd)
REAL :: sv(nz,maxsnd)
REAL :: sw(nz,maxsnd)
REAL :: stheta(nz,maxsnd)
REAL :: sqv(nz,maxsnd)
REAL :: spres(nz,maxsnd)
REAL :: shght(nz,maxsnd)
REAL :: selev(maxsnd)
INTEGER :: nlevs
!
!-----------------------------------------------------------------------
!
! Temporary work arrays
!
!-----------------------------------------------------------------------
!
REAL :: tem1(nx,ny,nz)
!
!-----------------------------------------------------------------------
!
! Functions called
!
!-----------------------------------------------------------------------
!
REAL :: aint2d
!
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k,kk,isnd
REAL :: w2,w3
REAL :: t1,t2,t3,hmid,tmid,qvsat,rh
INTEGER :: iorder
PARAMETER (iorder = 2)
REAL :: pntint2d
!
!-----------------------------------------------------------------------
!
! Function f_qvsat and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
REAL :: f_qvsat
!fpp$ expand (f_qvsat)
!!dir$ inline always f_qvsat
!*$* inline routine (f_qvsat)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
nlevs=nz-2
CALL setdrvy(nx,ny,nz, &
1,nx-1,1,ny-1,1,nz-1, &
dyfld,rdyfld,zp, &
slopey,alphay,betay)
DO isnd=1,nsnd
if ( is_good(isnd) .eq. 0 ) cycle
DO k=2,nz-1
shght(k,isnd)=pntint2d(nx,ny, &
1,nx-1,1,ny-1, &
iorder,xs,ys,xpt(isnd),ypt(isnd), &
ipt(isnd),jpt(isnd),zp(1,1,k), &
dxfld,dyfld,rdxfld,rdyfld, &
slopey(1,1,k),alphay(1,1,k),betay(1,1,k))
END DO
END DO
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)=ptbar(i,j,k)+ptprt(i,j,k)
END DO
END DO
END DO
CALL setdrvy(nx,ny,nz, &
1,nx-1,1,ny-1,1,nz-1, &
dyfld,rdyfld,tem1, &
slopey,alphay,betay)
DO isnd=1,nsnd
if ( is_good(isnd) .eq. 0 ) cycle
DO k=2,nz-1
stheta(k,isnd)=pntint2d(nx,ny, &
1,nx-1,1,ny-1, &
iorder,xs,ys,xpt(isnd),ypt(isnd), &
ipt(isnd),jpt(isnd),tem1(1,1,k), &
dxfld,dyfld,rdxfld,rdyfld, &
slopey(1,1,k),alphay(1,1,k),betay(1,1,k))
END DO
END DO
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)=pbar(i,j,k)+pprt(i,j,k)
END DO
END DO
END DO
CALL setdrvy(nx,ny,nz, &
1,nx-1,1,ny-1,1,nz-1, &
dyfld,rdyfld,tem1, &
slopey,alphay,betay)
DO isnd=1,nsnd
if ( is_good(isnd) .eq. 0 ) cycle
DO k=2,nz-1
spres(k,isnd)=pntint2d(nx,ny, &
1,nx-1,1,ny-1, &
iorder,xs,ys,xpt(isnd),ypt(isnd), &
ipt(isnd),jpt(isnd),tem1(1,1,k), &
dxfld,dyfld,rdxfld,rdyfld, &
slopey(1,1,k),alphay(1,1,k),betay(1,1,k))
END DO
END DO
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)=qvbar(i,j,k)+qvprt(i,j,k)
END DO
END DO
END DO
CALL setdrvy(nx,ny,nz, &
1,nx-1,1,ny-1,1,nz-1, &
dyfld,rdyfld,tem1, &
slopey,alphay,betay)
DO isnd=1,nsnd
if ( is_good(isnd) .eq. 0 ) cycle
DO k=2,nz-1
sqv(k,isnd)=pntint2d(nx,ny, &
1,nx-1,1,ny-1, &
iorder,xs,ys,xpt(isnd),ypt(isnd), &
ipt(isnd),jpt(isnd),tem1(1,1,k), &
dxfld,dyfld,rdxfld,rdyfld, &
slopey(1,1,k),alphay(1,1,k),betay(1,1,k))
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Get height at scalar points, since zp was defined at w points.
!
!-----------------------------------------------------------------------
!
DO isnd=1,nsnd
if ( is_good(isnd) .eq. 0 ) cycle
selev(isnd)=shght(2,isnd)
DO k=1,nz-1
shght(k,isnd)=0.5*(shght(k+1,isnd)+shght(k,isnd))
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Form total u wind component at scalar points
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)=0.5*(ubar(i,j,k)+ubar(i+1,j,k))+ &
0.5*(uprt(i,j,k)+uprt(i+1,j,k))
END DO
END DO
END DO
CALL setdrvy(nx,ny,nz, &
1,nx-1,1,ny-1,1,nz-1, &
dyfld,rdyfld,tem1, &
slopey,alphay,betay)
DO isnd=1,nsnd
if ( is_good(isnd) .eq. 0 ) cycle
DO k=2,nz-1
su(k,isnd)=pntint2d(nx,ny, &
1,nx-1,1,ny-1, &
iorder,xs,ys,xpt(isnd),ypt(isnd), &
ipt(isnd),jpt(isnd),tem1(1,1,k), &
dxfld,dyfld,rdxfld,rdyfld, &
slopey(1,1,k),alphay(1,1,k),betay(1,1,k))
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Form total v wind component at scalar points
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)=0.5*(vbar(i,j,k)+vbar(i,j+1,k)) + &
0.5*(vprt(i,j,k)+vprt(i,j+1,k))
END DO
END DO
END DO
CALL setdrvy(nx,ny,nz, &
1,nx-1,1,ny-1,1,nz-1, &
dyfld,rdyfld,tem1, &
slopey,alphay,betay)
DO isnd=1,nsnd
if ( is_good(isnd) .eq. 0 ) cycle
DO k=2,nz-1
sv(k,isnd)=pntint2d(nx,ny, &
1,nx-1,1,ny-1, &
iorder,xs,ys,xpt(isnd),ypt(isnd), &
ipt(isnd),jpt(isnd),tem1(1,1,k), &
dxfld,dyfld,rdxfld,rdyfld, &
slopey(1,1,k),alphay(1,1,k),betay(1,1,k))
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Form total w wind component at scalar points
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)=0.5*(wbar(i,j,k)+wbar(i,j,k+1)) + &
0.5*(wprt(i,j,k)+wprt(i,j,k+1))
END DO
END DO
END DO
CALL setdrvy(nx,ny,nz, &
1,nx-1,1,ny-1,1,nz-1, &
dyfld,rdyfld,tem1, &
slopey,alphay,betay)
DO isnd=1,nsnd
if ( is_good(isnd) .eq. 0 ) cycle
DO k=2,nz-1
sw(k,isnd)=pntint2d(nx,ny, &
1,nx-1,1,ny-1, &
iorder,xs,ys,xpt(isnd),ypt(isnd), &
ipt(isnd),jpt(isnd),tem1(1,1,k), &
dxfld,dyfld,rdxfld,rdyfld, &
slopey(1,1,k),alphay(1,1,k),betay(1,1,k))
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Get a value at the surface, by extrapolating from the
! 2nd and third levels.
!
!-----------------------------------------------------------------------
!
DO isnd=1,nsnd
if ( is_good(isnd) .eq. 0 ) cycle
shght(1,isnd)=selev(isnd)
w3=(shght(1,isnd)-shght(2,isnd)) &
/(shght(3,isnd)-shght(2,isnd))
w2=(1.-w3)
su(1,isnd)=w2* su(2,isnd) + w3* su(3,isnd)
sv(1,isnd)=w2* sv(2,isnd) + w3* sv(3,isnd)
sw(1,isnd)=0.5*sw(2,isnd)
IF(stheta(3,isnd) > stheta(2,isnd)) THEN
stheta(1,isnd)=w2*stheta(2,isnd) + w3*stheta(3,isnd)
ELSE
stheta(1,isnd)=stheta(2,isnd)
END IF
!
!-----------------------------------------------------------------------
!
! Integrate downward to get the pressure at level 1.
!
!-----------------------------------------------------------------------
!
t3=stheta(3,isnd)*(spres(3,isnd)/100000.)**rddcp
t2=stheta(2,isnd)*(spres(2,isnd)/100000.)**rddcp
hmid=0.5*(shght(2,isnd)+shght(1,isnd))
tmid=t3+((shght(3,isnd)-hmid)/ &
(shght(3,isnd)-shght(2,isnd)))*(t2-t3)
spres(1,isnd)=spres(2,isnd)* &
EXP(g*(shght(2,isnd)-shght(1,isnd))/(rd*tmid))
!
!-----------------------------------------------------------------------
!
! Use constant RH to extrapolate qv to level 1.
!
!-----------------------------------------------------------------------
!
qvsat = f_qvsat( spres(2,isnd), t2 ) ! saturated S.H.
rh=AMIN1((sqv(2,isnd)/qvsat),1.) ! R.H.
t1=stheta(1,isnd)*(spres(1,isnd)/100000.)**rddcp
qvsat = f_qvsat( spres(1,isnd), t1 )
sqv(1,isnd)=rh*qvsat
END DO
RETURN
END SUBROUTINE colint
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE CNVSND ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE cnvsnd(su,sv,sw,stheta,spres,sqv,slon, & 4,2
sdrct,ssped,somega,stemp,sdewp,nlevs)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Converts units of data extracted from ARPS history data to
! those required of sounding data. Determines direction and
! speed from u and v wind components.
!
! Dew-point formula from Bolton, 1980, MWR pp 1046-1053.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster
! April 1992.
!
! MODIFICATION HISTORY:
!
! October, 1992 (K. Brewster)
! Conversion to ARPS 3.0.
!
! 10/28/1992 (K. Brewster)
! Special allowance for low qv-to-dew pt
!
! 04/10/2005 (K. Brewster)
! Addition of vertical velocity processing.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! su Sounding u wind component. (m/s)
! sv Sounding v wind component. (m/s)
! sw Sounding v wind component. (m/s)
! stheta Sounding potential temperature (K).
! spres Sounding pressure. (Pascals)
! sqv Sounding specific humidity
! slon Sounding longitude
! nlevs Number of above-ground sounding levels.
!
! OUTPUT:
!
! spres Sounding pressure. (mb)
! sdrct Sounding wind direction (degrees from north)
! ssped Sounding wind speed (m/s)
! somega Sounding omega vertical velocity (Pa/s)
! stemp Sounding temperature (degrees C)
! sdewp Sounding dew point temperature (degrees C)
!
!-----------------------------------------------------------------------
!
! Variable declarations
!
!-----------------------------------------------------------------------
!
! Input arguments
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nlevs ! Number of above-ground sounding levels
REAL :: su (nlevs) ! Sounding u wind component (m/s)
REAL :: sv (nlevs) ! Sounding v wind component (m/s)
REAL :: sw (nlevs) ! Sounding w wind component (m/s)
REAL :: stheta(nlevs) ! Sounding potential temperature (K)
REAL :: spres (nlevs) ! Sounding pressure. (Pascals)
REAL :: sqv (nlevs) ! Sounding specific humidity (g/g)
REAL :: slon ! Sounding longitude (degrees E)
!
!-----------------------------------------------------------------------
!
! Output arguments
!
!-----------------------------------------------------------------------
!
REAL :: sdrct(nlevs) ! Sounding wind direction
! (degrees from north)
REAL :: ssped(nlevs) ! Sounding wind speed (m/s)
REAL :: somega(nlevs) ! Sounding verticel velocity (Pa/s)
REAL :: stemp(nlevs) ! Sounding temperature (degrees C)
REAL :: sdewp(nlevs) ! Sounding dew point temperature
! (degrees C)
!
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: k
REAL :: smix,e,bige,alge,rho
!
DO k=1,nlevs
!
!-----------------------------------------------------------------------
!
! Convert u,v to direction and speed
!
!-----------------------------------------------------------------------
!
CALL uvrotdd(1,1,slon,su(k),sv(k),sdrct(k),ssped(k))
!
!-----------------------------------------------------------------------
!
! Convert pressure from Pascals to mb
!
!-----------------------------------------------------------------------
!
spres(k)=spres(k)*0.01
!
!-----------------------------------------------------------------------
!
! Convert theta to temperature in degrees C.
!
!-----------------------------------------------------------------------
!
stemp(k)=stheta(k)*((spres(k)/1000.)**rddcp)
stemp(k)=stemp(k)-273.15
!
!-----------------------------------------------------------------------
!
! Convert w to omega. For simplicity use hydrostatic approximation.
! Use w=dz/dt omega=(dz/dt)*(dp/dz) and hydrostatic relation for dp/dz
! dp/dz=-rho*g and rho=p/(rd*Tv) Tv=T*(1.0+0.61*qv)
!
!-----------------------------------------------------------------------
!
rho=(100.0*spres(k))/(rd*(stemp(k)+273.15)*(1.0+0.61*sqv(k)))
somega(k)=-sw(k)*rho*g
!
!-----------------------------------------------------------------------
!
! Convert specific humidity to dew-point in degrees C.
!
!-----------------------------------------------------------------------
!
IF( sqv(k) > 0.) THEN
smix=sqv(k)/(1.-sqv(k))
e=(spres(k)*smix)/(0.62197 + smix)
bige=e/( 1.001 + ( (spres(k) - 100.) / 900.) * 0.0034)
alge = ALOG(bige/6.112)
sdewp(k)= (alge * 243.5) / (17.67 - alge)
ELSE
sdewp(k)= stemp(k)-30.
END IF
END DO
!
RETURN
END SUBROUTINE cnvsnd