funcs.f
************************************************************************
* Collection of handy functions and subroutines
* Just include in makefile to have available everywhere
* In and output/return is shortly described for each procedure
5: * The file is roughly subdivided in the following sections:
* - string manipulation
* - unit convertion
* - sorting
* - calendar functionality
10: * - random number generation
* - statistical functionality
* - grid transformation and interpolation
************************************************************************
15:
SUBROUTINE wordCaseConvert (string)
************************************************************************
* Function to convert an string (e.g. all uppercase) into lower case (with
* every word starting with uppercase letter).
20: * NOTE: Only normal alphabet characters supported.
*
* In+Output:
* - string (String): string to convert/converted string
************************************************************************
25: CHARACTER*(*) string
INTEGER strlen
CHARACTER*1 letter
INTEGER ishift, i
30: LOGICAL nextUp
nextUp = .true.
ishift=ICHAR('a')-ICHAR('A') !case shift difference
do i=1,strlen(string)
35: letter=string(i:i)
if ( ('a' <= letter .and. letter <= 'z')
c .or. ('A' <= letter .and. letter <= 'Z' ) ) then ! it's a character
if (nextUp .and. 'a' <= letter .and. letter <= 'z') then ! lower case that needs to go up
40: string(i:i)=CHAR(ICHAR(letter)-ishift)
elseif (.not.nextUp .and. 'A'<=letter .and. letter<='Z') then ! upper case that needs to go down
string(i:i)=CHAR(ICHAR(letter)+ishift)
else ! stays at is is
string(i:i) = letter
45: endif
nextUp = .false.
else
nextUp = .true.
endif
50: enddo
RETURN
END
55: INTEGER FUNCTION strlen (str)
************************************************************************
* Function to retrieve the length of a string w/o trailing blanks
*
* Input:
60: * - str (String): string whose length is to determine
* Return (Integer): string length
************************************************************************
IMPLICIT none
CHARACTER*(*) str
65: INTEGER i
DO 15, i = LEN(str), 1, -1
IF (str(i:i) .NE. ' ') GO TO 20
15 CONTINUE
20 strlen = i
70: RETURN
END
REAL FUNCTION tempC2F (p_tempC)
75: ************************************************************************
* Converts a celsius temperature to fahrenheit
* Input: p_tempC (Real): temperature in Celsius
* Return (Real): temperature in Fahrenheit
************************************************************************
80: REAL p_tempC
tempC2F = (9.0 / 5.0) * p_tempC + 32.0
END
REAL FUNCTION tempF2C (p_tempF)
85: ************************************************************************
* Converts a fahrenheit temperature to celsius
* Input: p_tempF (Real): temperature in Fahrenheit
* Return (Real): temperature in Celsius
************************************************************************
90: REAL p_tempF
tempF2C = (5.0 / 9.0) * (p_tempF - 32.0)
END
REAL FUNCTION tempK2C (p_tempK)
95: ************************************************************************
* Converts a kelvin temperature to celsius
* Input: p_tempK (Real): temperature in Kelvin
* Return (Real): temperature in Celsius
************************************************************************
100: REAL p_tempK
tempK2C = p_tempK - 273.15
END
REAL FUNCTION tempC2K (p_tempC)
105: ************************************************************************
* Converts a Celsius temperature to Kelvin
* Input: p_tempC (Real): temperature in Celsius
* Return (Real): temperature in Kelvin
************************************************************************
110: REAL p_tempC
tempC2K = p_tempC + 273.15
END
115: SUBROUTINE shell(n,a)
************************************************************************
* Sorts an array a(1:n) into ascending numerical order by Shell’s method
* (diminishing increment sort). n is input; a is replaced on output by its
* sorted rearrangement.
120: * SOURCE: Numerical Recipes in Fortran 77 Chapter 8.1
* Input:
* - n (Integer): array size
* In+Output:
* - a (Real[n]): unsorted input array/sorted output array
125: ************************************************************************
INTEGER n
REAL a(n)
INTEGER i,j,inc
130: REAL v
inc=1 !starting increment
1 inc=3*inc+1
if(inc.le.n)goto 1
2 continue !Loop over the partial sorts.
135: inc=inc/3
do i=inc+1,n ! Outer loop of straight insertion.
v=a(i)
j=i
3 if(a(j-inc).gt.v)then ! Inner loop of straight insertion.
140: a(j)=a(j-inc)
j=j-inc
if(j.le.inc) goto 4
goto 3
endif
145: 4 a(j)=v
enddo
if(inc.gt.1) goto 2
return
END
150:
SUBROUTINE createYearsArray (p_a, p_yearOffset, o_years)
************************************************************************
* Subroutine to create an array with years according to the year offset
155: * Input:
* - p_a (Integer): Number of years to create
* - p_yearOffset (Integer): Offset year for the years array
* Output:
* - o_years (Integer[p_a]): years array
160: ************************************************************************
INTEGER p_a, p_yearOffset
INTEGER o_years(p_a)
do i=1, p_a
o_years(i) = p_yearOffset + i-1
165: enddo
END
INTEGER FUNCTION daysInMonth(p_year, p_month)
170: ************************************************************************
* Function to finds the number of days in a certain month
* Input:
* - p_year (Integer): The year in which the month is (for leap years)
* - p_month (Integer): The month for which to find the number of days
175: * Return (Integer): number of days in the month
************************************************************************
INTEGER p_year, p_month
select case (p_month)
180: case(1, 3, 5, 7, 8, 10, 12) !Jan, Mar, May, Jul, Aug, Oct, Dec
daysInMonth = 31
return
case(4, 6, 9, 11) !Apr, Jun, Sep, Nov
daysInMonth = 30
185: return
case(2) !Feb
* if the year is divisable by 400 or divisable by 4 and not divisble by 100 it is a leap year
if (
c (p_year / 400. - int(p_year / 400.) .eq. 0) .or.
190: c (
c (p_year / 4. - int(p_year / 4.) .eq. 0) .and.
c (p_year / 100. - int(p_year / 100.) .gt. 0)
c ) ) then
daysInMonth = 29
195: else
daysInMonth = 28
endif
return
end select
200: END
FUNCTION ran1(idum)
************************************************************************
205: * 'Minimal' random number generator of Park and Miller with Bays-Durham shuffle and
* added safeguards. Returns a uniform random deviate between 0.0 and 1.0 (exclusive of
* the endpoint values). Call with idum a negative integer to initialize; thereafter, do not
* alter idum between successive deviates in a sequence. RNMX should approximate the largest
* floating value that is less than 1.
210: * SOURCE: Numerical Recipes in Fortran 77 Chapter 7.1
* Input: idum (Integer): initialization integer
* Return (Real) ]0;1[ : random number
*
************************************************************************
215: INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
REAL ran1,AM,EPS,RNMX
PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
c NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
220: INTEGER j,k,iv(NTAB),iy
SAVE iv,iy
DATA iv /NTAB*0/, iy /0/
if (idum.le.0.or.iy.eq.0) then ! Initialize.
idum=max(-idum,1) !Be sure to prevent idum = 0.
225: do j=NTAB+8,1,-1 !Load the shuffle table (after 8 warm-ups)
k=idum/IQ
idum=IA*(idum-k*IQ)-IR*k
if (idum.lt.0) idum=idum+IM
if (j.le.NTAB) iv(j)=idum
230: enddo
iy=iv(1)
endif
k=idum/IQ !Start here when not initializing.
idum=IA*(idum-k*IQ)-IR*k !Compute idum=mod(IA*idum,IM) without overflows by Schrage's method
235: if (idum.lt.0) idum=idum+IM
j=1+iy/NDIV !Will be in the range 1:NTAB.
iy=iv(j) !Output previously stored value and refill the shuffle table
iv(j)=idum
ran1=min(AM*iy,RNMX) !Because users don’t expect endpoint values.
240: return
END
SUBROUTINE corCoef (p_n, p_x, p_y, o_r)
245: ************************************************************************
* Subroutine to calculate the linear correlation coefficient of two same
* size samples p_x and p_y (simplified for these needs from the original
* routine)
* SOURCE (modified): Numerical Recipes in Fortran 77 Chapter 14.5
250: * Input:
* - p_n (Integer): size of the sample arrays
* - p_x (Real[p_n]): sample data set 1
* - p_y (Real[p_n]): sample data set 2
* Output:
255: * - o_r (Real) [-1;1]: linear correlation coefficient
*
************************************************************************
IMPLICIT none
INTEGER p_n ! size of variables
260: REAL tin ! for regularization of unusal case of complete correlation
PARAMETER(tin=1.e-20)
REAL p_x(p_n), p_y(p_n) ! data variables
REAL o_r ! correlation coeffiction output
INTEGER i ! process var
265: REAL ax, ay, sxx, syy, sxy, xt, yt
ax=0.
ay=0.
do i=1,p_n
ax=ax+p_x(i)
270: ay=ay+p_y(i)
enddo
ax=ax/real(p_n)
ay=ay/real(p_n)
sxx=0.
275: syy=0.
sxy=0.
do i=1, p_n
xt=p_x(i)-ax
yt=p_y(i)-ay
280: sxx=sxx+xt**2
syy=syy+yt**2
sxy=sxy+xt*yt
enddo
o_r=sxy/(sqrt(sxx*syy)+tin)
285: return
END
SUBROUTINE transform2RegGrid(p_oldX, p_oldY,
290: c p_lats, p_longs, p_var,
c p_newX, p_newY,
c o_lats, o_longs, o_var)
************************************************************************
* Transformation of geographical data from an irregular to a regular grid
295: * (i.e. 1-dimensional longitudes and lattitudes) by interpolation. The size
* of the new grid needs to be passed in but the lattitude and longitude
* extrema will be automatically calculated by the routine (i.e. no matter what
* the shape of the original data grid is, the routine tries to cut out a
* large rectangle for the regular grid).
300: *
* Input:
* - p_oldX (Integer): old data grid size
* - p_oldY (Integer): old data grid size
* - p_lats (Real[p_oldX, p_oldY]): latitudes of the data points
305: * - p_longs (Real[p_oldX, p_oldY]): longitudes of the data points
* - p_var (Real[p_oldX, p_oldY]): data points
* - p_newX (Integer): new data grid size
* - p_newY (Integer): new data grid size
* Output:
310: * - o_lats (Real[p_newY]): new 1-Dimensional Lattitudes of regular grid
* - o_longs (Real[p_newX]]: new1-Dimensional Longitudes of regular grid
* - o_var (Real[p_newX, p_newY]): interpolated data on new (regular) grid
************************************************************************
IMPLICIT none
315:
INTEGER p_oldX, p_oldY !old grid size
INTEGER p_newX, p_newY !new grid size
REAL p_longs(p_oldX, p_oldY), p_lats(p_oldX, p_oldY) ! coordinates
320: REAL p_var(p_oldX, p_oldY)
REAL o_lats(p_newY), o_longs(p_newX) ! transformed regular coordinates
REAL o_var(p_newX, p_newY)
325: REAL minLat, maxLat, minLon, maxLon ! coordinate system scaling
REAL iVar ! interpolate value
REAL interpolatePt ! interpolate function
INTEGER i,j ! Looping vars
INTEGER cc ! not interpolatable vars
330:
! find new grid lat/lon minima and maxima:
minLat = 1000.
maxLat = -1000.
minLon = 1000.
335: maxLon = -1000.
! find extreme points of earth grid
do i=1,p_oldX ! looping data arrays, thus data grid vars (p_oldX, p_oldY)
do j=1,p_oldY
if (p_lats(i,j) + p_longs(i,j) < minLat + minLon ) then
340: minLat = p_lats(i,j)
minLon = p_longs(i,j)
endif
if (p_lats(i,j) + p_longs(i,j) > maxLat + maxLon ) then
maxLat = p_lats(i,j)
345: maxLon = p_longs(i,j)
endif
enddo
enddo
350: ! restrict maximum longitude to one within the existent range
minLon = 1000.
maxLon = -1000.
minLat = minLat + 4 ! move frame slightly up
maxLat = maxLat ! extend frame towards the north (7.5 to have scandinavia included on HC derived grid)
355: do i=1,p_oldX ! looping data arrays, thus data grid vars (p_oldX, p_oldY)
do j=1,p_oldY
if (p_lats(i,j) <= minLat .and. p_longs(i,j) > maxLon)
c maxLon = p_longs(i,j)
if (p_lats(i,j) <= minLat .and. p_longs(i,j) < minLon)
360: c minLon = p_longs(i,j)
enddo
enddo
! offset for now
365: minLat = minLat + 0
maxLat = maxLat - 0
minLon = minLon + 0
maxLon = maxLon - 0
370: ! fill net cdf arrays
do j=1,p_newY
o_lats(j)=minLat +
c (j-1) * (maxLat-minLat)/real(p_newY-1)
enddo
375:
do i=1, p_newX
o_longs(i) =minLon +
c (i-1) * (maxLon-minLon)/real(p_newX-1)
enddo
380:
! find values through interpolation
cc = 0
do j=1, p_newY
do i=1, p_newX
385: iVar = interpolatePt(o_longs(i), o_lats(j),
c p_oldX, p_oldY, p_longs, p_lats, p_var,
c 1., cc)
!if (.not. (iVar == 0) ) iVar = 1./iVar
o_var(i,j) = iVar
390: enddo
enddo
print *, "new grid: min/max lon ", minLon, maxLon
print *, "new grid: min/max lat ", minLat, maxLat
395: print *, cc, " points not valid for new grid"
END
400: REAL FUNCTION interpolatePt (p_lon, p_lat,
c p_oldX, p_oldY, p_longs, p_lats, p_var,
c p_fillVal, p_count )
************************************************************************
* Interpolation of one single point on a grid from another grid and its data
405: * Input:
* - p_lon (Real): longitude of the point to interpolate
* - p_lat (Real): latitude of the point to interpolate
* - p_oldX (Integer): grid size of the grid from which to take the data to interpolate
* - p_oldY (Integer): grid size of the grid from which to take the data to interpolate
410: * - p_longs (Real[p_oldX, p_oldY]): longitudes of the available data points
* - p_lats (Real[p_oldX, p_oldY]): latitudes of the available data points
* - p_var (Real[p_oldX, p_oldY]): available data points
* - p_fillVal (Real): value that should be used for any points that cannot
* successfully be interpolated from the available data
415: * In+Output: p_count(Integer): is increased by one if the fillValue needs to be used
* Return (Real): The interpolated data value of the point
************************************************************************
IMPLICIT none
420: REAL p_lon, p_lat ! coordinates (on new grid) to interpolate
INTEGER p_oldX, p_oldY !old grid size
REAL p_longs(p_oldX, p_oldY), p_lats(p_oldX, p_oldY) ! old coordinates
REAL p_var(p_oldX, p_oldY) ! variable values
425: REAL p_fillVal ! value for grid points that are not interpolatable
INTEGER p_count ! how many points not gridable
INTEGER minDistLoc(6,2)
REAL minDist(6), dist ! 4 points for interpolation
430: REAL a, b, c, iLat1, iLat2 ! interpolation vars
REAL iVar1, iVar2, iVar ! interpolation vars
INTEGER quadrant ! quadrant ID
INTEGER l, k !looping vars
435:
minDist(1) = -1. ! distance in Q1
minDist(2) = -1. ! distance in Q2
minDist(3) = -1. ! distance in Q3
minDist(4) = -1. ! distance in Q4
440: minDist(5) = -1. ! distance on pos y axis
minDist(6) = -1. ! distance on neg y axis
! find closest 4 grid points in the 4 quadrants
do l=1, p_oldY
445: do k=1, p_oldX
if (p_lon == p_longs(k,l) .and. p_lat == p_lats(k,l) ) then
! point matches directly, no interpolation needed
interpolatePt = p_var(k,l)
return
450: endif
dist = (p_lon-p_longs(k,l))**2+
c (p_lat-p_lats(k,l))**2
if (p_lats(k,l) > p_lat ) then
if (p_longs(k,l) > p_lon ) then
455: quadrant = 1 ! upper right quadrant
else if ( p_longs(k,l) == p_lon ) then
quadrant = 5 ! on pos y axis
else
quadrant = 2 ! upper left quadrant
460: endif
else
if (p_longs(k,l) > p_lon ) then
quadrant = 4 ! lower right quadrant
else if ( p_longs(k,l) == p_lon ) then
465: quadrant = 6 ! on neg y axis
else
quadrant = 3 ! lower left quadrant
endif
endif
470: if ( (dist .lt. minDist(quadrant)) .or.
c (minDist(quadrant) .eq. -1.) ) then
minDist(quadrant) = dist
minDistLoc(quadrant,1) = k
minDistLoc(quadrant,2) = l
475: endif
enddo
enddo
! check if sufficient points exist
480: ! either 5 or 1 and 2, and either 6 or 3 and 4 must exist
if (
c ( minDist(5) == -1. .and.
c ( minDist(1) == -1. .or. minDist(2) == -1. ) ) .or.
c ( minDist(6) == -1. .and.
485: c ( minDist(3) == -1. .or. minDist(4) == -1. ) ) ) then
p_count = p_count + 1
interpolatePt = p_fillVal
return
endif
490: * In the following interpolation calculations, the suffix 0 identifies
* attributes of the netCdfGrid point (i.e. the origin of the calculations)
* where as 1 through 4 are the four nearest points in the four quadrants
* counterclockwise from +lon/+lat to -lon/+lat
* i prefix identifies interpolated values (intermediate and final results
495: * of the interpolation)
if ( .not. (minDist(5) == -1.) .and.
c ( ( minDist(1) == -1. .or. minDist(2) == -1. ) .or.
c ( minDist(5) <= minDist(1) .and. minDist(5) <= minDist(2))))
c then
500: iLat1 = p_lats(minDistLoc(5,1), minDistLoc(5,2))
iVar1 = p_var(minDistLoc(5,1), minDistLoc(5,2))
else ! no better distance on the axis
! a = (lon0 - lon2)/(lon1 - lon2)
a = (p_lon - p_longs(minDistLoc(2,1), minDistLoc(2,2)))/
505: c (p_longs(minDistLoc(1,1), minDistLoc(1,2)) -
c p_longs(minDistLoc(2,1), minDistLoc(2,2)))
! iLat1 = lat2 + a*(lat1-lat2)
iLat1 = p_lats(minDistLoc(2,1), minDistLoc(2,2)) +
c a*(p_lats(minDistLoc(1,1), minDistLoc(1,2)) -
510: c p_lats(minDistLoc(2,1), minDistLoc(2,2)))
! iVal1 = val2 + a*(val1-val2)
iVar1 = p_var(minDistLoc(2,1), minDistLoc(2,2)) +
c a*(p_var(minDistLoc(1,1), minDistLoc(1,2)) -
c p_var(minDistLoc(2,1), minDistLoc(2,2)))
515: endif
if ( .not. (minDist(6) == -1.) .and.
c ( ( minDist(3) == -1. .or. minDist(4) == -1. ) .or.
c ( minDist(6) <= minDist(3) .and. minDist(6) <= minDist(4))))
520: c then
iLat2 = p_lats(minDistLoc(6,1), minDistLoc(6,2))
iVar2 = p_var(minDistLoc(6,1), minDistLoc(6,2))
else ! no better distance on the axis
! b = (lon0 - lon3)/(lon4 - lon3)
525: b = (p_lon - p_longs(minDistLoc(3,1), minDistLoc(3,2)))/
c (p_longs(minDistLoc(4,1), minDistLoc(4,2)) -
c p_longs(minDistLoc(3,1), minDistLoc(3,2)))
! iLat2 = lat3 + b*(lat4-lat3)
iLat2 = p_lats(minDistLoc(3,1), minDistLoc(3,2)) +
530: c b*(p_lats(minDistLoc(4,1), minDistLoc(4,2)) -
c p_lats(minDistLoc(3,1), minDistLoc(3,2)))
! iVal2 = val3 + b*(val4-val3)
iVar2 = p_var(minDistLoc(3,1), minDistLoc(3,2)) +
c b*(p_var(minDistLoc(4,1), minDistLoc(4,2)) -
535: c p_var(minDistLoc(3,1), minDistLoc(3,2)))
endif
! c = (lat0 - iLat2)/(iLat1 - iLat2)
c = (p_lat - iLat2) /
540: c (iLat1 - iLat2)
! iVal = iVal2 + c * (iVal1 - iVal2)
iVar = iVar2 +
c c*(iVar1 - iVar2)
interpolatePt = iVar
545: END
SUBROUTINE polcoe(x,y,n,cof)
************************************************************************
550: * Coefficients of interpolating polynomials
* Given arrays x(1:n) and y(1:n) containing a tabulated function yi=f(xi),
* this routine returns an array of coefficients cof(1:n), such that
* yi = sum cofj xi^j-1
* SOURCE: Numerical Recipes in Fortran 77 Chapter 3.5
555: * Input:
* - x (Real[n]): x values of tabulated function
* - y (Real[n]): y values of tabulated function
* - n (Integer): number of data points/coefficients to be searched
* Output:
560: * - cof (Real[n]): coefficients of the interpolating polynomial
************************************************************************
INTEGER n, NMAX
REAL cof(n), x(n), y(n)
PARAMETER (NMAX=15)
565: INTEGER i,j,k
REAL b,ff,phi,s(NMAX)
do i=1,n
s(i) = 0.
cof(i) = 0.
570: enddo
s(n) = -x(1)
do i=2,n
do j=n+1-i,n-1
s(j)=s(j)-x(i)*s(j+1)
575: enddo
s(n)=s(n)-x(i)
enddo
do j=1,n
phi=n
580: do k=n-1,1,-1
phi=k*s(k+1)+x(j)*phi
enddo
ff=y(j)/phi
b=1.
585: do k=n,1,-1
cof(k)=cof(k)+b*ff
b=s(k)+x(j)*b
enddo
enddo
590: return
end
Info Section
Warning: externals (function calls) may not be acurate
back to top
f2html v0.3 (C) 1997,98 Beroud Jean-Marc. Fri Aug 11 17:54:58 CEST 2006