PROGRAM arpsextsnd,133
!
!##################################################################
!##################################################################
!###### ######
!###### 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).
!
! To create skewt plot of the generated sounding:
! bin/skewtpost(or skewtncar) -sfc -hodo data_file
!
! where data_file is the output of this program (a name specified
! in the input of this program). An PostScript or 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.
!
! 07/26/2006 (Yunheng Wang)
! Added ROAB sounding format and extracting multiple ARPS grids
! sounding.
!
!-----------------------------------------------------------------------
!
! 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 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 (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)
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
INTEGER :: scondition
REAL :: valuethres
INTEGER :: outfmt
INTEGER :: nhisfile_max,nhisfile,nfile
PARAMETER (nhisfile_max=200)
CHARACTER (LEN=256) :: hisfile(nhisfile_max)
CHARACTER (LEN=256) :: grdbasfn
CHARACTER (LEN=256) :: ofilehead
CHARACTER (LEN=80) :: sndtime,snddate,sndlocation
CHARACTER (LEN=256) :: oldfile(maxsnd)
INTEGER :: oldfmt
NAMELIST /message_passing/ nproc_x,nproc_y,readsplit
NAMELIST /input/ ustorm,vstorm,locopt,scondition,valuethres, &
nsnd,slat,slon,xpt,ypt,stid,istnm,oldfmt,oldfile
NAMELIST /output/ dirname, ofilehead, outfmt
INTEGER :: nlevs, istatus
INTEGER :: nsndx, nsndy
INTEGER :: xptint, yptint
INTEGER :: ijpt, i1pt,i2pt,i3pt
INTEGER :: ilg, jlg
INTEGER :: FIRST_SND, LAST_SND
INTEGER :: countsnd
CHARACTER(LEN=22), PARAMETER :: sconstr(2) = &
(/ 'Vert. Integ Condensate', 'Composite Reflectivity' /)
CHARACTER(LEN=256) :: ofile, tmpstr
LOGICAL :: GEMPAKSND, ARPSSND, ROABSND
REAL :: zhght
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
CALL mpinit_proc
IF(myproc == 0) WRITE(6,'(/8(/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,'(1x,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
!
!-----------------------------------------------------------------------
!
! Get the names of the input data files.
!
!-----------------------------------------------------------------------
!
IF(myproc == 0) THEN
CALL get_input_file_names(hinfmt,grdbasfn,hisfile,nhisfile)
lengbf = LEN_trim(grdbasfn)
WRITE(6,'(1x,a)') 'Successfully read in ARPS history file names.'
END IF
CALL mpupdatec(grdbasfn,256)
CALL mpupdatec(hisfile,256*nhisfile_max)
CALL mpupdatei(nhisfile,1)
CALL mpupdatei(hinfmt,1)
CALL mpupdatei(lengbf,1)
xpt(:) = -999.0
ypt(:) = -999.0
oldfmt = 0
oldfile(:) = ' '
IF (myproc == 0 ) THEN
READ(5,input)
WRITE(6,'(1x,a)') 'Successfully read in namelist block input.'
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 mpupdater(ustorm,1)
CALL mpupdater(vstorm,1)
CALL mpupdatei(locopt,1)
CALL mpupdatei(scondition,1)
CALL mpupdater(valuethres,1)
CALL mpupdatei(nsnd,1)
CALL mpupdater(slat,maxsnd)
CALL mpupdater(slon,maxsnd)
CALL mpupdater(xpt,maxsnd)
CALL mpupdater(ypt,maxsnd)
CALL mpupdatec(stid,8*maxsnd)
CALL mpupdatei(istnm,maxsnd)
CALL mpupdatei(oldfmt,1)
IF (oldfmt == 1) CALL mpupdatec(oldfile,256*maxsnd)
dirname = './'
outfmt = 0
ofilehead = 'SNLIST'
IF (myproc == 0) THEN
READ(5,output)
WRITE(6,'(1x,a)') 'Successfully read in namelist block output.'
lenfil = LEN_TRIM(dirname)
IF(lenfil > 0) THEN
IF(dirname(lenfil:lenfil) /= '/') THEN
dirname(lenfil+1:lenfil+1) = '/'
lenfil = lenfil + 1
END IF
ELSE
dirname = './'
END IF
END IF
IF (oldfmt == 1) outfmt = 0
CALL mpupdatec(dirname,256)
CALL mpupdatei(outfmt,1)
CALL mpupdatec(ofilehead,256)
GEMPAKSND = .FALSE.
ARPSSND = .FALSE.
ROABSND = .FALSE.
IF (outfmt == 0) THEN
GEMPAKSND = .TRUE.
ARPSSND = .TRUE.
ELSE IF (outfmt == 1) THEN
GEMPAKSND = .TRUE.
ELSE IF (outfmt == 2) THEN
ARPSSND = .TRUE.
ELSE IF (outfmt == 3) THEN
ROABSND = .TRUE.
END IF
!
!-----------------------------------------------------------------------
!
! Obtain the grid dimensions from input data.
!
!-----------------------------------------------------------------------
!
IF (mp_opt > 0 .AND. readsplit == 0) THEN
tmpstr = grdbasfn
WRITE(grdbasfn,'(a,a,2i2.2)') tmpstr(1:lengbf),'_',loc_x,loc_y
lengbf = lengbf + 5
DO nfile = 1,nhisfile
lenfil = LEN_TRIM(hisfile(nfile))
tmpstr = hisfile(nfile)
WRITE(hisfile(nfile),'(a,a,2i2.2)') tmpstr(1:lenfil),'_',loc_x,loc_y
lenfil = lenfil + 5
END DO
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,'(1x,4(a,i5))') 'nx =',nx,', ny=',ny,', nz=',nz,', nzsoil=',nzsoil
IF (locopt == 3) THEN
xptint = INT(xpt(1))
yptint = INT(ypt(1))
nsndx = 0
DO i = 2,nx-2
ilg = (nx-3)*(loc_x-1) + i
IF ( MOD((ilg-2),xptint) == 0 ) THEN
nsndx = nsndx + 1
ipt(nsndx) = i
END IF
END DO
nsndy = 0
DO j = 2,ny-2
jlg = (ny-3)*(loc_y-1)+j
IF ( MOD((jlg-2),yptint) == 0 ) THEN
nsndy = nsndy + 1
jpt((nsndy-1)*nsndx+1) = j
END IF
END DO
nsnd = nsndx*nsndy
IF( nsnd > maxsnd ) then
WRITE(6,'(/1x,a,/2(a,i5)/)') &
'The number of sounding locations to be extracted exceeded maximum ',&
' allowed. nsnd = ', nsnd,' maxsnd = ',maxsnd
CALL arpsstop('Please increase maxsnd in src/arpsextsnd/arpsextsnd.f90.',1)
ENDIF
jpt(2:nsndx) = jpt(1)
DO j = 1,nsndy
ipt((j-1)*nsndx+1) = ipt(1)
END DO
DO j = 2,nsndy
DO i = 2,nsndx
isnd = (j-1)*nsndx + i
ipt(isnd) = ipt(i)
jpt(isnd) = jpt((j-1)*nsndx+1)
END DO
END DO
END IF
!WRITE(0,'(/,a,I2.2,3(a,I4),6(a,I3),/)') &
! 'Rank = ',myproc,': nsnd = ',nsndx,' x ',nsndy,' = ',nsnd, &
! ', (',ipt(1),', ',jpt(1),') -- (',ipt(nsnd),',',jpt(nsnd), &
! '). xptint,yptint = ',xptint,', ',yptint
!
!-----------------------------------------------------------------------
!
! 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
!
!-----------------------------------------------------------------------
!
DO nfile = 1, nhisfile
lenfil=len_trim(hisfile(nfile))
WRITE(6,'(/a,a,a)') &
' Data set ', trim(hisfile(nfile)),' to be read.'
CALL dtaread(nx,ny,nz,nzsoil,nstyps, &
hinfmt,nchin,grdbasfn(1:lengbf),lengbf, &
hisfile(nfile)(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)
tem3(:,:,:) = 0.0
IF (locopt == 3 .AND. scondition == 1) THEN
CALL cal_vic(tem3,qc,qr,qi,qs,qh,rhobar,zp,nx,ny,nz,tem1)
ELSE IF (locopt == 3 .AND. scondition == 2) THEN
CALL temper (nx,ny,nz,ptbar, ptprt, pprt ,pbar,tem2)
CALL reflec_ferrier(nx,ny,nz, rhobar, qr, qs, qh, tem2, tem1)
! CALL reflec(nx,ny,nz, rhobar, qr, qs, qh, tem1)
CALL cal_rfc(nx, ny, nz, tem1, tem3)
END IF
!
! 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)
FIRST_SND = nsnd
LAST_SND = 1
countsnd = 0
WRITE(6,*)
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 IF (locopt == 2) THEN
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))
ELSE IF (locopt == 3) THEN
xpt(isnd) = xs(ipt(isnd))
ypt(isnd) = ys(jpt(isnd))
xsndmap=xpt(isnd)+xll
ysndmap=ypt(isnd)+yll
CALL xytoll(1,1,xsndmap,ysndmap,slat(isnd),slon(isnd))
istnm(isnd) = jpt(isnd)*10000+ipt(isnd)
ilg = (nx-3)*(loc_x-1) + ipt(isnd)
jlg = (ny-3)*(loc_y-1) + jpt(isnd)
ijpt = ilg-2 + (jlg-2)*(nxlg-3)
i1pt = mod(ijpt,26) + ICHAR('A')
i2pt = mod(ijpt/26,26) + ICHAR('A')
i3pt = mod(ijpt/(26*26),26) + ICHAR('A')
WRITE(stid(isnd),'(3a)') CHAR(i3pt),CHAR(i2pt),CHAR(i1pt)
ELSE
WRITE(6,'(/,a,i2,/)') 'ERROR: unknown option locopt = ',locopt
CALL arpsstop('Please check parameter locopt in the namelist input.',1)
END IF
! 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,'(/,2x,a,a6,a,/)') 'Station ',stid(isnd),' is outside of the grid.'
is_good(isnd)=0
CYCLE
END IF
!-----------------------------------------------------------------------
!
! Find location in ARPS grid.
!
!-----------------------------------------------------------------------
!
!
! It is requried that valuethres must > 0 and
! tem3(:,:,:) = 0 when scondition = 0
!
IF (locopt == 3) THEN
xsndmap = tem3(ipt(isnd),jpt(isnd),1)
IF (xsndmap > valuethres) THEN
WRITE(6,'(3(a,I6),4a,G10.2)') '- isnd =',isnd, &
' at (',ipt(isnd),',',jpt(isnd),') was skipped ', &
'because ',sconstr(scondition),' = ',xsndmap
is_good(isnd) = 0
END IF
ELSE
CALL setijloc(1,1,nx,ny,xpt(isnd),ypt(isnd), &
xs,ys,ipt(isnd),jpt(isnd))
END IF
IF (is_good(isnd) == 1) THEN
countsnd = countsnd + 1
WRITE(6,'(a,i6,a,I8.8,a,a3,2(a,F7.2),a)') &
'+ isnd =',isnd,', stnm = ',istnm(isnd),', stid = ',stid(isnd), &
', xpt = ',xpt(isnd)*0.001,' km, ypt = ', ypt(isnd)*0.001,' km'
IF (locopt == 3) THEN
WRITE (6,'(16x,2(a,i6),2(a,f7.2))') &
'at (',ipt(isnd),',',jpt(isnd), &
'), lat = ',slat(isnd),', lon = ',slon(isnd)
ELSE
WRITE (6,'(2x,4(a,f12.2),/2X,a,i6,a,i6,a)') &
' location x= ',(0.001*xpt(isnd)),', y= ',(0.001*ypt(isnd)), &
', lat= ',slat(isnd),', lon= ',slon(isnd), &
' found near point (',ipt(isnd),',',jpt(isnd),')'
WRITE (6,'(12x,4(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))
END IF
FIRST_SND = MIN(FIRST_SND,isnd)
LAST_SND = MAX(LAST_SND,isnd)
END IF
END DO
WRITE(6,'(/,1x,a,I6,/)') 'Total valid soundings = ',countsnd
!
!-----------------------------------------------------------------------
!
! Interpolate (in the horizontal) for the whole vertical column.
!
!-----------------------------------------------------------------------
!
IF (locopt == 3) THEN
CALL coextract(nx,ny,nz,nsnd,zp,ipt,jpt, &
uprt, vprt, wprt, ptprt, pprt, qvprt, &
ubar, vbar, wbar, ptbar, pbar, qvbar, is_good, &
su,sv,somega,stheta,spres,shght,sqv, &
selev,nlevs)
ELSE
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)
END IF
!
!-----------------------------------------------------------------------
!
! 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)
WRITE(ofile,'(a,I4.4,2i2.2,a,3i2.2)') &
TRIM(ofilehead),iyr,imo,iday,'_',ihr,imin,isec
tmpstr = ofile
IF (locopt == 1 .OR. locopt == 2) THEN
! append statid for locopt = 1/2
! because data for each station is in one file
IF (oldfmt == 1) THEN ! Use explicit file names
ofile=oldfile(isnd)
ELSE
WRITE(ofile,'(3a)') TRIM(tmpstr),'_',TRIM(stid(isnd))
END IF
ELSE IF (locopt == 3 .AND. mp_opt > 0) THEN
! append proc. num. for mp_opt >0
! because we want to distinguished file names
WRITE(ofile,'(2a,2I2.2)') TRIM(tmpstr),'_',loc_x,loc_y
END IF
IF (GEMPAKSND) THEN
!
!-----------------------------------------------------------------------
!
! 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
!
!-----------------------------------------------------------------------
!
IF(locopt /= 3 .OR. isnd == FIRST_SND) THEN
IF (oldfmt == 1) THEN
OPEN(3,file=trim(ofile),STATUS='unknown')
ELSE
OPEN(3,FILE=trim(dirname)//trim(ofile)//'.sounding',STATUS='unknown')
END IF
END IF
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
IF (locopt /= 3 .OR. isnd == LAST_SND) THEN
CLOSE(3)
WRITE(6,'(1x,3a,/)') 'Sounding file ',trim(ofile)//'.sounding', &
' was produced.'
END IF
END IF
!-----------------------------------------------------------------------
!
! 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.
!
!-----------------------------------------------------------------------
!
IF (ARPSSND) THEN
! 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)
IF (locopt /= 3 .OR. isnd == FIRST_SND) THEN
IF (oldfmt == 1) THEN
OPEN(4,FILE=trim(ofile)//'forARPS',STATUS='unknown')
ELSE
OPEN(4,FILE=trim(dirname)//trim(ofile)//'.sound',STATUS='unknown')
ENDIF
END IF
WRITE(4,*) '1-D Sounding Input for ARPS'
WRITE(4,*) 'supercell storm'
WRITE(4,'(1x,a)') sndtime
WRITE(4,'(a)') snddate
WRITE(4,'(1x,a)') sndlocation
WRITE(4,*) '''height'' ''potential temperature'' ''specific humidity'' ''uv'' '
WRITE(4,*) shght(1,isnd), spres(1,isnd)*100
WRITE(4,*) nlevs
WRITE(4,*) ' height potential temperature SH U V'
!
DO k=nlevs,1,-1
WRITE(4,'(1x,6(F14.5))') &
shght(k,isnd), &
stheta(k,isnd),sqv(k,isnd), &
su(k,isnd),sv(k,isnd)
END DO
IF (locopt /= 3 .OR. isnd == LAST_SND) THEN
CLOSE(4)
WRITE(6,'(1x,3a,/)') &
'Sounding file ',trim(ofile)//'.sound',' was produced.'
END IF
END IF
!-----------------------------------------------------------------------
IF (ROABSND) THEN
IF(locopt /= 3 .OR. isnd == FIRST_SND) THEN
OPEN(15,FILE=trim(dirname)//trim(ofile)//'.snd',STATUS='unknown')
END IF
WRITE (15,'(i12.8,i12,f11.4,f15.4,f15.0,5x,a3)') &
istnm(isnd),nlevs,slat(isnd),slon(isnd),selev(isnd),stid(isnd)
!
!-----------------------------------------------------------------------
!
! Print out of sounding
!
!-----------------------------------------------------------------------
!
DO k=1,nlevs
zhght = zp(ipt(isnd),jpt(isnd),k)-zp(ipt(isnd),jpt(isnd),2)
IF (locopt == 3 .AND. scondition == 3 .AND. zhght <= valuethres) CYCLE
WRITE(15,'(1x,F12.3,5(F12.4))') &
shght(k,isnd),spres(k,isnd), &
stemp(k,isnd),sdewp(k,isnd), &
sdrct(k,isnd),ssped(k,isnd)
END DO
IF (locopt /= 3 .OR. isnd == LAST_SND) THEN
CLOSE(15)
WRITE(6,'(1x,3a,/)') 'Sounding file ',trim(ofile)//'.snd', &
' was produced.'
END IF
END IF
END DO ! the number of sounding
ELSE
WRITE(6,'(1x,2a,/)') 'Error reading data file ',hisfile(nfile)
END IF
END DO
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, & 9,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
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE COEXTRACT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE coextract(nx,ny,nz,nsnd,zp,ipt,jpt, & 1,2
uprt, vprt, wprt, ptprt, pprt, qvprt, &
ubar, vbar, wbar, ptbar, pbar, qvbar, is_good, &
su,sv,sw,stheta,spres,shght,sqv, &
selev,nlevs)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Extract a column of data from ARPS history data in the horizontal
! located at arrray index ipt, jpt.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yunheng Wang
! July 2006.
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx,ny,nz Dimensions of ARPS grids.
!
! ipt i index of grid point for desired sounding
! jpt j index of grid point for desired sounding
!
! 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 not satified
! the condition.
!
! OUTPUT:
!
! su Extracted u wind component. (m/s)
! sv Extracted v wind component. (m/s)
! sw Extracted vertical velocity. (m/s)
! stheta Extracted potential temperature (K).
! spres Extracted pressure. (Pascals)
! shght Extracted height (meters)
! sqv Extracted water vapor mixing ratio (kg/kg).
! selev extracted surface elevation (m)
! nlevs Number of above-ground sounding levels.
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
! Arguments -- location data
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER, INTENT(IN) :: nx,ny,nz ! Dimensions of ARPS grids.
INTEGER, INTENT(IN) :: nsnd
REAL, INTENT(IN) :: zp(nx,ny,nz) ! z coordinate of grid points in
! physical space (m)
INTEGER, INTENT(IN) :: ipt(nsnd) ! i index to the west of desired
! location
INTEGER, INTENT(IN) :: jpt(nsnd) ! j index to the south of desired
! location
INTEGER, INTENT(IN) :: is_good(nsnd) ! 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)
!-----------------------------------------------------------------------
!
! Arguments -- Extracted sounding variables
!
!-----------------------------------------------------------------------
!
REAL, INTENT(OUT) :: su(nz,nsnd)
REAL, INTENT(OUT) :: sv(nz,nsnd)
REAL, INTENT(OUT) :: sw(nz,nsnd)
REAL, INTENT(OUT) :: stheta(nz,nsnd)
REAL, INTENT(OUT) :: sqv(nz,nsnd)
REAL, INTENT(OUT) :: spres(nz,nsnd)
REAL, INTENT(OUT) :: shght(nz,nsnd)
REAL, INTENT(OUT) :: selev(nsnd)
INTEGER, INTENT(OUT) :: nlevs
!
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: k,isnd
REAL :: w2,w3
REAL :: t1,t2,t3,hmid,tmid,qvsat,rh
!
!-----------------------------------------------------------------------
!
! 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
DO isnd=1,nsnd
if ( is_good(isnd) .eq. 0 ) cycle
DO k=2,nz-1
shght (k,isnd)= zp(ipt(isnd),jpt(isnd),k)
stheta(k,isnd)= ptbar(ipt(isnd),jpt(isnd),k)+ptprt(ipt(isnd),jpt(isnd),k)
spres (k,isnd)= pbar(ipt(isnd),jpt(isnd),k)+ pprt(ipt(isnd),jpt(isnd),k)
sqv (k,isnd)= qvbar(ipt(isnd),jpt(isnd),k)+qvprt(ipt(isnd),jpt(isnd),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=2,nz-2
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 isnd=1,nsnd
if ( is_good(isnd) .eq. 0 ) cycle
DO k=2,nz-1
su(k,isnd) = 0.5*(ubar(ipt(isnd),jpt(isnd),k)+ubar(ipt(isnd)+1,jpt(isnd),k))+ &
0.5*(uprt(ipt(isnd),jpt(isnd),k)+uprt(ipt(isnd)+1,jpt(isnd),k))
sv(k,isnd) = 0.5*(vbar(ipt(isnd),jpt(isnd),k)+vbar(ipt(isnd),jpt(isnd)+1,k))+ &
0.5*(vprt(ipt(isnd),jpt(isnd),k)+vprt(ipt(isnd),jpt(isnd)+1,k))
sw(k,isnd) = 0.5*(wbar(ipt(isnd),jpt(isnd),k)+wbar(ipt(isnd),jpt(isnd),k+1))+ &
0.5*(wprt(ipt(isnd),jpt(isnd),k)+wprt(ipt(isnd),jpt(isnd),k+1))
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 coextract