PROGRAM arpstern,13
!##################################################################
!##################################################################
!###### ######
!###### PROGRAM ARPSTERN ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! A stand alone program created for initializing the ARPS model with
! a terrain data base obtained from NCAR. This program
! will set up the horizontal grid according to the latitude and
! longitude of a point provided by the user for the domain of the
! model run. There is a choice of three different map projections
! at the time this code was prepared. The terrain elevation data
! consists of 30 second, 5 minute and 1 degree resolutions obtained
! from NCAR data services. The structure of this program is as
! follows:
!
!
! 1. Determine the model domain and set up the analysis grid.
! (subroutine setgrid)
!
! 2. Read the user specified data sets and fill the 1 degree blocks
! with the highest resolution data avaliable. (subroutine getter
! and readter)
!
! 3. Apply a filter (Barnes, or other scheme) to the data.
! (subroutine barnes)
!
! 4. Determine and plot the shape/response functions for the given
! filter. (subroutine response)
!
! 5. Calculate and plot the RMS (root mean square) difference
! fields. (subroutine rmsdif)
!
! Note: ARPSTERN is the driver for the ARPS terrain
! initialization.
!
! To compile and link the program:
!
! make arpstern
!
! To execute the program:
!
! arpstern < arpstern.input >! arpstern.output
!
! To create a tar file of all relevant files (except for terrain
! data):
!
! make arpster.tar
!
!-----------------------------------------------------------------------
!
! AUTHOR: Dan Weber
! 1/12/94
! Modification history:
!
! M. Xue 2/7/94.
! User specification of terrain data directory. More error
! checking when opening/reading files.
!
! D. Weber 12/15/94.
! Included pre-set values of input data and print statements
! of variables and parameters.
!
! K. Brewster 1/30/95
! Write file before doing any plotting.
!
! D. Weber 9/8/95.
! Revised the definition of tema and temb in GETTER to include the
! absolute value. This allows the proper setting of the
! DO LOOP limits for use in analyzing southern hemispheric
! data.
!
! 1999/07/16 (Gene Bassett)
! Added patches for Great Lakes and other areas.
!
!
! 2000/11/28 (D. Weber)
! Added dynamic allocation for expanded f90 capabilites. This
! involved setting the nx and ny in the arpstern.input file
! and computing the ntx and nty for the job. Ntx and nty are
! no longer needed as well as terrain.inc.
!
!-----------------------------------------------------------------------
!
! Input read in from the input file:
!
!
! INPUT:
!
!
! COMMON BLOCK VARIABLES:
!
! nx number of grid points in the x direction.
! ny number of grid points in the y direction.
!
! analtype Type of analysis for use in filtering the data
! analtype = 1 Barnes scheme
! = 2 others (not yet implemented)
!
! mapproj Type of map projection used to setup the analysis grid.
! mapproj = 1 Polar Stereographic projection
! = 2 Mercator projection
! = 3 Lambert
!
! itertype Type of terrain data desired,
! itertype = 30 second
! = 5 minute
! = 1 degree
! rmsopt RMS difference control parameter,
! rmsopt = 0 rms difference and background field will
! not be calculated for a 1 pass analysis.
! = 1 rms difference and background field will
! be calculated for a 1 pass analysis.
! NOTE: for multiple pass scheme, the background field
! must be calculated!!!!!!
!
! comtype Computer type the arpstern11.f program will be run on
! comtype = 1 for IBM and other machine capable of
! declaring 2 byte integers
! = 4 for CRAY machines and other machines
! in which the minimum integer declaration
! is integer *8
!
! knot Influence factor for the barnes scheme
! gamma Shape factor for the multiple pass Barnes analysis
! ipass Number of passes for the Barnes routine
! trulat Latitude at which the map projection is true (degree N)
! trulon Longitude at which the map projection is true
! (degree E)
! sclfct Map scale factor (1/1000000=1000000=sclfct)
! dx Analysis grid spacing in the e-w direction (m)
! dy Analysis grid spacing in the n-s direction (m)
! ctrlat Latitude of the center of the analysis grid (deg. N)
! ctrlon Longitude of the center of the analysis grid (deg. E)
! tol Tolerance level for the weight function
! wdn Number of grid points representing the wavelength
! associated
! with the rdnot initial response.
! rdnot First pass or initial response for the wnd length wave.
!
! LOCAL VARIABLES:
!
! glab Graphic label for the terrain plots
! nlat Latitude of the northern most extent of the analysis
! grid (degrees n)
! slat Latitude of the southern most extent of the analysis
! grid (degrees n)
! elon Longitude of the eastern edge of the analysis grid
! (degrees E)
! wlon Longitude of the western edge of the analysis grid
! (degrees E)
!
! nbufx Number of additional terrain data points outside the
! analysis domain used on the calculation of the
! boundaries. East-west
! nbufy Number of additional terrain data points outside the
! analysis domain used on the calculation of the
! boundaries. North-south
! ntx Number of terrain data points in the east-west
! direction
! nty Number of terrain data points in the north-south
! direction
! pi The number pi.
! temak Working knot used by the program (user specified or
! calc. in setgrid) (m**2)
!
!
! h Initial data array (m)
!
!
! OUTPUT:
!
! hback Analyzed Background field at data points (m)
! hterain Analyzed terrain data array written out to
! arpstern.dat. (m)
!
!
! TEMPORARY ARRAYS:
!
!
! tem1,tem2
! tem3,tem4
!
! tem11,tem12,tem13,tem14,tem15,tem16
!
! tem21,tem22,tem23,tem24,tem25,tem26,tem27,tem28
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! INPUT:
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! COMMON BLOCK VARIABLES:
!
!-----------------------------------------------------------------------
!
INTEGER :: nx ! number of grid points in the
! x-direction
INTEGER :: ny ! number of grid points in the
! y-direction
INTEGER :: ndx ! number of data points in the
! x-direction for the desired terrain
! data
INTEGER :: ndy ! number of grid points in the
! y-direction for the desired terrain
! data
INTEGER :: analtype ! Type of analysis technique
INTEGER :: mapproj ! Type of map projection
INTEGER :: itertype ! Type of terrain data
INTEGER :: rmsopt ! RMS difference control parameter
INTEGER :: comtype ! Computer type, =1 for ibm, =4 for
! cray
CHARACTER (LEN=80 ) :: tdatadir ! Directory containing the terrain
! data
CHARACTER (LEN=80 ) :: terndir ! Directory containing the ascii
! terrain data
INTEGER :: ltdir ! Length of none-blank part of
! tdatadir.
REAL :: knot ! Influence parameter for the barnes
! scheme
REAL :: gamma ! Shape factor for the multiple pass
! barnes.
INTEGER :: ipass ! Number of passes for the Barnes
! routine
REAL :: trulat (2) ! Latitude at which the map projection
! is true (degees N)
REAL :: trulon ! Longitude at which the map
! projection is true (degrees E)
REAL :: sclfct ! Map scale factor (m)
REAL :: dx ! Analysis grid spacing in the e-w dir
REAL :: dy ! Analysis grid spacing in the n-s dir
REAL :: ctrlat ! Latitude of the center analysis grid
REAL :: ctrlon ! Longitude of the center of analysis
! grid
REAL :: tol ! Tolerance level for the weight
! function
REAL :: wdn ! Wavelength in terms of gridpoints
! for the rdnot response
REAL :: rdnot ! Initial response to wdn length wave.
!
!-----------------------------------------------------------------------
!
! LOCAL VARIABLES:
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=65) :: glab ! Graphic label for the terrain plots
REAL :: nlat ! Latitude of the northern most
! extent of the analysis grid
! (degrees N)
REAL :: slat ! Latitude of the southern most
! extent of the analysis grid
! (degrees N)
REAL :: elon ! Longitude of the eastern edge of
! analysis grid (degrees E)
REAL :: wlon ! Longitude of the western edge of
! analysis grid (degrees E)
INTEGER :: nbufx ! Number of buffer points in e-w
INTEGER :: nbufy ! Number of buffer points in n-s
INTEGER :: ntx ! Number of terrain grid points e-w
INTEGER :: nty ! Number of terrain grid points n-s
REAL :: pi ! The number pi
REAL :: temak ! Working knot (m**2)
REAL, ALLOCATABLE :: h (:,:) ! Initial terrain data array
! for use in the analysis (m)
!
!-----------------------------------------------------------------------
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: hback(:,:) ! Analyzed background field at data
! points (m)
REAL, ALLOCATABLE :: hterain(:,:) ! Analyzed terrain data array at
! analysis points for use in model.(m)
!
!-----------------------------------------------------------------------
!
! TEMPORARY ARRAYS:
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: tem1 (:)
REAL, ALLOCATABLE :: tem2 (:)
INTEGER*2, ALLOCATABLE :: tem3(:,:) ! Array for reading 1km, 5min, 1deg
! data
INTEGER*2, ALLOCATABLE :: tem4(:,:) ! Array used to store terrain data
REAL, ALLOCATABLE :: tem11(:,:)
REAL, ALLOCATABLE :: tem12(:,:)
REAL, ALLOCATABLE :: tem13(:,:)
REAL, ALLOCATABLE :: tem14(:,:)
INTEGER, ALLOCATABLE :: tem15(:,:)
INTEGER, ALLOCATABLE :: tem16(:,:)
REAL, ALLOCATABLE :: tem21(:,:)
REAL, ALLOCATABLE :: tem22(:,:)
REAL, ALLOCATABLE :: tem23(:,:)
REAL, ALLOCATABLE :: tem24(:,:)
INTEGER, ALLOCATABLE :: tem25(:,:)
INTEGER, ALLOCATABLE :: tem26(:,:)
REAL, ALLOCATABLE :: tem27(:,:)
REAL, ALLOCATABLE :: tem28(:,:)
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
REAL :: tema,temb,temc,temd,trulat1,trulat2,cint
INTEGER :: i,j,ibgn,iend,jbgn,jend
INTEGER :: ierr,istatus
INTEGER :: idummy,imin,jmin,imax,jmax
REAL :: rdummy,hmin,hmax
REAL :: lat, lon, terdx,terdy
REAL :: fix_wisc, fix_bahamas
!
!-----------------------------------------------------------------------
!
! Defining COMMON Blocks:
!
!-----------------------------------------------------------------------
COMMON /terrainc/ analtype,mapproj,itertype,rmsopt,comtype, &
tdatadir,ltdir
COMMON /barnesc/ knot,gamma,ipass
COMMON /projectc/ trulat,trulon,sclfct
COMMON /gridc/ dx,dy,ctrlat,ctrlon,tol,wdn,rdnot
!-----------------------------------------------------------------------
!
! Defining namelists:
!
!-----------------------------------------------------------------------
NAMELIST /terraind/analtype,mapproj,itertype,rmsopt,comtype, &
tdatadir,terndir
NAMELIST /grid/ nx,ny
NAMELIST /barnesd/ knot,gamma,ipass
NAMELIST /mapprojd/trulat1,trulat2,trulon,sclfct
NAMELIST /gridd/ dx,dy,ctrlat,ctrlon,tol,wdn,rdnot
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Set the Namelist variables before they are read in.
! This is done for the case of a read error in the namelist input
! file. Namelist variables will be set to values provided in the
! example in section 9.3 of the ARPS User's Guide.
!
nx = 67
ny = 67
analtype = 1
mapproj = 1
itertype = 1
rmsopt = 1
comtype = 1
tdatadir = 'arpstern.data'
terndir = 'arpstern.data'
knot = 0.0
gamma = 1.00
ipass = 4
trulat1 = 40.70
trulat2 = 0.0
trulon = 248.00
sclfct = 1.0
dx = 1000.0
dy = 1000.0
ctrlat = 40.70
ctrlon = 248.00
tol = 1.0E-7
wdn = 2.00
rdnot = 0.006
!
! Read and Print of Namelist variables and parameters.
!
!-----------------------------------------------------------------------
!
READ(5,grid)
WRITE(6,'(/5x,a,i4)') 'nx was: ',nx
WRITE(6,'(/5x,a,i4)') 'ny was: ',ny
READ(5,terraind)
WRITE(6,'(/5x,a,i4)') 'Analysis type was: ',analtype
WRITE(6,'(/5x,a,i4)') 'Projection type was: ',mapproj
WRITE(6,'(/5x,a,i4)') 'Terrain data type was: ',itertype
IF(itertype == 1)THEN
ndx = 120
ndy = 120
ELSEIF(itertype == 2)THEN
ndx = 12
ndy = 12
ELSEIF(itertype == 3)THEN
ndx = 2
ndy = 2
END IF
WRITE(6,'(/5x,a,i6)')'number of x data points per block was:',ndx
WRITE(6,'(/5x,a,i6)')'number of y data points per block was:',ndy
WRITE(6,'(/5x,a,i4)') 'Rmsopt was: ',rmsopt
WRITE(6,'(/5x,a,i4)') 'Computer type was: ',comtype
WRITE(6,'(/5x,a,a)') 'The source data directory was: ',tdatadir
WRITE(6,'(/5x,a,a)') 'The analyzed data directory was:',terndir
READ(5,barnesd)
WRITE(6,'(/5x,a,f6.3)') 'knot was: ',knot
WRITE(6,'(/5x,a,f6.3)') 'Gamma was: ',gamma
WRITE(6,'(/5x,a,i4)') 'Ipass was: ',ipass
READ(5,mapprojd)
WRITE(6,'(/5x,a,f8.2)') 'trulat1 was: ',trulat1
WRITE(6,'(/5x,a,f8.2)') 'trulat2 was: ',trulat2
WRITE(6,'(/5x,a,f8.2)') 'trulon was: ',trulon
WRITE(6,'(/5x,a,f12.1)') 'The scale factor was:',sclfct
READ(5,gridd)
WRITE(6,'(/5x,a,f9.2)') 'The x grid spacing was: ',dx
WRITE(6,'(/5x,a,f9.2)') 'The y grid spacing was: ',dy
WRITE(6,'(/5x,a,f8.2)') 'The center latitude was: ',ctrlat
WRITE(6,'(/5x,a,f8.2)') 'The center longitude was:',ctrlon
WRITE(6,'(/5x,a,f12.9)') 'The tolerance level was: ',tol
WRITE(6,'(/5x,a,f6.3)') 'The base line wavelength was:',wdn
WRITE(6,'(/5x,a,f6.3)') 'The base line response was:',rdnot
! allocate some arrays associated with nx and ny, ndx and ndy.
ALLOCATE(hterain(nx,ny),STAT=istatus)
hterain = 0
ALLOCATE(tem1(nx),STAT=istatus)
tem1 = 0
ALLOCATE(tem2(ny),STAT=istatus)
tem2 = 0
ALLOCATE(tem3(ndx,ndy),STAT=istatus)
tem3 = 0
ALLOCATE(tem4(ndx,ndy),STAT=istatus)
tem4 = 0
ALLOCATE(tem11(nx,ny),STAT=istatus)
tem11 = 0
ALLOCATE(tem12(nx,ny),STAT=istatus)
tem12 = 0
ALLOCATE(tem13(nx,ny),STAT=istatus)
tem13 = 0
ALLOCATE(tem14(nx,ny),STAT=istatus)
tem14 = 0
ALLOCATE(tem15(nx,ny),STAT=istatus)
tem15 = 0
ALLOCATE(tem16(nx,ny),STAT=istatus)
tem16 = 0
! initial allocation complete....
!
!-----------------------------------------------------------------------
!
! Test to see if the comtype is suitable....
!
!-----------------------------------------------------------------------
!
IF(comtype /= 1.AND.comtype /= 4)THEN
PRINT *,'comtype in terrain.input is not set to IBM or CRAY ', &
'type. Please reset to 1 for IBM or 4 for CRAY'
STOP
END IF
!
!-----------------------------------------------------------------------
!
! setting ctrlon in terms of degrees east
!
!-----------------------------------------------------------------------
!
IF(ctrlon < 0)THEN
ctrlon=360+ctrlon
END IF
!
!-----------------------------------------------------------------------
!
! test to see if reset value is within longitudinal limits..
!
!-----------------------------------------------------------------------
!
IF(ctrlon < 0.OR.ctrlon > 360)THEN
PRINT *,'ctrlon is incorrectly set in terrain.input, please ', &
'reset, use degrees east'
STOP
END IF
!
!-----------------------------------------------------------------------
!
! test to see if the ctrlat is within limits...
!
!-----------------------------------------------------------------------
!
IF(ctrlat < -90.OR.ctrlat > 90)THEN
PRINT *,'ctrlat in terrain.input is too large or too small. ', &
'must be set between -90. and +90. degree north.'
STOP
END IF
!
!-----------------------------------------------------------------------
!
! setting trulon in terms of degrees east
!
!-----------------------------------------------------------------------
!
IF(trulon < 0)THEN
trulon=360+trulon
END IF
!
!-----------------------------------------------------------------------
!
! test to see if reset value is within longitudinal limits..
!
!-----------------------------------------------------------------------
!
IF(trulon < 0.OR.trulon > 360)THEN
PRINT *,'trulon is incorrectly set in terrain.input, please ', &
'reset, use degrees east'
STOP
END IF
!
!-----------------------------------------------------------------------
!
! test to see if trulat is within limits...
!
!-----------------------------------------------------------------------
!
IF(trulat (1) < -90.OR.trulat (1) > 90)THEN
PRINT *,'trulat (1) in terrain.input is too large or too ', &
'small. must be set between -90 and +90 degree north.'
STOP
END IF
IF(trulat (2) < -90.OR.trulat (2) > 90)THEN
PRINT *,'trulat (2) in terrain.input is too large or too ', &
'small. must be set between -90 and +90 degree north.'
STOP
END IF
ltdir=LEN(tdatadir)
CALL strlnth
( tdatadir , ltdir)
IF( ltdir == 0 ) THEN
tdatadir = '.'
ltdir =1
END IF
IF( tdatadir(ltdir:ltdir) /= '/') THEN
ltdir=ltdir+1
tdatadir(ltdir:ltdir)='/'
END IF
PRINT*,tdatadir(ltdir:ltdir)
pi=4.0*ATAN(1.0)
trulat (1)=trulat1
trulat (2)=trulat2
temak=knot
CALL opngks
!
!-----------------------------------------------------------------------
!
! Set the base grid parameters.
!
!-----------------------------------------------------------------------
!
CALL setgrid(nx,ny,ndx,ndy,pi, &
ntx,nty,nbufx,nbufy,slat,nlat,wlon,elon,temak, &
tema,temb,temc,temd,terdx,terdy, &
tem1,tem2,tem11,tem12,tem13,tem14)
!
!-----------------------------------------------------------------------
!
! Allocate the space for the variables used in the analysis. Can't
! do this until the grid parameters are computed (above).
!
!-----------------------------------------------------------------------
!
ALLOCATE( h = 0
ALLOCATE(hback(ntx,nty),STAT=istatus)
hback = 0
ALLOCATE(tem21(ntx,nty),STAT=istatus)
tem21 = 0
ALLOCATE(tem22(ntx,nty),STAT=istatus)
tem22 = 0
ALLOCATE(tem23(ntx,nty),STAT=istatus)
tem23 = 0
ALLOCATE(tem24(ntx,nty),STAT=istatus)
tem24 = 0
ALLOCATE(tem25(ntx,nty),STAT=istatus)
tem25 = 0
ALLOCATE(tem26(ntx,nty),STAT=istatus)
tem26 = 0
ALLOCATE(tem27(ntx,nty),STAT=istatus)
tem27 = 0
ALLOCATE(tem28(ntx,nty),STAT=istatus)
tem28 = 0
!
!-----------------------------------------------------------------------
!
! Set the grid arrays for the model and the extended buffer area.
!
!-----------------------------------------------------------------------
!
call compgrid(nx,ny,ndx,ndy,pi, &
ntx,nty,nbufx,nbufy,slat,nlat,wlon,elon,temak, &
tema,temb,temc,temd,terdx,terdy, &
tem1,tem2,tem11,tem12,tem13,tem14,tem15,tem16, &
tem21,tem22,tem23,tem24,tem25,tem26)
!
!-----------------------------------------------------------------------
!
! Call the GETTER subroutine to read the terrain data files and
! fill the initial terrain array (h).
!
!-----------------------------------------------------------------------
!
CALL getter(ntx,nty,ndx,ndy,slat,nlat,wlon,elon, &
h, &
tem3,tem4)
print *,'initial terrain field'
WRITE (*,'(10f7.1)') ((h(i,j),i=1,10),j=1,10)
!
!-----------------------------------------------------------------------
!
! Fill in Great Lakes. The original data set corresponds to the
! bottom of the lakes.
!
!-----------------------------------------------------------------------
!
CALL fix_lake_eliv(h,tem23,tem24,ntx,nty)
!
!-----------------------------------------------------------------------
!
! Correct for bad terrain point in Wisc. & Bermuda.
!
!-----------------------------------------------------------------------
!
DO j=1,nty
DO i=1,ntx
lat = tem23(i,j)
lon = tem24(i,j) - 360.0
fix_wisc = 290.0
IF (lat > 43.0 .AND. lat < 44.0 .AND. lon > -90.5 .AND. &
lon < -89.0 .AND. h(i,j) > 500.0) THEN
PRINT 9901, 'ARPSTERN: Lowering Wisc: i,j,0,lat,lon,h: ', &
i,j,0,lat,lon,h(i,j),fix_wisc
h(i,j) = fix_wisc
END IF
fix_bahamas = 0.0
IF (lat >= 22.5 .AND. lat <= 24.0 .AND. lon >= -76.5 .AND. &
lon <= -75.0 .AND. ABS( PRINT 9901, 'ARPSTERN: Raising Bahamas: i,j,0,lat,lon,h: ', &
i,j,0,lat,lon,h(i,j),fix_bahamas
h(i,j) = fix_bahamas
END IF
END DO
END DO
9901 FORMAT (a,3I4,2F8.2,2F8.0)
!
!-----------------------------------------------------------------------
!
! This section creates an analytical function for testing purposes
!
!-----------------------------------------------------------------------
!
! DO 2 i=1,ntx
!
! DO 2 j=1,nty
! h(i,j)=20.*SIN(2.*pi*real((i-1.)/12.)) ! the denominator is the
! horizontal wavelength.....
! compute max and min of h...
hmax = -1000000.
hmin = 1000000.
DO j=1,nty
DO i=1,ntx
IF( imax = i
jmax = j
hmax = h(i,j)
END IF
IF( hmin = h(i,j)
imin = i
jmin = j
END IF
END DO
END DO
print *,'hmax is =',hmax,imax,jmax
print *,'hmin is =',hmin,imin,jmin
!
!-----------------------------------------------------------------------
!
! Perform the desired analysis on the data base h.
!
!-----------------------------------------------------------------------
IF (analtype == 1) THEN ! perform a Barnes analysis on the data.
DO i=1,ipass
CALL barnes(i,1,nx,1,ny,nx,ny,ntx,nty,nbufx,nbufy, &
temak,hback,h, &
hterain, &
tem11,tem12,tem21,tem22,tem15,tem16,tem13,tem14)
IF(ipass > 1.OR.rmsopt == 1) THEN ! perform a barnes analysis
! at the data points
ibgn=1+i*nbufx
iend=ntx-i*nbufx
jbgn=1+i*nbufy
jend=nty-i*nbufy
CALL barnes(i,ibgn,iend,jbgn,jend,ntx,nty,ntx,nty, &
nbufx,nbufy,temak,hback,h, &
hback, &
tem21,tem22,tem21,tem22,tem25,tem26,tem27,tem28)
END IF
END DO
!
!-----------------------------------------------------------------------
!
! Write the analyzed terrain data (model grid data) into file
! arpstern.dat.
!
!-----------------------------------------------------------------------
!
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile('arpstern.dat', '-F f77 -N ieee', ierr)
OPEN(11,FILE='arpstern.dat',FORM='unformatted', &
STATUS='unknown')
WRITE(11) nx,ny
idummy = 0
rdummy = 0.0
WRITE(11) analtype,mapproj,itertype,ipass,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy
WRITE(11) dx ,dy ,ctrlat,ctrlon,knot , &
gamma ,trulat1,trulat2,trulon,sclfct, &
tol ,wdn ,rdnot ,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy
WRITE(11) hterain
CLOSE(11)
!
!-----------------------------------------------------------------------
!
! Plotting the results.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Plot the analysis according to the user defined map proj.
!
!-----------------------------------------------------------------------
!
WRITE (glab,65) ipass,gamma
65 FORMAT('Barnes Analysis (m) ','Pass=',i2,' Gamma=',f4.2)
cint=0.
CALL plot(nx,ny,nx,mapproj,trulat,trulon,cint, &
tema,temb,temc,temd,glab,hterain)
!
!-----------------------------------------------------------------------
!
! Plot the initial data field using a cylindrical equidistant proj.
!
!-----------------------------------------------------------------------
!
PRINT *,'initial field, h'
WRITE (*,'(10f7.1)') ((h(i,j),i=1,10),j=1,10)
WRITE (glab,55) itertype
55 FORMAT('Initial Terrain Field (m) ',' Itertype=',i2)
cint=0.
CALL plot(ntx,nty,ntx,8,trulat,trulon,cint, &
tem23(1,nty),tem24(1,nty),tem23(ntx,1), &
tem24(ntx,1),glab,h)
!
!-----------------------------------------------------------------------
!
! Print the maximum initial data set value.
!
!-----------------------------------------------------------------------
!
tema=-1000000.
temb= 1000000.
DO j=1,nty
DO i=1,ntx
tema=MAX(tema,h(i,j))
temb=MIN(temb,h(i,j))
END DO
END DO
PRINT *,'the maximum terrain initial value is',tema
PRINT *,'the minimum terrain initial value is',temb
!
!-----------------------------------------------------------------------
!
! Print the maximum analysis value.
!
!-----------------------------------------------------------------------
!
hmax = -1000000.
hmin = 1000000.
DO j=1,ny
DO i=1,nx
IF(hterain(i,j) > hmax)THEN
imax = i
jmax = j
hmax = hterain(i,j)
END IF
IF(hterain(i,j) < hmax)THEN
hmin = hterain(i,j)
imin = i
jmin = j
END IF
END DO
END DO
print *,'hterain max is =',hmax,imax,jmax
print *,'hterain min is =',hmin,imin,jmin
tema=-100000.
temb= 100000.
DO j=1,ny
DO i=1,nx
tema=MAX(tema,hterain(i,j))
temb=MIN(temb,hterain(i,j))
END DO
END DO
PRINT *,'the maximum analysis value is',tema
PRINT *,'the minimum analysis value is',temb
!
!-----------------------------------------------------------------------
!
! Print the maximum background value.
!
!-----------------------------------------------------------------------
!
IF(ipass > 1.OR.rmsopt == 1)THEN
tema=-100.
DO j=1,nty
DO i=1,ntx
tema=MAX(tema,hback(i,j))
END DO
END DO
PRINT *,'the maximum background analysis value is',tema
!
!-----------------------------------------------------------------------
!
! Plot the background data field using a cylindrical equidistant
! proj.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Setting the lat/lon bounds for the plot routine.
!
!-----------------------------------------------------------------------
!
ibgn= 1+ipass*nbufx
iend= ntx-ipass*nbufx
jbgn= 1+ipass*nbufy
jend= nty-ipass*nbufy
WRITE (glab,56) ipass,gamma
56 FORMAT(1X,'Analyzed Terrain field (m)', &
' Pass=',i2,' Gamma=',f4.2)
cint=0.
CALL plot((iend-ibgn+1),(jend-jbgn+1),ntx, &
8,trulat,trulon,cint, &
tem23(ibgn,jend),tem24(ibgn,jend), &
tem23(iend,jbgn),tem24(iend,jbgn), &
glab,hback(ibgn,jbgn))
!-----------------------------------------------------------------------
!
! Call the RMS difference subroutine to compare the analysis of
! data field at the data points to the initial data field.
!
!-----------------------------------------------------------------------
PRINT *,'tem23,tem24 before rms'
PRINT *,tem23(ibgn,jend),tem24(ibgn,jend),tem23(iend,jbgn) &
,tem24(iend,jbgn)
WRITE (glab,57) ipass,gamma
57 FORMAT(1X,'RMS Difference at terrain points (m)',' Pass=' &
,i2,' Gamma=',f4.2)
CALL rmsdif(ibgn,iend,jbgn,jend,ntx,nty, &
tem23(ibgn,jend),tem24(ibgn,jend), &
tem23(iend,jbgn), &
tem24(iend,jbgn),glab,h,hback,tem27)
END IF
ELSE IF(analtype == 2) THEN !perform other type of analysis
WRITE (6,'(a)') 'Other type of analysis is implemented yet.'
!
!-----------------------------------------------------------------------
!
! Other analysis methods have not yet been developed and implemented
!
!-----------------------------------------------------------------------
!
END IF
CALL clsgks
STOP
END PROGRAM arpstern
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SETGRID ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
SUBROUTINE setgrid(nx,ny,ndx,ndy,pi, & 1,10
ntx,nty,nbufx,nbufy,slat,nlat,wlon,elon,temak, &
temc,temd,teme,temf, terdx,terdy, &
xg1d,yg1d,xg,yg,glat,glon)
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! To set the grid of the analysis domain for use by the chosen
! data analysis technique.
!
! The structure of this program is as follows:
!
! 1. Read the arpstern.input file for the following required info:
!
! - grid spacing in meters (dx,dy)
!
! - type of map projection (polar stereographic, mercator,
! lambert) and standard latitude(s) and longitude of the these
! projections.
!
! - scale factor of the projection (eg. 1/1000000=1000000)
!
! - central latitude and longitude of the model grid
!
!-----------------------------------------------------------------------
!
! AUTHOR: Dan Weber
! 1/12/94
!
!
! REVISED: Dan Weber
! 6/17/94
!
! Modified initial position to center of analysis grid.
! Code cleanup and added comtype option for computer specific
! record length requirements.
!
! MODIFIED: Dan Weber
! 11/28/2000
!
! Removed computation of grid arrays for f90 allocation purposes.
!
!-----------------------------------------------------------------------
!
!
!
! INPUT :
!
! nx Number of grid points for the analysis grid e-w dir.
! ny Number of grid points for the analysis grid n-s dir.
! ndx Number of data points in a 1 degree block (e-w)
! ndy Number of data points in a 1 degree block (n-s)
! pi The constant Pi
!
!-----------------------------------------------------------------------
!
! COMMON BLOCK VARIABLES:
!
!-----------------------------------------------------------------------
!
! analtype Analysis type (1=barnes,2=other,...)
! mapproj Type of map projection used to setup the analysis grid.
! mapproj = 1 Polar Stereographic projection
! = 2 Mercator projection
! = 3 Lambert
! itertype Type of terrain desired (1=30 second, 2= 5 minute,
! 3= 1 degree)
! rmsopt Option for calc. the rms differences. (0=no,1=yes)
! comtype Computer type the arpstern11.f program will be run on
! comtype = 1 for IBM and other machine capable of
! declaring 2 byte integers
! = 4 for CRAY machines and other machines
! in which the minimum integer declaration
! is integer *8
! knot Shape parameter for the barnes weight function.
! gamma Shape factor for multi pass scheme
! ipass Total number of passes thru the barnes scheme
! trulat Latitude(s) at which the map projection is true
! (degree N)
! trulon Longitude at which the map projection is true
! (degree E)
! sclfct Map scale factor (m)
!
! dx Analysis grid spacing in the x-direction east-west (m)
! dy Analysis grid spacing in the y-direction north-south(m)
! ctrlat Latitude of the center of the analysis grid (deg. N)
! ctrlon Longitude of the center of the analysis grid (deg. E)
! tol Tolerance for the weight function used in determining
! the cut-off radius or radius of influence
! wdn Wavelength in terms of grid points in which rdnot is
! valid
! rdnot Initial Response for wdn length wave
!
! LOCAL VARIABLES:
!
! dincx Analysis x-direction grid spacing normalized by the map
! scale dincx=dx/sclfct
! dincy Analysis y-direction grid spacing normalized by the map
! scale dincy=dy/sclfct
! mterdx Max. distance between e-w data pts. (m)
! terdx Distance in meters between data points in the e-w dir.
! (m)
! terdy Distance in meters between data points in the n-s dir.
! (m)
!
! OUTPUT:
!
! ntx Number of data points needed in the east-west direction
! nty Number of data points needed in the north-south
! direction
! nbufx Number of additional points needed to perform the
! Barnes
! analysis at the edges of the domain. e-w
! nbufy Number of additional points needed to perform the
! Barnes
! analysis at the edges of the domain. n-s
! nlat Latitude of the northern edge of the data area
! (degrees N)
! slat Latitude of the southern edge of the data area
! (degrees N)
! elon Longitude of the eastern edge of the data area
! (degrees E)
! wlon Longitude of the western edge of the data area
! (degrees E)
! temak Working value of knot (user or program specified)
! xg1d Analysis grid points e-w (1-d) in grid units
! (=meters/sclfct)
! yg1d Analysis grid points n-s (1-d) in grid units
! (=meters/sclfct)
! xg Analysis grid points in the e-w direction
! (in grid units)
! yg Analysis grid points in the n-s direction
! (in grid units)
! glat Latitude of the analysis data points
! glon Longitude of the analysis data points
!
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! INPUT:
!
!-----------------------------------------------------------------------
!
INTEGER :: nx ! Number of analysis grid points e-w
INTEGER :: ny ! Number of analysis grid points n-s
INTEGER :: ndx ! Number of data points (e-w)
INTEGER :: ndy ! Number of data points (n-s)
REAL :: pi ! The constant Pi
!
!-----------------------------------------------------------------------
!
! COMMON BLOCK VARIABLES:
!
!-----------------------------------------------------------------------
!
INTEGER :: analtype ! Analysis type
INTEGER :: mapproj ! Type of map projection used in
! analysis
INTEGER :: itertype ! Desired terrain resolution for
! analysis
INTEGER :: rmsopt ! Option for enabling the rms diff.
! routine.
INTEGER :: comtype ! Computer type, =1 for ibm, =4 for
! cray
CHARACTER (LEN=80 ) :: tdatadir ! Directory containing the terrain
! data
INTEGER :: ltdir ! Length of none-blank part of
! tdatadir.
REAL :: knot ! Influence parameter for the barnes
! scheme
REAL :: gamma ! Barnes shape factor for multi pass
! scheme
INTEGER :: ipass ! Total number of passes
REAL :: trulat (2) ! Latitude of true map projection
! (deg N)
REAL :: trulon ! Longitude of true map projection
! (deg E)
REAL :: sclfct ! Map scale factor (m)
REAL :: dx ! Analysis grid spacing e-w direction
! (m)
REAL :: dy ! Analysis grid spacing n-s direction
! (m)
REAL :: ctrlat ! Latitude of the center of the
! analysis grid (deg. n)
REAL :: ctrlon ! Longitude of the center of the
! analysis grid (deg. e)
REAL :: tol ! Tolerance level for the weight
! function
REAL :: wdn ! Wavelength in terms of analysis grid
! points for the rdnot response
REAL :: rdnot ! Initial Response for wdn length wave
!
!-----------------------------------------------------------------------
!
! LOCAL VARIABLES:
!
!-----------------------------------------------------------------------
!
REAL :: dincx ! Analysis x-direction grid
! spacing/map scale
REAL :: dincy ! Analysis y-direction grid
! spacing/map scale
REAL :: mterdx ! Max. distance between e-w data pts
REAL :: terdx ! Distance in grid meters between e-w
! data pts
REAL :: terdy ! Distance in grid meters between n-s
! data pts
!
!-----------------------------------------------------------------------
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
INTEGER :: ntx ! Actual number of data points needed
! for data domain e-w
INTEGER :: nty ! Actual number of data points needed
! for data domain n-s
INTEGER :: nbufx ! Number of buffer points in the x
! direction
INTEGER :: nbufy ! Number of buffer points in the y
! direction
REAL :: nlat ! Max. Latitude of the data area
! (degrees n)
REAL :: slat ! Min. Latitude of the data area
! (degrees n)
REAL :: elon ! Max. Longitude of the data area
! (degrees e)
REAL :: wlon ! Min. Longitude of the data area
! (degrees e)
REAL :: temak ! Working value of knot
REAL :: xg1d(nx) ! Analysis grid points e-w (1-d) in
! grid units
REAL :: yg1d(ny) ! Analysis grid points n-s (1-d) in
! grid units
REAL :: xg(nx,ny) ! 2-D analysis grid points east-west
REAL :: yg(nx,ny) ! 2-D analysis grid points north-south
REAL :: glat(nx,ny) ! Latitude of the analysis data points
REAL :: glon(nx,ny) ! Longitude of the analysis data
! points
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k,l,m,n,itema,itemb
REAL :: tema,temb,temc,temd,teme,temf
!
!-----------------------------------------------------------------------
!
! Defining COMMON blocks:
!
!-----------------------------------------------------------------------
COMMON /terrainc/ analtype,mapproj,itertype,rmsopt,comtype &
,tdatadir,ltdir
COMMON /projectc/ trulat,trulon,sclfct
COMMON /gridc/ dx,dy,ctrlat,ctrlon,tol,wdn,rdnot
COMMON /barnesc/ knot,gamma,ipass
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
sclfct=1.0/sclfct
dincx=dx*sclfct
dincy=dy*sclfct
PRINT *,'dincx,dincy'
PRINT *,dincx,dincy
!
!-----------------------------------------------------------------------
!
! Note IMPORTANT!!!!: dx and dy are in meters...and the grid is
! oriented so that the y-axis runs along a longitude line towards
! the northpole and the x-axis is perpendicular to the y-axis.
! Create the x,y grid in grid meters (multiplied by sclfct), the
! origin is the southwest corner of the analysis domain and is
! translated from the center point specified by the user (lat,lon).
!
!
!-----------------------------------------------------------------------
CALL setmapr(mapproj,sclfct,trulat,trulon)
CALL lltoxy(1,1,ctrlat,ctrlon,tema,temb)
!
!-----------------------------------------------------------------------
!
! translate the center point to the first physical point.
!
!-----------------------------------------------------------------------
!
tema=tema-((nx-3)*dincx)/2
temb=temb-((ny-3)*dincy)/2
!
!-----------------------------------------------------------------------
!
! translating from the first physical point to the first scalar
! point.
!
!-----------------------------------------------------------------------
!
tema=tema-dincx/2.0
temb=temb-dincy/2.0
!
!-----------------------------------------------------------------------
!
! Calculate the rest of the analysis grid points in earth
! meters*sclfct
!
!-----------------------------------------------------------------------
!
DO j=1,ny
yg1d(j)=temb+(j-1)*dincy
END DO
DO i=1,nx
xg1d(i)=tema+(i-1)*dincx
END DO
!
!-----------------------------------------------------------------------
!
! Print a sample of the analysis grid values...
!
!-----------------------------------------------------------------------
!
PRINT *,'xg'
WRITE(*,'(6f12.1)') (xg1d(i),i=1,6)
PRINT *,'yg'
WRITE(*,'(3f12.1)') (yg1d(j),j=1,3)
!-----------------------------------------------------------------------
!
! Need to determine the boundaries of the analysis grid domain in
! terms of latitude and longitude. This will be done by taking the
! southwest in terms of analysis grid values (in meters) and
! calculating the lat/lon.
!
!-----------------------------------------------------------------------
CALL xytoll(nx,ny,xg1d,yg1d,glat,glon)
!
!-----------------------------------------------------------------------
!
! Print a sample of the analysis grid lat/lons...
!
!-----------------------------------------------------------------------
!
PRINT *,'glat'
WRITE (*,'(10f7.2)') ((glat(i,j),i=1,10),j=1,3)
PRINT *,'glon'
WRITE (*,'(10f7.2)') ((glon(i,j),i=1,10),j=1,3)
!
!-----------------------------------------------------------------------
!
! Converting the longitudes obtained from the xytoll program to
! degrees east for use in the data sets.
!
!-----------------------------------------------------------------------
!
! commented out the follwoinig code testing 11/28/00
! DO j=1,ny
! DO i=1,nx
! IF(glon(i,j) < 0.0)glon(i,j)=glon(i,j)+360.
! start of new code
! IF(glon(i,j) >= 360.0)glon(i,j)=glon(i,j)-360.
! end of new code
! END DO
! END DO
! PRINT *,'glon adjusted'
! WRITE (*,'(10f7.2)') ((glon(i,j),i=1,nx),j=1,3)
!-----------------------------------------------------------------------
!
! Calculating the number of data points to be included in the
! buffered areas (nbufx,nbufy). But first we need to determine the
! smallest horizontal spacing of the user requested data set for use
! in calculating nbufx and nbufy.
!
! note: lat and lon are in degrees north and east...
!
!-----------------------------------------------------------------------
nlat=-91.
slat=91.
elon=-900.
wlon=900.
DO i=1,nx
nlat=MAX(glat(i,ny),nlat)
slat=MIN(glat(i,1),slat)
elon=MAX(glon(i,ny),elon)
wlon=MIN(glon(i,1),wlon)
END DO
DO j=1,ny
nlat=MAX(glat(1,j),nlat)
nlat=MAX(glat(nx,j),nlat)
slat=MIN(glat(1,j),slat)
slat=MIN(glat(nx,j),slat)
elon=MAX(glon(nx,j),elon)
wlon=MIN(glon(1,j),wlon)
END DO
!
!-----------------------------------------------------------------------
!
! Test to see if maxlat and minlat are within physical limits....
!
!-----------------------------------------------------------------------
!
IF (nlat > 90.0.OR.slat <= -90.0)THEN
PRINT *,'maxlat is greater than 90. or minlat is less than ', &
'-90. degrees latitude'
PRINT *,'maxlat,minlat',nlat,slat
STOP
END IF
!
!-----------------------------------------------------------------------
!
! Print the bound of the lat/lons of the analysis grid...
!
!-----------------------------------------------------------------------
!
PRINT *,'maxlat,minlat',nlat,slat
PRINT *,'maxlon,minlon',elon,wlon
!
!-----------------------------------------------------------------------
!
! Estimating the minumum data spacing in meters
! using the two furthest northward data points. These values
! will be used to determine the largest number of data points
! needed to determine the barnes analysis data buffer
!
!-----------------------------------------------------------------------
!
CALL lltoxy(1,1,nlat,glon(1,ny),tema,temb)
temb=glon(1,ny)+REAL(1./ndx)
CALL lltoxy(1,1,nlat,temb,temc,temd)
terdx=ABS(tema-temc)/sclfct
!
!-----------------------------------------------------------------------
!
! Calculating the maximum terdx at the southern most point in the domain
!
!-----------------------------------------------------------------------
!
CALL lltoxy(1,1,slat,glon(1,1),tema,temb)
temb=glon(1,1)+REAL(1./ndx)
CALL lltoxy(1,1,slat,temb,temc,temd)
mterdx=ABS(tema-temc)/sclfct
!
!-----------------------------------------------------------------------
!
! Calculating terdy at the northern most point in the domain
!
!-----------------------------------------------------------------------
!
CALL lltoxy(1,1,nlat,glon(1,ny),tema,temb)
tema=nlat+REAL(1./ndy)
CALL lltoxy(1,1,tema,glon(1,ny),temc,temd)
terdy=ABS(temb-temd)/sclfct
!
!-----------------------------------------------------------------------
!
! Pring the terrain data spacings...
!
!-----------------------------------------------------------------------
!
PRINT *,'terdx,terdy,mterdx',terdx,terdy,mterdx
PRINT *,'2*mterdx=',2*mterdx,'2*terdy=',2*terdy
!
!-----------------------------------------------------------------------
!
! Setting the Barnes influence factor and response to the desired
! wdn wavelength set by the user in arpstern.input.
!
!-----------------------------------------------------------------------
!
temb=MAX(dx,dy)
tema=wdn*temb*sclfct
IF (knot /= 0.0)THEN ! temak is knot read from
! arpstern.input
temak=knot
ELSE
temak=-(tema/pi)*(tema/pi)*LOG(rdnot)
END IF
!
!-----------------------------------------------------------------------
!
! Calculate and plot the barnes response function
!
!-----------------------------------------------------------------------
!
temb=tema/wdn
temc=temb/sclfct
CALL response(tema,temb,temc,temak,pi)
!
!-----------------------------------------------------------------------
!
! Calculate nbufx and nbufy...
!
!-----------------------------------------------------------------------
!
tema=SQRT(-LOG(tol)*temak)
nbufx=INT(tema/(terdx*sclfct)+1)
nbufy=INT(tema/(terdy*sclfct)+1)
!
!-----------------------------------------------------------------------
!
! Print buffer sizes, working knot and cutoff radius....
!
!-----------------------------------------------------------------------
!
PRINT *, 'nbufx,nbufy',nbufx,nbufy
PRINT *, 'knot,temak,rc',knot,temak,tema
!-----------------------------------------------------------------------
!
! Set up the limits of the latitude and longitude values for the
! data area.
!
!-----------------------------------------------------------------------
temc=REAL(1./ndy)
temd=REAL(1./ndx)
slat=slat-REAL((ipass+1)*nbufy*temc)
wlon=wlon-REAL((ipass+1)*nbufx*temd)
tema=slat-INT(slat)
temb=wlon-INT(wlon)
tema=INT(tema/temc)
temb=INT(temb/temd)
slat=REAL(INT(slat)+tema*temc)
wlon=REAL(INT(wlon)+temb*temd)
nlat=nlat+REAL((ipass+1)*nbufy*temc)
elon=elon+REAL((ipass+1)*nbufx*temd)
tema=nlat-INT(nlat)
temb=elon-INT(elon)
! tema=REAL(INT(tema/temc)+1.) ! TESTING ORIGINAL CODE
! temb=REAL(INT(temb/temd)+1.) ! TESTING ORIGINAL CODE
tema=INT(tema/temc) ! TESTING NEW CODE
temb=INT(temb/temd) ! TESTING NEW CODE
nlat=REAL(INT(nlat)+tema*temc)
elon=REAL(INT(elon)+temb*temd)
!
!-----------------------------------------------------------------------
!
! Print lat/lon bounds of the initial data field...
!
!-----------------------------------------------------------------------
!
PRINT *,'Buffer adjusted nlat,slat,wlon,elon'
PRINT *,nlat,slat,wlon,elon
!-----------------------------------------------------------------------
!
! Determine ntx, nty the total number of data points needed
! to perform the barnes analysis for the analysis domain.
!
!-----------------------------------------------------------------------
! original code
nty=INT(REAL(nlat-slat)*ndy)+1
ntx=INT(REAL(elon-wlon)*ndx)+1
! testing new code...
nty=INT(REAL(nlat-slat)*ndy)
ntx=INT(REAL(elon-wlon)*ndx)
PRINT *,'ntx,nty',ntx,nty
PRINT *,'the number of terrain points in the x dir is',ntx
PRINT *,'the number of terrain points in the y dir is',nty
RETURN
END SUBROUTINE setgrid
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE COMPGRID ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
SUBROUTINE compgrid(nx,ny,ndx,ndy,pi, & 1,1
ntx,nty,nbufx,nbufy,slat,nlat,wlon,elon,temak, &
temc,temd,teme,temf,terdx,terdy, &
xg1d,yg1d,xg,yg,glat,glon,closi,closj, &
xob,yob,doblat,doblon,closit,closjt)
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! To compute the grid point locations for use by the chosen
! data analysis technique. Added due to f90 changes.
!
!
!-----------------------------------------------------------------------
!
! AUTHOR: Dan Weber
! 11/28/2000
!
!-----------------------------------------------------------------------
!
!
!
! INPUT :
!
! nx Number of grid points for the analysis grid e-w dir.
! ny Number of grid points for the analysis grid n-s dir.
! ntx number of e-w data points
! nty number of n-s data points
! ndx Number of data points in a 1 degree block (e-w)
! ndy Number of data points in a 1 degree block (n-s)
! pi The constant Pi
!
!-----------------------------------------------------------------------
!
! COMMON BLOCK VARIABLES:
!
!-----------------------------------------------------------------------
!
! analtype Analysis type (1=barnes,2=other,...)
! mapproj Type of map projection used to setup the analysis grid.
! mapproj = 1 Polar Stereographic projection
! = 2 Mercator projection
! = 3 Lambert
! itertype Type of terrain desired (1=30 second, 2= 5 minute,
! 3= 1 degree)
! rmsopt Option for calc. the rms differences. (0=no,1=yes)
! comtype Computer type the arpstern11.f program will be run on
! comtype = 1 for IBM and other machine capable of
! declaring 2 byte integers
! = 4 for CRAY machines and other machines
! in which the minimum integer declaration
! is integer *8
! knot Shape parameter for the barnes weight function.
! gamma Shape factor for multi pass scheme
! ipass Total number of passes thru the barnes scheme
! trulat Latitude(s) at which the map projection is true
! (degree N)
! trulon Longitude at which the map projection is true
! (degree E)
! sclfct Map scale factor (m)
!
! dx Analysis grid spacing in the x-direction east-west (m)
! dy Analysis grid spacing in the y-direction north-south(m)
! ctrlat Latitude of the center of the analysis grid (deg. N)
! ctrlon Longitude of the center of the analysis grid (deg. E)
! tol Tolerance for the weight function used in determining
! the cut-off radius or radius of influence
! wdn Wavelength in terms of grid points in which rdnot is
! valid
! rdnot Initial Response for wdn length wave
!
! LOCAL VARIABLES:
!
! dincx Analysis x-direction grid spacing normalized by the map
! scale dincx=dx/sclfct
! dincy Analysis y-direction grid spacing normalized by the map
! scale dincy=dy/sclfct
! mterdx Max. distance between e-w data pts. (m)
! terdx Distance in meters between data points in the e-w dir.
! (m)
! terdy Distance in meters between data points in the n-s dir.
! (m)
!
! OUTPUT:
!
! ntx Number of data points needed in the east-west direction
! nty Number of data points needed in the north-south
! direction
! nbufx Number of additional points needed to perform the
! Barnes
! analysis at the edges of the domain. e-w
! nbufy Number of additional points needed to perform the
! Barnes
! analysis at the edges of the domain. n-s
! nlat Latitude of the northern edge of the data area
! (degrees N)
! slat Latitude of the southern edge of the data area
! (degrees N)
! elon Longitude of the eastern edge of the data area
! (degrees E)
! wlon Longitude of the western edge of the data area
! (degrees E)
! temak Working value of knot (user or program specified)
! xg1d Analysis grid points e-w (1-d) in grid units
! (=meters/sclfct)
! yg1d Analysis grid points n-s (1-d) in grid units
! (=meters/sclfct)
! xg Analysis grid points in the e-w direction
! (in grid units)
! yg Analysis grid points in the n-s direction
! (in grid units)
! glat Latitude of the analysis data points
! glon Longitude of the analysis data points
! closi Array storing the data i index value corresponding to
! the closest analysis grid point
! closj Array storing the data j index value corresponding to
! the closest analysis grid point
! xob Analysis grid value in the east-west direction for the
! data points (terrain data)
! yob Analysis grid value in the north-south direction for
! the data points (terrain data)
! doblat Latitude of the data points (degrees north)
! doblon Longitude of the data points (degrees east)
! closit Array storing the terrain i index value corresponding
! to the closest terrain point
! closjt Array storing the terrain j index value corresponding
! to the closest terrain point
!
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! INPUT:
!
!-----------------------------------------------------------------------
!
INTEGER :: nx ! Number of analysis grid points e-w
INTEGER :: ny ! Number of analysis grid points n-s
INTEGER :: ndx ! Number of data points (e-w)
INTEGER :: ndy ! Number of data points (n-s)
REAL :: pi ! The constant Pi
!
!-----------------------------------------------------------------------
!
! COMMON BLOCK VARIABLES:
!
!-----------------------------------------------------------------------
!
INTEGER :: analtype ! Analysis type
INTEGER :: mapproj ! Type of map projection used in
! analysis
INTEGER :: itertype ! Desired terrain resolution for
! analysis
INTEGER :: rmsopt ! Option for enabling the rms diff.
! routine.
INTEGER :: comtype ! Computer type, =1 for ibm, =4 for
! cray
CHARACTER (LEN=80 ) :: tdatadir ! Directory containing the terrain
! data
INTEGER :: ltdir ! Length of none-blank part of
! tdatadir.
REAL :: knot ! Influence parameter for the barnes
! scheme
REAL :: gamma ! Barnes shape factor for multi pass
! scheme
INTEGER :: ipass ! Total number of passes
REAL :: trulat (2) ! Latitude of true map projection
! (deg N)
REAL :: trulon ! Longitude of true map projection
! (deg E)
REAL :: sclfct ! Map scale factor (m)
REAL :: dx ! Analysis grid spacing e-w direction
! (m)
REAL :: dy ! Analysis grid spacing n-s direction
! (m)
REAL :: ctrlat ! Latitude of the center of the
! analysis grid (deg. n)
REAL :: ctrlon ! Longitude of the center of the
! analysis grid (deg. e)
REAL :: tol ! Tolerance level for the weight
! function
REAL :: wdn ! Wavelength in terms of analysis grid
! points for the rdnot response
REAL :: rdnot ! Initial Response for wdn length wave
!
!-----------------------------------------------------------------------
!
! LOCAL VARIABLES:
!
!-----------------------------------------------------------------------
!
REAL :: terdx ! Distance in grid meters between e-w
! data pts
REAL :: terdy ! Distance in grid meters between n-s
! data pts
!
!-----------------------------------------------------------------------
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
INTEGER :: ntx ! Actual number of data points needed
! for data domain e-w
INTEGER :: nty ! Actual number of data points needed
! for data domain n-s
INTEGER :: nbufx ! Number of buffer points in the x
! direction
INTEGER :: nbufy ! Number of buffer points in the y
! direction
REAL :: nlat ! Max. Latitude of the data area
! (degrees n)
REAL :: slat ! Min. Latitude of the data area
! (degrees n)
REAL :: elon ! Max. Longitude of the data area
! (degrees e)
REAL :: wlon ! Min. Longitude of the data area
! (degrees e)
REAL :: temak ! Working value of knot
REAL :: xg1d(nx) ! Analysis grid points e-w (1-d) in
! grid units
REAL :: yg1d(ny) ! Analysis grid points n-s (1-d) in
! grid units
REAL :: xg(nx,ny) ! 2-D analysis grid points east-west
REAL :: yg(nx,ny) ! 2-D analysis grid points north-south
REAL :: glat(nx,ny) ! Latitude of the analysis data points
REAL :: glon(nx,ny) ! Longitude of the analysis data
! points
INTEGER :: closi(nx,ny) ! Array storing the data i index value
! corresponding to the closest
! analysis pt.
INTEGER :: closj(nx,ny) ! Array storing the data j index value
! corresponding to the closest
! analysis pt.
REAL :: xob(ntx,nty) ! Data set grid point values east-west
REAL :: yob(ntx,nty) ! Data set grid point values
! north-south
REAL :: doblat(ntx,nty) ! Latitude of the data points
! (degrees north)
REAL :: doblon(ntx,nty) ! Longitude of the data points
! (degrees east)
INTEGER :: closit(ntx,nty) ! array storing the data i index value
! corresponding to the closest data
! point
INTEGER :: closjt(ntx,nty) ! array storing the data j index value
! corresponding to the closest data
! point
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k,l,m,n,itema,itemb
REAL :: tema,temb,temc,temd,teme,temf
!
!-----------------------------------------------------------------------
!
! Defining COMMON blocks:
!
!-----------------------------------------------------------------------
COMMON /terrainc/ analtype,mapproj,itertype,rmsopt,comtype &
,tdatadir,ltdir
COMMON /projectc/ trulat,trulon,sclfct
COMMON /gridc/ dx,dy,ctrlat,ctrlon,tol,wdn,rdnot
COMMON /barnesc/ knot,gamma,ipass
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
! Print lat/lon bounds of the initial data field...
!
!-----------------------------------------------------------------------
!
temc=REAL(1./ndy)
temd=REAL(1./ndx)
PRINT *,'in compgrid nlat,slat,wlon,elon'
PRINT *,nlat,slat,wlon,elon
DO j=1,nty
DO i=1,ntx
doblat(i,j)=slat+REAL((j-1.)*temc)
doblon(i,j)=wlon+REAL((i-1.)*temd)
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Print a sample of the initial data field lat/lon values...
!
!-----------------------------------------------------------------------
!
PRINT *,'doblon'
WRITE (*,'(8f9.3)') ((doblon(i,j),i=1,8),j=1,3)
PRINT *,'doblat'
WRITE (*,'(8f9.3)') ((doblat(i,j),i=1,8),j=1,3)
!-----------------------------------------------------------------------
!
! Calculate the grid point values for the data
!
!-----------------------------------------------------------------------
CALL lltoxy(ntx,nty,doblat,doblon,xob,yob)
!
!-----------------------------------------------------------------------
!
! Print a sample of the initial data field grid values...
!
!-----------------------------------------------------------------------
!
PRINT *,'xob'
WRITE (*,'(6f12.1)') ((xob(i,j),i=1,7),j=1,3)
PRINT *,'yob'
WRITE (*,'(6f12.1)') ((yob(i,j),i=1,7),j=1,3)
!
!-----------------------------------------------------------------------
!
! Setting the corners of the analysis grid in degrees lat/lon for the
! purpose of plotting later.( glat and glon will be overwritten in
! subroutine barnes)
!
!-----------------------------------------------------------------------
!
temc=glat(1,ny)
temd=glon(1,ny)
teme=glat(nx,1)
temf=glon(nx,1)
!
!-----------------------------------------------------------------------
!
! Looking for the data point which is closest to the analysis
! point (1,1). These indicies will be used to find the remaining points
!
!-----------------------------------------------------------------------
!
print *,yg1d(1),xg1d(1)
tema=10000000.
DO l=1,nty/2
DO k=1,ntx/2
temb=SQRT((xg1d(1)-xob(k,l))*(xg1d(1)-xob(k,l))+ &
(yg1d(1)-yob(k,l))*(yg1d(1)-yob(k,l)))
IF(temb < tema)THEN
m=k
n=l
tema=temb
END IF
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Calculate an increment which insures each closest point will be found
!
!-----------------------------------------------------------------------
!
itema=INT(dx/terdx+3)
itemb=INT(dy/terdy+3)
DO j=1,ny
DO i=1,nx
tema=10000000.
DO l=n-itemb,n+itemb
DO k=m-itema,m+itema
temb=SQRT((xg1d(i)-xob(k,l))*(xg1d(i)-xob(k,l))+ &
(yg1d(j)-yob(k,l))*(yg1d(j)-yob(k,l)))
IF (temb < tema)THEN
closi(i,j)=k
closj(i,j)=l
tema=temb
END IF
END DO
END DO
m=closi(i,j)
n=closj(i,j)
END DO
m=closi(1,j)
n=closj(1,j)
END DO
!
!-----------------------------------------------------------------------
!
! Print a sample of closi and closj...
!
!-----------------------------------------------------------------------
!
PRINT *,'closi'
WRITE (*,'(14i5)') ((closi(i,j),i=1,14),j=1,3)
PRINT *,'closj'
WRITE (*,'(14i5)') ((closj(i,j),i=1,14),j=1,3)
!
!-----------------------------------------------------------------------
!
! Set the closit, closjt for the terrain data points analysis.
!
!-----------------------------------------------------------------------
!
DO j=1,nty
DO i=1,ntx
closit(i,j)=i
closjt(i,j)=j
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Expand the 1 dimensional grid arrays into 2-D arrays for use
! in the barnes scheme.
!
!-----------------------------------------------------------------------
!
DO j=1,ny
DO i=1,nx
xg(i,j)=xg1d(i)
yg(i,j)=yg1d(j)
END DO
END DO
RETURN
END SUBROUTINE compgrid
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GETTER ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE getter(ntx,nty,ndx,ndy,slat,nlat,wlon,elon, & 1,3
h, &
nval,readdat)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! To retrieve elevation data from terrain data files and fill (if
! needed) areas which are not covered by the user specified data set
! . This program will read the desired terrain resolution header
! file, if the data is found the program will read the data in and
! load the terrain array. If the initial read does not find the
! data avaliable, the program will read the header file for the next
! best avaliable resolution set. This process will continue until
! data has been found. If no data exists for this particular 1
! degree block (this can happen!!!) the program will stop.
!
! Three terrain resolutions are currently supported by this code:
!
! 1. 30 second resolution data for the continental United
! States(approx. 1km resolution)
! 2. 5 minute resolution data for most of North America and
! some of Europe.
! (approx. 10km resolution)
! 3. 1 degree resolution (30 minute in some locations) for the
! entire globe. (minus a few locations including the south
! pole) (appros. 110km resolution)
!
! The structure of this program is as follows: (eg. 30 second data
! chosen)
! 1. Read 30 second header for a block match. If block is
! found, read the data file and continue to the next block.
! If no 30 second data exists for this block, the 5 minute
! header file will be read to find a block match. If there
! is a match, the data will be read, and expanded to fit the
! size of the 30 second data array. If there is no match
! with the 5 minute header file, the 1 degree header file
! will be read. If there is a match the 1 degree data will
! be read and expanded to fit the 30 second data structure.
! If no data exists for the block in the 1 degree header
! file the program will print a stop message and end.
!
!
!-----------------------------------------------------------------------
!
! AUTHOR: Dan Weber
! 1/12/94
!
! MODIFICATION HISTORY:
!
! 2/16/98, Dan Weber.
! Bug fix for tema and temb for southern hemisphere indexing.
!
! 2000/02/03 Gene Bassett.
! Moved Great Lakes terrain corrections into a subroutine.
!
! 11/29/2000 Dan Weber
! Bug fix to the loading of the first terrain value along a longitude
! circle.
! Bug fix for data sets that span the 0 longitude line.
!
!-----------------------------------------------------------------------
!
!
! INPUT :
!
! ntx Limit of terrain data points for use in the east-west dir.
! nty Limit of terrain data points for use in the north-south dir.
! ndx Number of terrain data points in the (east-west) direction
! for a 1 degree block of data
! ndy Number of terrain data points in the (north-south) direction
! for a 1 degree block of data
! nlat Northern extent of the buffered data set (degrees N)
! slat Southern extent of the buffered data set (degrees N)
! elon Eastern extent of the buffered data set (degrees E)
! wlon Western extent of the buffered data set (degrees E)
!
!
! COMMON BLOCK VARIABLES:
!
! analtype Type of analysis variable requested
! mapproj Type of projection desired
! itertype Type of terrain data requested
! rmsopt Rms error option
! comtype Type of computer this analysis is run on
!
! tdatadir Directory containing the terrain data
! ltdir Length of none-blank part of tdatadir
!
! LOCAL VARIABLES :
!
! id Array for reading header information from the data files
! idim East-west dimension of the readdat array in the readter sub.
! jdim North-south dimension of the readdat array in the readter sub.
! n Record number for use in readter
! nfile2 Character string storing the data file name
! rlength Record length of the data set
!
! OUTPUT:
!
! h Terrain data for use in the Barnes analysis scheme (m)
!
!
! TEMPORARY ARRAYS:
!
! nval Temporary array which contains the terrain data (m)
! readdat Temporary array used to read the data file
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! INPUT:
!
!-----------------------------------------------------------------------
!
INTEGER :: ntx ! Limit of terrain data points e-w dir.
INTEGER :: nty ! Limit of terrain data points n-s dir.
INTEGER :: ndx ! Number of terrain data points in the e-w
INTEGER :: ndy ! Number of terrain data points in the n-s
REAL :: nlat ! Northern extent of the buffered data set (N)
REAL :: slat ! Southern extent of the buffered data set (N)
REAL :: elon ! Eastern extent of the buffered data set (E)
REAL :: wlon ! Western extent of the buffered data set (E)
!
!-----------------------------------------------------------------------
!
! COMMON BLOCK VARIABLES:
!
!-----------------------------------------------------------------------
!
INTEGER :: analtype ! Type of analysis variable requested
INTEGER :: mapproj ! Type of projection desired
INTEGER :: itertype ! =1,2,or 3 for 30 sec,5 min, or 1 deg data
INTEGER :: rmsopt ! Rms error option
INTEGER :: comtype ! Computer type
CHARACTER (LEN=80 ) :: tdatadir ! Directory containing the terrain data
INTEGER :: ltdir ! Length of none-blank part of tdatadir.
!
!-----------------------------------------------------------------------
!
! LOCAL VARIABLES:
!
!-----------------------------------------------------------------------
!
INTEGER :: id(2) ! Data file header array
INTEGER :: idim ! East-west dimension of the readdat array
! in the readter subroutine
INTEGER :: jdim ! North-south dimension of the readdat array
! in the readter subroutine
INTEGER :: n ! Record number for use in readter
CHARACTER (LEN=80) :: nfile2 ! Character storing the data file name
INTEGER :: rlength ! Length of the record in a data set
!
!-----------------------------------------------------------------------
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
REAL :: h(ntx,nty) ! Terrain data for use in the analysis scheme,
! includes the buffered area. (m)
!
!-----------------------------------------------------------------------
!
! TEMPORARY ARRAYS:
!
!-----------------------------------------------------------------------
!
INTEGER*2 nval(ndx,ndy) ! Temporary array which contains the data (m)
INTEGER*2 readdat(ndx,ndy)! Temporary array which used to read the file
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
REAL :: newslat,newnlat,newelon,newwlon,teme,temf
INTEGER :: i,j,k,l,ll,inew,acount
INTEGER :: tema,temb,temc,temd
INTEGER :: ii,jj,kk,kbgn,kend,lbgn,lend,istat
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
! COMMON block definition:
!
!-----------------------------------------------------------------------
!
COMMON /terrainc/ analtype,mapproj,itertype,rmsopt,comtype &
,tdatadir,ltdir
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
PRINT *,'entering getter'
!-----------------------------------------------------------------------
!
! Setting the internal dimensions of the readdat array specific to the
! type of data read.
!
! NOTE: If the desired data resolution is 5 minutes, the 30 second data
! file will not be accessed. If the desired data resolution is 1 degree
! the 5 minute and 30 second data sets will not be accessed.
!
!-----------------------------------------------------------------------
PRINT *,'itertype=',itertype
IF(itertype > 3.OR.itertype < 1)THEN
PRINT *,'itertype is not supported, program will stop'
STOP
END IF
OPEN(12,FILE=tdatadir(1:ltdir)//'dir30sec.hdr', &
FORM='formatted',STATUS='old',IOSTAT=istat)
IF( istat /= 0) THEN
WRITE(6,'(/a,/a/)') &
' Error occured when opening terrain header file ', &
tdatadir(1:ltdir)//'dir30sec.hdr',' Program stopped in GETTER.'
STOP
END IF
OPEN(13,FILE=tdatadir(1:ltdir)//'dir5min.hdr', &
FORM='formatted',STATUS='old',IOSTAT=istat)
IF( istat /= 0) THEN
WRITE(6,'(/a,/a/)') &
' Error occured when opening terrain header file.', &
tdatadir(1:ltdir)//'dir5min.hdr',' Program stopped in GETTER.'
STOP
END IF
OPEN(14,FILE=tdatadir(1:ltdir)//'dir1deg.hdr', &
FORM='formatted',STATUS='old',IOSTAT=istat)
IF( istat /= 0) THEN
WRITE(6,'(/a,/a/)') &
' Error occured when opening terrain header file.', &
tdatadir(1:ltdir)//'dir1deg.hdr',' Program stopped in GETTER.'
STOP
END IF
!
! new code added to fix a southern hemisphere indexing bug...
! (if statements)
IF(slat >= 0.0)THEN ! compute the index for the southern lat.
tema=INT((slat-INT(slat))*ndy)+1
ELSE ! compute the index for the southern lat.
IF( ABS(slat-INT(slat)) < 0.000001)THEN
tema = 1
ELSE
tema=ABS(INT((1.+slat-INT(slat))*ndy))+1
END IF
END IF
IF(nlat >= 0.0)THEN ! compute the index for the northern lat.
temb=INT((nlat-INT(nlat))*ndy)+1
ELSE ! compute the index for the northern lat.
IF(abs(nlat-INT(nlat)) < 0.000001) THEN
temb = 1
ELSE
temb=ABS(INT((1.+nlat-INT(nlat))*ndy))+1
END IF
END IF
! added/modified the following to fix a bug DBW, 11/29/2000
IF(wlon < 0.0)THEN
IF(abs(wlon-INT(wlon)) < 0.000001) THEN
print *,'got here a'
temc = 1
ELSE
temc=ABS(INT((1.+wlon-INT(wlon))*ndx))+1
print *,'got here b'
END IF
ELSE
temc=INT((wlon-INT(wlon))*ndx)+1
END IF
IF(elon < 0.0)THEN
IF(abs(elon-INT(elon)) < 0.000001) THEN
temd = 1
ELSE
temd=ABS(INT((1.+elon-INT(elon))*ndx))+1
END IF
ELSE
temd=INT((elon-INT(elon))*ndx)+1
END IF
PRINT *,'tema,temb,temc,temd',tema,temb,temc,temd
!
!-----------------------------------------------------------------------
!
! Starting the main search and read loop.....
!
!-----------------------------------------------------------------------
!
jj=0
IF(slat < 0.0)THEN ! testing for negative latitude limits...
newslat=slat
teme = slat-INT(slat)
IF(ABS(teme) > 0.000001 )newslat=slat-1.0
ELSE
newslat=slat
END IF
IF(nlat < 0.0)THEN ! testing for negative latitude limits...
newnlat=nlat
teme = nlat-INT(nlat)
IF(ABS(teme) > 0.000001 )newnlat=nlat-1.0
ELSE
newnlat=nlat
END IF
IF(wlon < 0.0)THEN ! testing for negative longitude limits...
newwlon=wlon
teme = wlon-INT(wlon)
IF(ABS(teme) > 0.000001 )newwlon=wlon-1.0
ELSE
newwlon=wlon
END IF
IF(elon < 0.0)THEN ! testing for negative longitude limits...
newelon=elon
teme = elon-INT(elon)
IF(ABS(teme) > 0.000001 )newelon=elon-1.0
ELSE
newelon=elon
END IF
DO j=INT(newslat),INT(newnlat)
jj=jj+1
ii=0
acount = 0
DO i=INT(newwlon),INT(newelon)
ii=ii+1
IF(i >= 360.AND.acount == 0)THEN ! need to rewind the header files...
acount = 1
IF(itertype == 3)THEN ! rewind dir1deg.hdr only...
REWIND(14)
ELSE IF(itertype == 2)THEN ! rewind dir5min and dir1deg.hdrs
REWIND(13)
REWIND(14)
ELSE IF(itertype == 1)THEN ! rewind all header files...
REWIND(12)
REWIND(13)
REWIND(14)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Testing for values of longitude greater than 359.
!
!-----------------------------------------------------------------------
!
inew=i
IF(i >= 360)THEN
inew=inew-360
END IF
! start of bug fix area for data spanning the 0 longitude line...
! DBW 11/29/2000
IF(i < 0)THEN
inew=inew+360
ENDIF
IF(int(wlon) < 0.AND.(i >= 0.AND.i-1 < 0))THEN
IF(itertype.eq.3)then ! rewind dir1deg.hdr only...
rewind(14)
ELSE IF(itertype == 2)then ! rewind dir5min and dir1deg.hdrs
rewind(13)
rewind(14)
ELSE IF(itertype == 1)then ! rewind all header files...
rewind(12)
rewind(13)
rewind(14)
ENDIF
ENDIF
! end of bug fix area DBW 11/29/2000.
PRINT *,'looking for latitude, longitude',j,inew
IF(itertype == 1)THEN ! reading the 30 second header file
1 READ(12,'(3i6)',END=2) id(1),id(2),n
IF(id(1) == j.AND.id(2) == inew)THEN ! testing for a match
nfile2=tdatadir(1:ltdir)//'dir30sec.dat '
PRINT *,'we have a match reading ',nfile2
idim=120
jdim=120
rlength=28800*comtype
CALL readter(n,idim,jdim,ndx,ndy,nfile2,rlength, &
readdat,nval)
GO TO 99 ! Data has been read, add to the h array...
ELSE ! no match yet read header file again...
GO TO 1
END IF
END IF
2 PRINT *,'block not found in dir30sec.hdr or data not desired'
REWIND (12)
IF(itertype == 2.OR.itertype == 1)THEN ! testing for correct itertype
3 READ(13,'(3i6)',END=4) id(1),id(2),n
IF(id(1) == j.AND.id(2) == inew)THEN ! testing for a match (5 minute)
nfile2=tdatadir(1:ltdir)//'dir5min.dat '
PRINT *,'reading ',nfile2,n
idim=12
jdim=12
rlength=288*comtype
CALL readter(n,idim,jdim,ndx,ndy,nfile2,rlength, &
readdat,nval)
GO TO 99
ELSE
GO TO 3
END IF
END IF
4 PRINT *,'block not found in dir5min.hdr or data not desired'
REWIND (13)
!
!-----------------------------------------------------------------------
!
! Reading the dir1deg header file.
!
!-----------------------------------------------------------------------
!
PRINT *,'reading dir1deg.hdr file'
5 READ(14,'(3i6)',END=6) id(1),id(2),n
IF(id(1) == j.AND.id(2) == inew)THEN ! testing for a match...
nfile2=tdatadir(1:ltdir)//'dir1deg.dat'
PRINT *,'reading ',nfile2
idim=2
jdim=2
rlength=8*comtype
CALL readter(n,idim,jdim,ndx,ndy,nfile2,rlength, &
readdat,nval)
print *,nval(1,1),nval(2,1),nval(1,2),nval(2,2)
print *,readdat(1,1),readdat(2,1),readdat(1,2),readdat(2,2)
GO TO 99 ! update the terrain array.....
ELSE
GO TO 5 ! read the header file again......
END IF
6 PRINT *,'block not found in dir1deg.hdr,no data was found ', &
'for this block, the program will stop'
STOP
!
!-----------------------------------------------------------------------
!
! Update the terrain array (h), unpack, and convert to meters.
!
!-----------------------------------------------------------------------
!
99 PRINT *,'adding nval to the h array'
!
!-----------------------------------------------------------------------
!
! Setting up the ranges for adding each nval to the h array....
!
!-----------------------------------------------------------------------
!
IF(jj == 1.AND.j == INT(newnlat))THEN
lbgn=tema
lend=temb
ll=0
ELSE IF(jj == 1)THEN
lbgn=tema
lend=ndy
ll=0
ELSE IF (j == INT(newnlat))THEN
lbgn=1
lend=temb
ll=ndy-tema+1+(jj-2)*ndy
ELSE
lbgn=1
lend=ndy
ll=ndy-tema+1+(jj-2)*ndy
END IF
DO l=lbgn,lend
ll=ll+1
IF(ii == 1.AND.i == INT(newelon))THEN
kbgn=temc
kend=temd
kk=0
ELSE IF(ii == 1)THEN
kbgn=temc
kend=ndx
kk=0
ELSE IF(i == INT(newelon))THEN
kbgn=1
kend=temd
kk=ndx-temc+1+(ii-2)*ndx
ELSE
kbgn=1
kend=ndx
kk=ndx-temc+1+(ii-2)*ndx
END IF
DO k=kbgn,kend
kk=kk+1
h(kk,ll)=REAL(nval(k,l)*20-4000) ! units are in meters...
! print *,'kk ll and h = ',kk,ll,h(kk,ll),nval(k,l),k,l
END DO
END DO
print *,'kbgn,kend,lbgn,lend ',kbgn,kend,lbgn,lend
END DO
END DO
CLOSE(12)
CLOSE(13)
CLOSE(14)
RETURN
END SUBROUTINE getter
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE READTER ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE readter(n,idim,jdim,ndx,ndy,nfile2,rlength, & 3
readdat,nval)
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! To retrieve terrain data from the selected data files
! This program will read either the 30 second, 5 minute, or 1 degree
! terrain data base files in a direct access format.
! (original sources, DMA_ELEV.DAT or elev.dat data files obtain from NCAR)
!
! NOTE: these original data files have been modified to run more
! efficiently (put in direct access format by dir30sec.f, dir5min.f and
! dir1deg.f).
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable definitions
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! n value of the record number
! idim x dimension of the readdat array
! jdim y dimension of the readdat array
! ndx Number of terrain data points in the e-w direction
! for the DESIRED data resolution
! ndy Number of terrain data points in the n-s direction
! for the DESIRED data resolution
! nfile2 Filename of the data set to be accessed
! rlength Length of the record to be read from the data set
! readdat Array which will read the data file
!
! OUTPUT:
!
! nval Array in which the terrain values will be stored
!
!-----------------------------------------------------------------------
!
! Variable declarations
!
!-----------------------------------------------------------------------
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! INPUT:
!
!-----------------------------------------------------------------------
!
INTEGER :: n ! Specific value of the record number
INTEGER :: idim ! East-west readdat dimension
INTEGER :: jdim ! North-south readdat dimension
INTEGER :: ndx ! East-west nval dimension
INTEGER :: ndy ! North-south nval dimension
CHARACTER (LEN=* ) :: nfile2 ! Data set filename
INTEGER :: rlength ! Length of the data record
INTEGER*2 readdat(idim,jdim) ! Array which will read the data file
!
!-----------------------------------------------------------------------
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
INTEGER*2 nval(ndx,ndy) ! Array of terrain data for the block
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: itema,itemb
REAL :: tema,temb,temc,temd
INTEGER :: k,l,istat
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
OPEN(11,FILE=nfile2,ACCESS='direct',FORM='unformatted', &
STATUS='old',RECL=rlength,IOSTAT=istat)
IF( istat /= 0) THEN
WRITE(6,'(/a,/a,a/)') &
'Error occured when opening terrain data file.', &
nfile2,', Program stopped in READTER.'
STOP
END IF
READ(11,REC=n,ERR=991) readdat
CLOSE (11)
!
!-----------------------------------------------------------------------
!
! Filling the nval array
!
!-----------------------------------------------------------------------
!
IF(idim < ndx)THEN ! if true, nval will be filled stepwise
temc=REAL(idim/REAL(ndx))
temd=REAL(jdim/REAL(ndy))
DO l=1,ndy
DO k=1,ndx
tema=REAL((k-1)*temc)
temb=REAL((l-1)*temd)
itema=INT(tema)
itemb=INT(temb)
nval(k,l)=readdat(1+itema,1+itemb)
END DO
END DO
ELSE ! direct fill...
DO l=1,ndy
DO k=1,ndx
nval(k,l)=readdat(k,l)
END DO
END DO
END IF
RETURN
991 PRINT*,'Error occurred when reading data ',nfile2, &
' in READTER. Program stopped in READTER.'
STOP
END SUBROUTINE readter
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE BARNES ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE barnes(ipnum,ibgn,iend,jbgn,jend,nx,ny,ntx,nty, & 2
nbufx,nbufy,temak,zbt,za, &
zga, &
xg,yg,xob,yob,closi,closj,tem1,tem2)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Perform a Barnes analysis on a data base for a prespecified
! knot and gamma (multipass scheme). The user can specific knot, the
! influence factor, in arpstern.input or allow the program to calculate
! one given specified parameters (in setgrid). The cutoff radius is
!
! rc=(ln(rdnot) * KNOT)**0.5
!
! The Barnes weight function is w = EXP(-r**2/knot)
! where r is the distance between the analysis point and the data
! point.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Dan Weber
! 1/13/94
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! ipnum Current pass number
! ibgn Starting index for the analysis (e-w)
! iend Ending index for the analysis (e-w)
! jbgn Starting index for the analysis (n-s)
! jbgn Ending index for the analysis (n-s)
! nx Number of analysis grid points in the east-west direction
! ny Number of analysis grid points in the north-south direction
! ntx Limit of number of data points e-w
! nty Limit of number of data points n-w
! nbufx Number of additional data points outside the model
! domain used on the calculation of the boundaries. east-west
! nbufy Number of additional data points outside the model
! domain used on the calculation of the boundaries. north-south
! temak Working value of knot (set in setgrid subroutine)
!
! zbt Analyzed data at analysis points obtained
! from performing a barnes analysis at the data points.
! za Data for the buffered and analysis areas (m)
!
! xg Analysis grid points east-west direction (grid units)
!
! yg Analysis grid points north-south direction (grid units)
!
! xob Data set grid point values east-west direction
! yob Data set grid point values north-south direction
!
! closi Array storing the data i index value corresponding to
! the closest analysis grid point
! closj Array storing the data j index value corresponding to
! the closest analysis grid point
!
! COMMON block variables:
!
! knot Shape parameter for the barnes weight function.
! gamma Response factor for multiple passes of the Barnes analysis
! ipass Total number of passes to be performed
!
! OUTPUT:
!
! zga Analyzed data field. (m)
!
!
! TEMPORARY ARRAYS
!
! tem1 Temporary working array (barnes weight function)
! tem2 Temporary working array
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!-----------------------------------------------------------------------
!
! INPUT:
!
!-----------------------------------------------------------------------
!
INTEGER :: ipnum ! Current pass number
INTEGER :: ibgn ! starting index for the analysis (e-w)
INTEGER :: iend ! ending index for the analysis (e-w)
INTEGER :: jbgn ! starting index for the analysis (n-s)
INTEGER :: jend ! ending index for the analysis (n-s)
INTEGER :: nx ! Number of analysis grid points (e-w)
INTEGER :: ny ! Number of analysis grid points (n-s)
INTEGER :: ntx ! Limit of data points e-w
INTEGER :: nty ! Limit of data poitns n-s
INTEGER :: nbufx ! Number of buffer data points e-w
INTEGER :: nbufy ! Number of buffer data points n-s
REAL :: temak ! Working value of knot
REAL :: zbt(ntx,nty) ! Background data field (not required)
! at the data points. (m)
REAL :: za(ntx,nty) ! Data base for use in the Barnes
! analysis for the buffered area. (m)
REAL :: xg(nx,ny) ! Analysis grid points east-west
! direction (analysis grid units)
REAL :: yg(nx,ny) ! Analysis grid points north-south
! direction (analysis grid units)
REAL :: xob(ntx,nty) ! Data set grid point values east-west
! direction
REAL :: yob(ntx,nty) ! Data set grid point values north-south
! direction
INTEGER :: closi(nx,ny) ! array storing the data i index value
! corresponding to the closest grid pt
INTEGER :: closj(nx,ny) ! array storing the data j index value
! corresponding to the closest grid pt
!
!-----------------------------------------------------------------------
!
! COMMON block variables:
!
!-----------------------------------------------------------------------
!
REAL :: knot ! Shape parameter for the barnes scheme
REAL :: gamma ! Response factor for multiple pass Barnes
INTEGER :: ipass ! Total number of passes to be performed
!
!-----------------------------------------------------------------------
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
REAL :: zga(nx,ny) ! Analyzed data base for use by the user (m)
!
!-----------------------------------------------------------------------
!
! TEMPORARY ARRAYS
!
!-----------------------------------------------------------------------
!
REAL :: tem1(nx,ny) ! Temporary working array (Barnes Weights)
REAL :: tem2(nx,ny) ! Temporary working array
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,m,n,in,jn
REAL :: w,tema,temb,temc
!
!-----------------------------------------------------------------------
!
! COMMON block definition:
!
!-----------------------------------------------------------------------
!
COMMON /barnesc/ knot,gamma,ipass
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
IF(ipnum == 1)THEN
tema=1.0
DO j=jbgn,jend
DO i=ibgn,iend
zga(i,j)=0.0
END DO
END DO
DO j=1,nty
DO i=1,ntx
zbt(i,j)=0.0
END DO
END DO
ELSE
tema=gamma
END IF
DO j=jbgn,jend
DO i=ibgn,iend
tem1(i,j)=0.0
tem2(i,j)=0.0
END DO
END DO
!
!-----------------------------------------------------------------------
!
! BEGIN THE BARNES ANALYSIS***********
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Print the bounds of the barnes analysis and current pass number.
!
!-----------------------------------------------------------------------
!
PRINT *,'ibgn,iend,jbgn,jend',ibgn,iend,jbgn,jend
PRINT *,'current barnes pass ',ipnum
PRINT *,'closi(1,ny),closj(1,ny),nbufx,nbufy', &
closi(1,ny),closj(1,ny),nbufx,nbufy
tema=temak*tema
DO n=-nbufy,nbufy
DO m=-nbufx,nbufx
DO j=jbgn,jend
DO i=ibgn,iend
in=closi(i,j)+m
jn=closj(i,j)+n
temb=xg(i,j)-xob(in,jn)
temc=yg(i,j)-yob(in,jn)
w=EXP(-(temb*temb+temc*temc)/tema)
tem2(i,j)=(za(in,jn)-zbt(in,jn))*w+tem2(i,j)
tem1(i,j)=tem1(i,j)+w
END DO
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Summing up the weights from each observation inside the rc
! and dividing by the total weight for each analysis grid point.
!
!-----------------------------------------------------------------------
!
DO j=jbgn,jend
DO i=ibgn,iend
zga(i,j)=zga(i,j)+tem2(i,j)/tem1(i,j)
END DO
END DO
PRINT *,'leaving the barnes subroutine'
RETURN
END SUBROUTINE barnes
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE RMSDIF ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE rmsdif(ibgn,iend,jbgn,jend,ntx,nty,tema,temb,temc, & 1,1
temd,glab,f1,f2,diff)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate and plot the root mean square difference between
! two data fields.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Dan Weber
! 1/12/94
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! ibgn Starting i index for the differencing loop
! iend Ending i index for the differencing loop
! jbgn Starting j index for the differencing loop
! jend Ending j index for the differencing loop
!
! ntx Dimension in the e-w direction for the f1,f2, diff arrays
! nty Dimension in the n-s direction for the f1,f2, diff arrays
!
! tema Latitude of the northwest corner of the area
! temb Longitude of the northwest corner of the area
! temc Latitude of the southwest corner of the area
! temd Longitude of the southwest corner of the area
! glab Label for the difference field plot
!
! f1 First data array
! f2 Second data array
!
!
! OUTPUT:
!
! diff Root mean square difference field between f1 and f2
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: ibgn ! Starting i index for the differencing loop
INTEGER :: iend ! Ending i index for the differencing loop
INTEGER :: jbgn ! Starting j index for the differencing loop
INTEGER :: jend ! Ending j index for the differencing loop
INTEGER :: ntx ! Dimension in the e-w direction for f1 and f2
INTEGER :: nty ! Dimension in the n-s direction for f1 and f2
REAL :: tema ! Lat of the northwest corner of the area
REAL :: temb ! Lon of the northwest corner of the area
REAL :: temc ! Lat of the southwest corner of the area
REAL :: temd ! Lon of the southwest corner of the area
CHARACTER (LEN=65) :: glab ! Label for the difference field plot
REAL :: f1(ntx,nty) ! Array containing field number 1
REAL :: f2(ntx,nty) ! Array containing field number 2
REAL :: diff(ntx,nty) ! Root mean square difference field
! diff=sqrt((f1-f2)*(f1-f2))
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j
REAL :: cint
!
!-----------------------------------------------------------------------
!
! Defining COMMON blocks:
!
!-----------------------------------------------------------------------
!
REAL :: trulat (2) ! Latitude of true map projection
! (deg N)
REAL :: trulon ! Longitude of true map projection
! (deg E)
REAL :: sclfct ! Map scale factor (m)
COMMON /projectc/ trulat,trulon,sclfct
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
PRINT *,'ibgn,iend,jbgn,jend in rmsdiff'
PRINT *,ibgn,iend,jbgn,jend
DO j=jbgn,jend
DO i=ibgn,iend
diff(i,j)=SQRT((f1(i,j)-f2(i,j))*(f1(i,j)-f2(i,j)))
END DO
END DO
cint=-4.
CALL plot((iend-ibgn)+1,(jend-jbgn)+1,ntx,8,trulat,trulon,cint, &
tema,temb,temc,temd,glab,diff(ibgn,jbgn))
RETURN
END SUBROUTINE rmsdif
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE RESPONSE ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE response(tema,temb,temc,temak,pi) 1
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate the final response function/curve for the Barnes analysis
! Details of this method of calculating the response function
! can be found in Koch et.al. JCAM Sept. 1983.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Dan Weber
! 12/6/93
!
!-----------------------------------------------------------------------
!
! INPUT :
!
!
! tema Wdn*max grid spacing(dx or dy)*sclfct
! temb max grid spacing*sclfct
! temc max grid spacing
! temak Working value of knot (set in setgrid)
! pi The value pi
!
!
! COMMON BLOCK VARIABLES:
!
! knot Barnes response shape factor
! gamma Multipass shape adjustment
! ipass Total passes for the Barnes scheme
!
!
! OUTPUT:
!
! d1star response function for the barnes analysis
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! INPUT:
!
!-----------------------------------------------------------------------
!
REAL :: tema ! Wdn*max spacing*sclfct
REAL :: temb ! Max spacing*sclfct
REAL :: temc ! max data or grid spacing
REAL :: temak ! Working value of knot (set in setgrid)
REAL :: pi ! The value pi
INTEGER :: inum ! Number of points used to define the response
! curve.
PARAMETER (inum=35)
REAL :: dnot(inum) ! response function for the first pass of
! the barnes analysis
REAL :: xplot(inum) ! x axis for response curve plot
!
!-----------------------------------------------------------------------
!
! COMMON BLOCK VARIABLES:
!
!-----------------------------------------------------------------------
!
REAL :: knot ! Barnes response shape factor
REAL :: gamma ! Multipass shape adjustment
INTEGER :: ipass ! Total passes for the Barnes scheme
!
!-----------------------------------------------------------------------
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
REAL :: dnstar(inum) ! response function for the barnes analysis
!
!-----------------------------------------------------------------------
!
! TEMPROARY ARRAYS:
!
!-----------------------------------------------------------------------
!
REAL :: tem1(inum)
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,n
REAL :: alamnot,akstar,alam,alamstr,tem
CHARACTER (LEN=75) :: glab
!
!-----------------------------------------------------------------------
!
! COMMON BLOCK DEFINITION:
!
!-----------------------------------------------------------------------
COMMON /barnesc/ knot,gamma,ipass
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
alamnot=tema
akstar=temak/(alamnot*alamnot)
PRINT *,'tema,temak in response subroutine'
PRINT *,tema,temak
!
!-----------------------------------------------------------------------
!
! Calculate the initial response....
!
!-----------------------------------------------------------------------
!
DO n=1,inum
alam=REAL(n)*temb
alamstr=alam/alamnot
tem= -akstar*(pi/alamstr)*(pi/alamstr)
dnot(n)=EXP(tem)
xplot(n)=n
END DO
IF(ipass == 1)THEN ! set dnot = dnstar and plot
DO n=1,inum
dnstar(n)=dnot(n)
END DO
ELSE IF(ipass == 2)THEN ! calculate D1* = dnstar
DO n=1,inum
tem=dnot(n)**gamma
dnstar(n)=dnot(n)+tem-tem*dnot(n)
END DO
ELSE ! calculate Dn* = dnstar
DO n=1,inum
dnstar(n)=0.0
tem1(n)=1.0
END DO
DO i=2,ipass-1
tem=gamma*(i-1)
DO n=1,inum
tem1(n)=(1.0-dnot(n)**tem)*tem1(n)
END DO
tem=i*gamma
DO n=1,inum
dnstar(n)=(dnot(n)**tem)*tem1(n)+dnstar(n)
END DO
END DO
DO n=1,inum
dnstar(n)=dnot(n)+(1.0-dnot(n))*((dnot(n)**gamma)+dnstar(n))
END DO
END IF
!-----------------------------------------------------------------------
!
! Print the response for the first and last passes.
!
!-----------------------------------------------------------------------
PRINT *,'wavelength (dx) first pass(dnot) last pass(dnotstar)'
DO n=1,inum
PRINT *,n,dnot(n),dnstar(n)
END DO
!-----------------------------------------------------------------------
!
! Set up the labels for the plot axis
!
!-----------------------------------------------------------------------
CALL agsetc('LABEL/NAME.','T')
CALL agseti('LINE/NUMBER.',100)
CALL agseti('LINE/MAXIMUM.',70)
CALL agsetc('LABEL/NAME.','B')
CALL agseti('LINE/NUMBER.',-100)
CALL agsetc('LINE/TEXT.','Wavelength (x*max grid spacing)')
CALL agsetc('LABEL/NAME.','L')
CALL agseti('LINE/NUMBER.',100)
CALL agsetc('LINE/TEXT.','Response (g/f)')
WRITE(glab,70) ipass,gamma,temc
70 FORMAT('Barnes Response Function',' Pass # =',i2, &
' Gamma =',f5.2,' max. spacing=',f8.1)
CALL ezxy(xplot,dnstar,inum,glab)
RETURN
END SUBROUTINE response