C******************* MARCH 2000 ADDITION *******************
      SUBROUTINE THERMIK(TREAD,SREAD,PREAD,RALPHA,RBETA)
C    CALCULATES THE THERMAL EXPANSION COEFFICIENT (RALPHA)
C                      [ (degrees C)**-1 ]
C    CALCULATES THE SALINITY CONTRACTION COEFFICIENT (RBETA)
C                       [ pss**-1 ]
C    Input and Output parameters are REAL
C    ALL COMPUTATIONS in DOUBLE PRECISION
C     **REFERENCE**
C    UNESCO(1987); INTERNATIONAL OCEANOGRAPHIC TABLES 	
C    NO.40
C    JA found ERRATA
C    last term in equation (18) at page 104 must be as followed
C     ... (c1 + c2*t) * S**1.5 (NO POWERS for CONSTANTS c1 and c2!!!)
C
C-----------------------------------------------------------------------        
C     THIS SUBROUTINE CALCULATES THE DENSITY OF SEAWATER USING THE NEW          
C     STANDARD EQUATION OF STATE RECOMMENDED BY UNESCO(1981).                   
C     **UNITS**                                                                 
C       IN-SITU TEMPERATURE(T): DEGREES CENTIGRADE                              
C       SALINITY(S): PERMIL                                                     
C       PRESSURE(P): DECIBARS                                                   
C       DENSITY(RHO): KILOGRAMS PER CUBIC METER                                 
C     **REFERENCE**                                                             
C       UNESCO(1981); TENTH REPORT OF THE JOINT PANEL ON OCEANOGRAPHIC          
C       TABLES AND STANDARDS. UNESCO TECHNICAL PAPERS IN MARINE SCI.            
C       NO.36. UNESCO,PARIS.                                                    
C-----------------------------------------------------------------------        
C                                                                               
      DOUBLE PRECISION RW,RSTO,RHO,KW,KSTO,KSTP                                                         
      DOUBLE PRECISION DAW,DA,DBW,DB,DKW,DKSTO,DKSTP,KMON,ALPHA
      DOUBLE PRECISION BETA,DRSTOS,DAS,DBS,DKSTOS,DKSTPS
C                                                                               
      T=DBLE(TREAD)
      S=DBLE(SREAD)
      P=DBLE(PREAD*0.1)                                                               
      RW =     999.842594D0 + 6.793952D-2*T - 9.095290D-3*T**2               
     1        + 1.001685D-4*T**3 - 1.120083D-6*T**4 + 6.536332D-9*T**5          
C                                                                               
      RSTO = RW + (0.824493D0 - 4.0899D-3*T + 7.6438D-5*T**2                    
     1        - 8.2467D-7*T**3 + 5.3875D-9*T**4) * S                            
     2       + (-5.72466D-3 + 1.0227D-4*T - 1.6546D-6*T**2) * S**1.5            
     3       + 4.8314D-4 * S**2                                                 
C                                                                               
      KW =19652.21D0+148.4206D0*T-2.327105D0*T**2 + 1.360477D-2*T**3         
     1        - 5.155288D-5*T**4                                                
C                                                                               
      KSTO =KW+(54.6746D0-0.603459D0*T+1.09987D-2*T**2                     
     1        - 6.1670D-5*T**3) * S                                             
     2       + (7.944D-2 + 1.6483D-2*T - 5.3009D-4*T**2) * S**1.5               
C                                                                               
      KSTP =KSTO+(3.239908D0+1.43713D-3*T+1.16092D-4*T**2                
     1        - 5.77905D-7*T**3) * P                                            
     2       + (2.2838D-3 - 1.0981D-5*T - 1.6078D-6*T**2) * P * S               
     3       + 1.91075D-4 * P * S**1.5                                          
     4       + (8.50935D-5 - 6.12293D-6*T + 5.2787D-8*T**2) * P**2              
     5       + (-9.9348D-7 + 2.0816D-8*T + 9.1697D-10*T**2) * P**2 * S          
C                                                                               
      RHO =    RSTO/(1.0-P/KSTP)                                                
C
C JA ADDITION 03/17/00, THERMAL EXPANSION COEFFICIENT
C                                                                               
      DRW = 6.793952D-2 - 1.819058D-2*T + 3.005055D-4*T**2
     1     - 4.480332D-6*T**3 + 3.268166D-8*T**4
      DRSTO = DRW + (-4.0899D-3 + 1.52876D-4*T - 2.470401D-6*T**2 
     1       +2.155D-8*T**3 )*S + (1.0227D-4 - 3.3092D-6*T)*S**1.5
