getsigprop.f
by
Webmaster Legos
—
last modified
Jul 05, 2012 12:34 PM
getsigprop.f
—
Fortran source code,
2 kB (2,704 bytes)
File contents
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