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