C
      DAW = 1.43713D-3 + 2.32184D-4*T - 1.733715D-6*T**2 
      DA  = DAW + (-1.0981D-5 - 3.2156D-6*T) * S
      DBW = -6.12293D-6 + 1.05574D-7*T
      DB  = DBW + (2.0816D-8 + 1.83394D-9*T) * S
      DKW = 148.4206D0 - 4.6542D0*T + 4.081431D-2*T**2
     1     -2.0621152D-4*T**3
      DKSTO = DKW +(-0.603459D0 + 2.19974D-2*T - 1.8501D-4*T**2)*S
     1       + (1.6483D-2 - 1.06018D-3*T) * S**1.5
      DKSTP = DKSTO + DA * P + DB * P**2
C
      KMON=1.D0-P/KSTP
      ALPHA= DRSTO/KMON - ((RSTO*P*DKSTP)/KSTP/KMON)/KSTP/KMON
      ALPHA = -1.D0*ALPHA/RHO
      RALPHA=REAL(ALPHA)
C
C SALINITY CONTRACTION COEFFICIENT
C
      DRSTOS = 8.24493D-1 - 4.0899D-3 * T + 7.6438D-5 * T**2  
     1        - 8.2467D-7 * T**3 + 5.3875D-9 * T**4
     2        + 1.5D0 * (-5.72466D-3 + 1.0227D-4 * T
     3        - 1.6546D-6 * T**2) * S**0.5
     4        + 9.6628D-4 * S

      DAS = 2.2838D-3 - 1.0981D-5*T - 1.6078D-6*T**2
     1     +2.866125D-4*S**0.5
      DBS = -9.9348D-7 + 2.0816D-8*T + 9.1697D-10*T**2
      DKSTOS = 54.6746D0 - 0.603459D0*T + 1.09987D-2*T**2
     1        - 6.1670D-5*T**3 + 1.5D0*(7.944D-2
     2        + 1.6483D-2*T - 5.3009D-4*T**2)*S**0.5
      DKSTPS = DKSTOS + DAS * P + DBS * P**2
C
      BETA = DRSTOS/KMON - ((RSTO*P*DKSTPS)/KSTP/KMON)/KSTP/KMON
      BETA = BETA/RHO
      RBETA = REAL(BETA)

      
      RETURN                                                                    
      END                                                                       
C******************* JULY 1997 ADDITION ********************
c     real function atg(s,t,p)
C adiabatic temperature gradient DEG C per DECIBAR
C Bryden, H., 1973, Deep-sea Res., 20, 401-408
C   PRESSURE     P    DECIBARS
C   TEMPERATURE  T    DEG C (IPTS-68)
C   SALINITY     S    (PSS-78)
C   ADIABATIC    ATG   DEG C/DECIBAR
Code after UNESCO #44 (1983), page 40
C***********************************************************
c     ds = s - 35.0
c     atg = (((-2.1687e-16*t + 1.8676e-14)*t-4.6206e-13)*p
c    1      ((2.7759e-12*t-1.1351e-10)*ds+((-5.4481e-14*t
c    2      +8.733e-12)*t-6.7795e-10)*t+1.8741e-8))*p
c    3      +(-4.2393e-8*t+1.8932e-6)*ds
c    4      +((6.6228e-10*t-6.836e-8)*t+8.5258e-6)*t+3.5803e-5
c     return
c     end
C******************** JULY 1997 ADDITION *********************
c***********************************************************
      function deplat(x,xlat)
c pressure in DECIBARS from depth in METERS using
c LEROY's method J.Acoust.Soc.Amer.,1969,46,N 1(part 2),216-226
c  depth     X      in  meters
c  latitude  XLAT   in  degrees
c  pressure  DEPLAT in  decibars ( 100 db =  MPa)
      y= sin(xlat)
      y=y*y
      y=1.02E-1 + (1 + 5.28E-3*y)* 1.00524E-2* x
      deplat=1.E2*(y + 2.475E-8* x* x)
      return
      end
c***********************************************************
      function depkal(x,xlat)
c pressure in DECIBARS from depth in METERS using
c KALASHNIKOV(1985) formula
c
      sqlat=xlat*xlat
      y=(1.00574E-2 - 6.5279E-8*xlat + 2.1889E-8*sqlat - 
     1   1.5709E-10*xlat*sqlat)*x + 2.3486E-8*x*x
      depkal=1.E2*y
      return
      end
      SUBROUTINE DIS(T,S,PREAD,RHO)                                             
