Vous êtes ici : Accueil / Alexandre Ganachaud / hydrosys / getsigprop.f

getsigprop.f

Par Webmaster Legos Dernière modification 05/07/2012 12:34

Fortran source code icon getsigprop.f — Fortran source code, 2 kB (2704 bytes)

Contenu du fichier

subroutine getsigprop(nd,np,dep,prop,botdep,nl,sigdep,sigprop)
c
c synopsis :
 
c  depths must be positive
c

c description : 
 
c   interpolates linearily the property prop to values along at
c   certain depths (following an isopycnal typically)
c   given by sigdep
c
c INPUTS
c   dep(nd)     :standard depths for the data
c   prop(nd,np)  :property to interpolate, for each pair
c   botdep(np)  :bottom depth for each pair
c
c   sigdep(np,nl):depth of the isopycnal for each pair, and
c                 for each isopycnal layer
c
c OUTPUTS (I/O)
c   sigprop(np,nl):interpolated property
c

c uses :

c side effects : if we are outside the data range 
c                sigprop is linearily extrapolated
c                if under bottom, sigprop is the one at the bottom.
c author : A.Ganachaud, Nov 96

c see also :

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      integer*4 nd,np,nl
      real*8    dep(nd),prop(nd,np),botdep(np)
      real*8    sigdep(np,nl),sigprop(np,nl)

c     indices: id = depth indices
c              ip = pair indices
c              il = sigma-layer indices


c     FOR EACH LAYER
      do il=1,nl
         
c        FOR EACH PAIR
         do ip=1,np
            xd=sigdep(ip,il)
c           FIND SURROUNDING STANDARD DEPTHS
            do id=1,nd
               if (xd.le.dep(id)) then
                  goto 100
               endif
            enddo
 100        continue
            if (id.eq.1) then   
c              sigma at the surface
               if (dep(2).lt.botdep(ip)) then
                  sigprop(ip,il)=prop(1,ip) + (xd-dep(1))*
     &                 (prop(2,ip)-prop(1,ip))/(dep(2)-dep(1))

               else 
c                 THERE IS ONLY ONE POINT
                  sigprop(ip,il)=prop(1,ip)
               endif

            elseif (dep(id).gt.botdep(ip)) then
               if (xd.gt.botdep(ip)) then 
c                 UNDERGROUND CASE
C                 THE SIGPROP TAKES THE EXTRAPOLATED VALUE 
C                 AT THE BOTTOM
                  xd=botdep(ip)
               endif
C              EXTRAPOLATION
               sigprop(ip,il)=prop(id-1,ip) + (xd-dep(id-1))*
     &              (prop(id-1,ip)-prop(id-2,ip))/(dep(id-1)-dep(id-2))

            else                
c              REGULAR INTERPOLATION
               sigprop(ip,il)=(prop(id-1,ip)*(dep(id)-xd) -
     &              prop(id,ip)*(dep(id-1)-xd))/(dep(id)-dep(id-1))
               
            endif

         enddo                  
c        LOOP ON EACH PAIR

      enddo                     
c     LOOP ON EACH LAYER

      return
      end

Actions sur le document

logo cnes logo IRD Logo université de Toulouse Logo université Paul Sabatier Logo CNRS
Logo bibliothèque OBS Logo Observatoire Midi Pyrénées