C                                                                               
C-----------------------------------------------------------------------        
C     THIS SUBROUTINE CALCULATES THE DENSITY OF SEAWATER USING THE NEW          
C     STANDARD EQUATION OF STATE RECOMMENDED BY UNESCO(1981).                   
C     **UNITS**                                                                 
C       IN-SITU TEMPERATURE(T): DEGREES CENTIGRADE                              
C       SALINITY(S): PERMIL                                                     
C       PRESSURE(P): DECIBARS                                                   
C       DENSITY(RHO): KILOGRAMS PER CUBIC METER                                 
C     **REFERENCE**                                                             
C       UNESCO(1981); TENTH REPORT OF THE JOINT PANEL ON OCEANOGRAPHIC          
C       TABLES AND STANDARDS. UNESCO TECHNICAL PAPERS IN MARINE SCI.            
C       NO.36. UNESCO,PARIS.                                                    
C-----------------------------------------------------------------------        
C                                                                               
      REAL KW,KSTO,KSTP                                                         
C                                                                               
      P=PREAD*0.1                                                               
C                                                                               
      RW =     999.842594 + 6.793952E-2*T - 9.095290E-3*T**2                    
     1        + 1.001685E-4*T**3 - 1.120083E-6*T**4 + 6.536332E-9*T**5          
C                                                                               
      RSTO =   RW + (0.824493 - 4.0899E-3*T + 7.6438E-5*T**2                    
     1        - 8.2467E-7*T**3 + 5.3875E-9*T**4) * S                            
     2       + (-5.72466E-3 + 1.0227E-4*T - 1.6546E-6*T**2) * S**1.5            
     3       + 4.8314E-4 * S**2                                                 
C                                                                               
      KW =     19652.21 + 148.4206*T - 2.327105*T**2 + 1.360477E-2*T**3         
     1        - 5.155288E-5*T**4                                                
C                                                                               
      KSTO =   KW + (54.6746 - 0.603459*T + 1.09987E-2*T**2                     
     1        - 6.1670E-5*T**3) * S                                             
     2       + (7.944E-2 + 1.6483E-2*T - 5.3009E-4*T**2) * S**1.5               
C                                                                               
      KSTP =   KSTO + (3.239908 + 1.43713E-3*T + 1.16092E-4*T**2                
     1        - 5.77905E-7*T**3) * P                                            
     2       + (2.2838E-3 - 1.0981E-5*T - 1.6078E-6*T**2) * P * S               
     3       + 1.91075E-4 * P * S**1.5                                          
     4       + (8.50935E-5 - 6.12293E-6*T + 5.2787E-8*T**2) * P**2              
     5       + (-9.9348E-7 + 2.0816E-8*T + 9.1697E-10*T**2) * P**2 * S          
C                                                                               
      RHO =    RSTO/(1.0-P/KSTP)                                                
C                                                                               
      RETURN                                                                    
      END                                                                       
C                                                                               
      SUBROUTINE PTS(T,S,P,PT)                                                  
C                                                                               
C-----------------------------------------------------------------------        
C     THIS SUBROUTINE CALCULATES THE POTENTIAL TEMPERATURE REFERRED TO          
C     THE SURFACE. IT USES BRYDEN'S(1973) POLYNOMIAL FORMULA.                   
C     **UNITS**                                                                 
C       IN-SITU TEMPERATURE(T): DEGREES CENTIGRADE                              
C       SALINITY(S): PERMIL                                                     
C       PRESSURE(P): DECIBARS                                                   
C       POTENTIAL TEMPERATURE(PT): DEGREES CENTIGRADE                           
C     ** REFERENCE**                                                            
C       H.L.BRYDEN(1973); NEW POLYNOMIALS FOR THERMAL EXPANSION, ADIA-          
C       BATIC TEMPERATURE GRADIENT AND POTENTIAL TEMPERATURE OF SEA             
C       WATER. DEEP-SEA RESEARCH, VOL.20, PP.401-408.                           
C-----------------------------------------------------------------------        
C                                                                               
      PT = T - (0.36504E-4+0.83198E-5*T-0.54065E-7*T**2+0.40274E-9*T**3)        
     1        *P                                                                
     2       - (0.17439E-5-0.29778E-7*T)*P*(S-35.0)                             
     3       + 0.41057E-10*P**2*(S-35.0)                                        
     4       - (0.89309E-8-0.31628E-9*T+0.21987E-11*T**2)*P**2                  
     5       - (-0.16056E-12+0.50484E-14*T)*P**3                                
C                                                                               
      RETURN                                                                    
      END                                                                       
      SUBROUTINE SIGMA0(T,S,P,PD)                                               
C                                                                               
C-----------------------------------------------------------------------        
C     THIS SUBROUTINE CALCULATES THE POTENTIAL DENSITY OF A WATER PARCEL        
C     REFERRED TO THE SURFACE. IT USES BRYDEN'S(1973) FORMULA FOR SUR-          
C     FACE POTENTIAL TEMPERATURE AND UNESCO(1981) EQUATION OF STATE FOR         
C     SEA WATER.                                                                
C     **UNITS**                                                                 
C       IN-SITU TEMPERATURE(T): DEGREES CENTIGRADE                              
C       SALINITY(S): PERMIL                                                     
C       PRESSURE(P): DECIBARS                                                   
C       POTENTIAL DENSITY(PD): GIVEN AS DENSITY ANOMALY (POT.DENSITY            
C         - 1000.0 KG. PER CUBIC METER)                                         
C-----------------------------------------------------------------------        
C                                                                               
      CALL PTS(T,S,P,PT)                                                        
      CALL DIS(PT,S,0.0,RHO)                                                    
C                                                                               
      PD=RHO-1000.0                                                             
C                                                                               
      RETURN                                                                    
      END
      subroutine basinOCLsubALL(f,idim,jdim,kdim,h,filemask,
     * bOCL,numdeps)
      parameter (kdimloc=33) 
      dimension ibas(100),nstart(100)
      dimension f(idim,jdim,kdim), m1(idim,jdim)
      dimension h(kdim),hloc(kdimloc)
      character*80 filename, filemask, bOCL
      data nstart/ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     *             1, 1, 8,10,14,16,19,19,23,23,
     *            25,26,26,28,28,28,25,28,29,29,
     *            29,29,29,29,29,29,29,29,29,30,
     *            30,30,30,30,30,30,30,30,31,31,
     *            30,31, 1,14,14, 1, 6,32, 1, 1,
     *             1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     *             1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     *             1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     *             1, 1, 1, 1, 1, 1, 1, 1, 1, 1/
      data hloc/0.,10.,20.,30.,50.,75.,100.,125.,150.,200.,250.,
     *      300.,400.,500.,600.,700.,800.,900.,1000.,1100.,
     *      1200.,1300.,1400.,1500.,1750.,2000.,2500.,3000.,
     *      3500.,4000.,4500.,5000.,5500./

      call addCend(filemask,80,0)
      open(19,file=filemask,form='formatted',status='old')
      write(*,*) filemask
      do i=1,100
        read(19,*,end=500) ibas(i)
      enddo
 500  continue
      numbsns=i-1
      close(19)
      ibasinmask=-1
      call addCend(bOCL,80,0)
      call extraname('masks'//CHAR(0),bOCL,filename)
      call fileopen(filename,'rb+'//CHAR(0),ibasinmask)
      do k=1,numdeps
       call readleveli(1,1,1,1,1,idim,jdim,idim,jdim,idim,jdim,
     *                 k,k,ibasinmask,m1)
       do i=1,idim
       do j=1,jdim
          f(i,j,k)=float(m1(i,j))
       enddo
       enddo
      enddo
      call onefileclose(ibasinmask)
     
      do ii=1,numbsns
      if(ibas(ii).eq.0)then
c       kk=nstart(ii)
        xk=hloc(nstart(ii))+0.1
        do 75 k2=1,numdeps
75       if (h(k2) .le. xk ) kk=k2
        do k=kk,numdeps
        do j=1,jdim
        do i=1,idim
         ifijk=int(f(i,j,k))
         if(ifijk.eq.ii .or. ifijk.gt.numbsns) f(i,j,k)=-100.
        enddo
        enddo
        enddo
      endif
      enddo
      return
      end
       subroutine basin_se(pot,n,xin,win,now,nrms,tmean,semiss)
C 06/11/01
C INPUT ARRAYS should be without MISSING values
C
C 03/27/00 if POT<=0.0 TMEAN is UNTRIMMED mean
C
c computes standard error of basin
c trimmed/untrimmed mean/total for input unsorted data array
c
c      n - number of input values
c      x(n) - array of ODSQ standard errors (ODSQse)
c      w(n) - array of corresponding weights (area,volume, etc.)
c      pot  - proportion of data trimmed off each end of x*w
c             or (x*w)**2
c             e.g., pot=5 means 5% of the smallest and 
c             5% of the biggest input values will not be used
c             to compute the trimmed mean 
c      now  - if 1 then w=1 (for unweighted basin total)
c                2 then w=1/n (for unweighted basin mean)
c                3 then w are read-in weights
c      nrms - if 1 then se of basin mean is sum(w*ODSQse)
c                2 then se of basin mean is
c                                   sqrt(sum((w*ODSQse)**2))

       dimension xin(*),win(*)
       double precision s1,s2,w1,w2,sum,wsum,vk
       double precision x(64800),w(64800)
       integer uend

       if(now.eq.1)then
         do i=1,n
          w(i)=1.d0
         enddo
       elseif(now.eq.2)then
         do i=1,n
          w(i)=1.d0/real(n)
         enddo
       else
         do i=1,n
          w(i)=dble(win(i))
         enddo
       endif

       if(nrms.eq.1)then
         do i=1,n
          x(i)=dble(xin(i)) * w(i)
         enddo
       elseif(nrms.eq.2)then
         do i=1,n
          x(i)=dble(xin(i)) * w(i) * dble(xin(i)) * w(i)
         enddo
       endif

       if(pot.le.0.0)then
          s1=0.d0
          s2=0.d0
          w1=0.d0
          w2=0.d0
          goto 100
       endif

       p=0.01*pot

       lend=nint(p*n)
       call findk_se(lend,x,w,n,vk)
       s1=0.d0
       w1=0.d0
       do i=1,lend
         s1=s1+x(i)
         w1=w1+w(i)
       enddo

       uend=nint((1.-p)*n)
       call findk_se(uend,x,w,n,vk)
       s2=0.d0
       w2=0.d0
       do i=uend+1,n
         s2=s2+x(i)
         w2=w2+w(i)
       enddo

 100   sum=0.d0
       wsum=0.d0
       do i=1,n
        sum=sum+x(i)
        wsum=wsum+w(i)
       enddo

        sum=sum-s1-s2
        wsum=wsum-w1-w2
        if (wsum .gt. 0.d0)then
           if(now.le.2) wsum=1.d0
           if(nrms.eq.1)then
              tmean = real( sum / wsum )
           elseif(nrms.eq.2)then
              tmean = real( dsqrt(sum) / wsum )
           else
              tmean = semiss
           endif
        else
           tmean = semiss
        endif

       return
       end
************************************************************
       subroutine findk_se(k,a,a1,n,vkth)
       double precision a(*),a1(*),xx,ww,vkth
       integer r

       l=1
       r=n

       do 1000 while(l .lt. r)
        xx=a(k)
        i=l
        j=r
        do 1010 irept=1,1000000
         do while(a(i) .lt. xx)
         i=i+1
         enddo
         do while(xx .lt. a(j))
         j=j-1
         enddo
         if(i .le. j)then
          ww=a(i)
          a(i)=a(j)
          a(j)=ww
 
          ww=a1(i)
          a1(i)=a1(j)
          a1(j)=ww
 
          i=i+1
          j=j-1
         endif
         if(i .gt. j) goto 2000
1010  continue
2000  if(j .lt. k) l=i
      if(k .lt. i) r=j
1000  continue
      vkth=a(k)
      return
      end
      subroutine boxintegral01(f,n1,n2,bmiss,fout,fw,h,kdim)
C revised 06/05/02 see  4th note
C revised 06/01/01
C revised 12/10/00
C JIANT, 04/12/00
C
C Program computes vertical integral of standard depth (SD) level data.
C Data at each SD level weighted with half thickness of the layer defined
C as depth difference between previous and next SD levels except
C first (n1) and last (n2) depth data have a weight of half distance
C to the corresponding nearest SD level within n1,n2 range.
C
C NOTES:
C(1) If value at 1st depth is missing then output is BMISS
C
C(2) Integral will be evaluated from 1st depth to last depth with data
C    which could be less than h(n2)
C
C(3) If value at 2nd depth is missing then
C    output are fw=0.5*(h(n1+1)-h(n1)) and fout=f(n1)*fw
C
C(4) For n1=n2: (06/05/02) fout=f(n1)*0.5*(h(n1+1)-h(n1-1))
C
      dimension f(kdim),h(kdim)

      nsum=0
      fout=bmiss
      fw=bmiss
      hmax=h(kdim)+(h(kdim)-h(kdim-1))

      fval=f(n1)
      if(fval .gt. bmiss)then
         if ( n1 .lt. kdim ) then
          wh=0.5*(h(n1+1)-h(n1))
         else
          wh=0.5*(hmax-h(n1))
         endif
         sum=wh*fval
         wsum=wh
         nsum=nsum+1
      else
         fout=bmiss
         fw=bmiss
         return
      endif

      if(n1.eq.n2)then
       if(n1.gt.1)then
         if ( n1 .lt. kdim ) then
          fw=0.5*(h(n1+1)-h(n1-1))
         else
          fw=0.5*(hmax-h(n1-1))
         endif
         fout=fval*fw
         return
       else
         fout=sum
         fw=wsum
         return
       endif
      endif

      do k=n1+1,n2-1
         if ( k .le. kdim ) then
          fval=f(k)
         else
          fval=f(kdim)
         endif
         if(fval .gt. bmiss)then
          if ( k .lt. kdim ) then
	   wh= 0.5*(h(k+1)-h(k-1)) 
          else
	   wh= 0.5*(hmax-h(k-1)) 
          endif
	  sum=sum + wh*fval 
          wsum=wsum + wh
          nsum=nsum+1
         elseif(nsum.eq.1)then
          fout=sum
          fw=wsum
          return
         else
          if ( k .lt. kdim ) then
           wh= 0.5*(h(k)-h(k-1))
          else
           wh= 0.5*(hmax-h(k-1))
          endif
          fout=sum - wh*f(k-1)
          fw=wsum - wh
          return
         endif
      enddo

      fval=f(n2)
      wh= 0.5*(h(n2)-h(n2-1))
      if(fval .gt. bmiss)then
        fout=sum + wh*fval
        fw=wsum + wh
      else
        fout=sum - wh*f(n2-1)
        fw=wsum - wh
      endif

      return
      end
      subroutine boxintegralSE(f,n1,n2,bmiss,fout,fw)
C
C 06/18/07 updated for n1=n2 > 1 (to avoid negative FOUT)
C revised 06/05/02 see 4th note
C based on sub_depth_integral.f 
C JIANT, 06/04/01
C
C Computes Root_Mean_Square errors of vertically integrated data. 
C
C  Definition of vertical integration:
C Data at each SD level weighted with half thickness of the layer defined
C as depth difference between previous and next SD levels except
C first (n1) and last (n2) depth data have a weight of half distance
C to the corresponding nearest SD level within n1,n2 range.
C
C NOTES:
C(1) If value at 1st depth is missing then output is BMISS
C
C(2) Integral will be evaluated from 1st depth to last depth with data
C    which could be less than h(n2)
C
C(3) If value at 2nd depth is missing then
C    output are fw=0.5*(h(n1+1)-h(n1)) and fout=f(n1)*fw
C
C(4) For n1=n2: (06/05/02) see sub_depth_integral.f
C
      dimension f(33)
      double precision dpp,wh,sum,wsum,h(34)
      data h/0.d0,10.d0,20.d0,30.d0,50.d0,75.d0,100.d0,
     *      125.d0,150.d0,200.d0,250.d0,300.d0,400.d0,500.d0,
     *      600.d0,700.d0,800.d0,
     *      900.d0,1000.d0,1100.d0,1200.d0,1300.d0,1400.d0,
     *      1500.d0,1750.d0,2000.d0,2500.d0,3000.d0,
     *      3500.d0,4000.d0,4500.d0,5000.d0,5500.d0,6000.d0/

      nsum=0
      fout=bmiss
      fw=bmiss

      fval=f(n1)
      if(fval .gt. bmiss)then
         wh=0.5d0*(h(n1+1)-h(n1))
         sum=wh*wh*dble(fval*fval)
         wsum=wh*wh
         nsum=nsum+1
      else
         fout=bmiss
         fw=bmiss
         return
      endif

      if(n1.eq.n2)then
       if(n1.gt.1)then
         fw= 0.5*(h(n1+1)-h(n1-1))  
         fout= fval*fw
C line below added on 06/18/07
         fout=sqrt(fout*fout)
         return
       else
         fout=real( dsqrt(sum) )
         fw=real( dsqrt(wsum) )
         return
       endif
      endif

      do k=n1+1,n2-1
         fval=f(k)
         if(fval .gt. bmiss)then
	  wh= 0.5d0*(h(k+1)-h(k-1)) 
	  sum=sum + wh*wh*dble(fval*fval) 
          wsum=wsum + wh*wh
          nsum=nsum+1
         elseif(nsum.eq.1)then
          fout=real( dsqrt(sum) )
          fw=real( dsqrt(wsum) )
          return
         else
          wh= 0.5d0*(h(k)-h(k-1))
          dpp=  sum - wh*wh*dble(f(k-1)*f(k-1))
          fout= real( dsqrt(dpp) )
          fw= real( dsqrt(wsum - wh*wh) )
          return
         endif
      enddo

      fval=f(n2)
      wh= 0.5d0*(h(n2)-h(n2-1))
      if(fval .gt. bmiss)then
        dpp= sum + wh*wh*dble(fval*fval) 
        fout= real( dsqrt(dpp) )
        fw= real( dsqrt(wsum + wh*wh) )
      else
        dpp= sum - wh*wh*dble(f(n2-1)*f(n2-1))
        fout= real( dsqrt(dpp) )
        fw= real( dsqrt(wsum - wh*wh) )
      endif

      return
      end
      subroutine janameit(n,infl,c1,outfl,c)
      character*1 c1
      character*(*) infl
      character*(*) outfl
      character*(*) c(*)
      ic=index(infl,c1)-1
      outfl=infl(1:ic)
      do i=1,n
         itrim=index(outfl,' ')-1
         outfl=outfl(1:itrim)//c(i)
      enddo
      return
      end
      subroutine geo2grid(ra,ro,ia,io,ca,co)
C converts geographical coordinates of the center of 1x1 square
C ( SD[DD].5 where S is + for EAST(NORTH) or - for WEST(SOUTH)
C                         INTO
C corresponding grid numbers (1 to 360 for LONG and 1 to 180
C for LAT) and 
C its character images (e.g., for outfile name construction)
      character*3 ca
      character*4 co

C convert into characters
C      LATITUDE
       if(ra.lt.0.)then
          ix=int(-1.*ra)
          ca(3:3)='s'
       else
          ix=int(ra)
          ca(3:3)='n'
       endif

       if(ix.le.9)then
          ca(1:1)='0'
          write(ca(2:2),'(i1)') ix
       else
          write(ca(1:2),'(i2)') ix
       endif
       
C      LONGITUDE
       if(ro.lt.0.)then
          ix=int(-1.*ro)
          co(4:4)='w'
       else
          ix=int(ro)
          co(4:4)='e'
       endif

       if(ix.le.9)then
          co(1:2)='00'
          write(co(3:3),'(i1)') ix
       elseif(ix.le.99)then
          co(1:1)='0'
          write(co(2:3),'(i2)') ix
       else
          write(co(1:3),'(i3)') ix
       endif

C convert into grid numbers

       ra=ra+90.5
       ia=int(ra)
       if(ro.le.0.)then
          io=nint(ro+360.)
       else
          io=nint(ro)
       endif
      
      return
      end
C*************************************************************
      subroutine ch12dep(cx)
      character*12 cx(33)
      cx(1) ='          0m'
      cx(2) ='         10m'
      cx(3) ='         20m'
      cx(4) ='         30m'
      cx(5) ='         50m'
      cx(6) ='         75m'
      cx(7) ='        100m'
      cx(8) ='        125m'
      cx(9) ='        150m'
      cx(10)='        200m'
      cx(11)='        250m'
      cx(12)='        300m'
      cx(13)='        400m'
      cx(14)='        500m'
      cx(15)='        600m'
      cx(16)='        700m'
      cx(17)='        800m'
      cx(18)='        900m'
      cx(19)='       1000m'
      cx(20)='       1100m'
      cx(21)='       1200m'
      cx(22)='       1300m'
      cx(23)='       1400m'
      cx(24)='       1500m'
      cx(25)='       1750m'
      cx(26)='       2000m'
      cx(27)='       2500m'
      cx(28)='       3000m'
      cx(29)='       3500m'
      cx(30)='       4000m'
      cx(31)='       4500m'
      cx(32)='       5000m'
      cx(33)='       5500m'
      return
      end
      subroutine rinit3D(f,idim,jdim,kdim,val)
      dimension f(idim,jdim,kdim)
        do k=1,kdim
        do j=1,jdim
        do i=1,idim
           f(i,j,k)=val
        enddo
        enddo
        enddo
      return
      end
      subroutine rinit2D(f,idim,jdim,val)
      dimension f(idim,jdim)
        do j=1,jdim
        do i=1,idim
           f(i,j)=val
        enddo
        enddo
      return
      end
      subroutine rinit1D(f,idim,val)
      dimension f(idim)
        do i=1,idim
           f(i)=val
        enddo
      return
      end
      subroutine init3D(f,idim,jdim,kdim,ival)
      integer f(idim,jdim,kdim)
        do k=1,kdim
        do j=1,jdim
        do i=1,idim
           f(i,j,k)=ival
        enddo
        enddo
        enddo
      return
      end
      subroutine init2D(f,idim,jdim,ival)
      integer f(idim,jdim)
        do j=1,jdim
        do i=1,idim
           f(i,j)=ival
        enddo
        enddo
      return
      end
      subroutine init1D(f,idim,ival)
      integer f(idim)
        do i=1,idim
           f(i)=ival
        enddo
      return
      end
      subroutine nameTIME(itime,iter,ct1,ct2)
      character*2 ct1,ct2
      character*1 az(9)
      data az/'A','B','C','D','E','F','G','H','I'/
       if(iter.le.9)then
        ct1(1:1)='0'
        write(ct1(2:2),'(i1)') iter
       elseif(iter.le.99)then
        write(ct1(1:2),'(i2)') iter
       else
        nMM=iter-100
        nRR=int(nMM/10)
        nMM=nMM-10*nRR
        ct1(1:1)=az(nRR+1)
        write(ct1(2:2),'(i1)') nMM
       endif

       inext=iter+itime
       if(inext.le.9)then
        ct2(1:1)='0'
        write(ct2(2:2),'(i1)') inext
       elseif(inext.le.99)then
        write(ct2(1:2),'(i2)') inext
       else
        nMM=inext-100
        nRR=int(nMM/10)
        nMM=nMM-10*nRR
        ct2(1:1)=az(nRR+1)
        write(ct2(2:2),'(i1)') nMM
       endif
      return
      end
      subroutine std_dev_mn(y,n,ax,sig)
      real y(*)
      ax=0.
      sx=0.
      do i=1,n
         ax=ax+y(i)
      enddo
      ax=ax/n
      do i=1,n
         x=y(i)-ax
         sx=sx+x**2
      enddo
      sig=sqrt(sx/(n-1))
      return
      end
       subroutine trimmer(pot,n,x,w,now,tmean,bmiss)
C 05/21/07 corrected computation of UEND
c 06/12/01 basin total computation is added
C 03/27/00 if POT<=0.0 TMEAN is UNTRIMMED mean
c computes trimmed mean for input unsorted data array
c      n - number of input values
c      x(n) - input data array
c      w(n) - array of corresponding weights
c      pot  - proportion of data trimmed off each end of x(n)
c             e.g., pot=5 means 5% of the smallest and 
c             5% of the biggest input values will not be used
c             to compute the trimmed mean 
c      now  - if 1 unweighted basin total
c                2 unweighted basin mean
c                3 weighted basin mean (read-in weights)

       dimension x(*),w(*)
       integer uend

       if(now.le.2)then
         do i=1,n
          w(i)=1.
         enddo
       endif

       do i=1,n
         x(i)=x(i)*w(i)
       enddo

       if(pot.le.0.0)then
          s1=0.
          s2=0.
          w1=0.
          w2=0.
          goto 100
       endif

       p=0.01*pot

       lend=nint(p*n)
       call findk2(lend,x,w,n,vk)
       s1=0.
       w1=0.
       do i=1,lend
         s1=s1+x(i)
         w1=w1+w(i)
       enddo

c      uend=nint((1.-p)*n)
       uend=n-lend+1
       call findk2(uend,x,w,n,vk)
       s2=0.
       w2=0.
c      do i=uend+1,n
       do i=uend,n
         s2=s2+x(i)
         w2=w2+w(i)
       enddo

 100   sum=0.
       wsum=0.
       do i=1,n
        sum=sum+x(i)
        wsum=wsum+w(i)
       enddo

        sum=sum-s1-s2
        wsum=wsum-w1-w2
        if (wsum .ne. 0.)then
          if(now.eq.1)then
           tmean=sum
          else
           tmean=sum/wsum
          endif
        else
           tmean=bmiss
        endif

       return
       end
************************************************************
       subroutine findk2(k,a,a1,n,vkth)
       dimension a(*),a1(*)
       integer r

       l=1
       r=n

       do 1000 while(l .lt. r)
        x=a(k)
        i=l
        j=r
        do 1010 irept=1,1000000
         do while(a(i) .lt. x)
         i=i+1
         enddo
         do while(x .lt. a(j))
         j=j-1
         enddo
         if(i .le. j)then
          w=a(i)
          a(i)=a(j)
          a(j)=w
c
          w=a1(i)
          a1(i)=a1(j)
          a1(j)=w
c
          i=i+1
          j=j-1
         endif
         if(i .gt. j) goto 2000
1010  continue
2000  if(j .lt. k) l=i
      if(k .lt. i) r=j
1000  continue
      vkth=a(k)
      return
      end
