program elength ! calculates equivalent length from tracer data !============================================================================== ! This code may be freely distributed and altered as required. ! For details of equivalent length and its calculation see: ! Haynes and Shuckburgh 2000a,b, J. Geophys. Res, 105, 22777-22810. ! Expression calculated is $\kappa_{\rm eff}=f/g$, ! where $f=\frac{1}{r^2\cos(\phi_e)\frac{\partial}{\partial\phi_e}\int_{c\le C}|\nabla C|^2dA$, ! and where $g=(\frac{1}{r^2\cos(\phi_e)}\frac{\partial C}{\partial\phi_e})^2$ !============================================================================== !!Section 1: Initialisation !SETTING ARRAY DIMENSIONS !++++++++++++++++++++++++ real*8 pi !tracer data !----------- integer*8 mm, jg, mg, jg2, n ! gaussian latitude grid, T"mm", "jg2"=#lats, "mg"=#lons, "n=#grid-pts integer*8 tdim ! "tdim="#days real, allocatable, dimension(:,:) :: tr, tra, sr real*4, allocatable, dimension(:) :: atr ! "tr"=tracer data c ! "atr"=zonal-mean tracer data !equivalent latitude grid !------------------------ real*8, allocatable, dimension(:,:) :: arlat ! "arlat"=area of segment of band either side a latitude real*8, allocatable, dimension(:) :: da, aq, arcap ! "arcap"=area in half-latitude cap ! "da"=arlat ! "aq"=area in cap defined by gaussian latitude !equivalent length results !------------------------- integer*8 na ! "na"=#equiv-lats integer*8 iaq, i, iskip, incr, incr2 ! "iaq", "i", "iskip", "incr", "incr2" = indexes integer*8 ilataN, ilataS, nalats, nn ! "[ilataN,ilataS]"=equiv-lat-band of calculation, "nalats"=#lats in band, "nn"=#grid-pts in band real*8, allocatable, dimension(:,:) :: ltracer, dtrdlo ! "ltracer"=$|\nabla c|^2=\chi$, "dtrdlo"=differentiated c real*8, allocatable, dimension(:) :: q, qa, chi, chia ! "q" is the tracer value ! "qa" is the tracer value equivalent to area aq ! "chia" is the chi value equivalent to area aq real*8 area real*8 chiint ! "chiint"=$\int\chi dA$, where $chi=|\nabla c|^2$ real*8, allocatable, dimension(:) :: chihat ! "chihat"=area-averaged quantity real*8, allocatable, dimension(:) :: le, mle ! "le"=equivalent length, L_{\rm eq} integer nband ! "nband"=no of grid boxes to mark as anomalous due to tracer contours being cut off at edge of band !index variables !--------------- integer*8, allocatable, dimension(:) :: index,inos ! "index" = indices of the ascending tracer values ! "inos" = the number of bits of area added between area markers !grid dependent variables !------------------------ real*8, allocatable, dimension(:) :: lats,lons real*4, allocatable, dimension(:) :: plons,plats integer*8, allocatable, dimension(:) :: jgval,mgval,iflag integer*2 nplon,nplat,idd,imm,iyy ! "lats" pole -> pole ! "jgval" = latitude value of ith n value ! "mgval" = longitude value of ith n value !SETTING PARAMETERS !++++++++++++++++++ pi=4.d0*datan(1.d0) mm=170 jg=(3*mm+4)/4 mg=4*jg jg2=jg*2 tdim=2190 na=jg2 n=jg2*mg nband=7 !INPUT FILES !+++++++++++ open(11,file='../T170-DATA/o3-170-mls.dat',form='unformatted') ! tracer data open(12,file='latitudes',form='formatted') ! gaussian latitudes open(15,file='../N2O/mls-o3.hrzm',form='formatted') ! zonal-mean of smoothed tracer data !OUTPUT FILES !++++++++++++ open(13,file='L-mls-170.dat',form='formatted') !ALLOCATE ARRAYS !--------------- !tracer data !----------- allocate(tr(mg,jg2),tra(mg,jg*2),atr(jg2)) !equivalent latitude grid !------------------------ allocate(arlat(mg,jg2),da(n),aq(na+1),arcap(jg2)) !equivalent length results !------------------------- allocate(mle(jg2)) !allocate grid dependent variables !--------------------------------- allocate(lons(mg),lats(jg2),iflag(jg2)) allocate(plons(mg),plats(jg2)) !SET UP GRID !+++++++++++ !lats and lons (lats are Nth(90) -> Eq(0) -> Sth(-90)) !----------------------------------------------------- write(14,*)mm write(14,*)tdiss write(6,*)'mm=',mm read(12,*)lats lats(:)=dble(pi*lats(:)/180.d0) close(12) write(6,*)'latitudes=',lats do i=1,mg lons(i)=dble(2.d0*pi*(i-1)/real(mg)) enddo !to calculate area in band +/- half-latitude circle, and in latitude cap !----------------------------------------------------------------------- do i=1,jg2-1 arcap(i)=dble(2.d0*pi*(1.d0-dsin((lats(i)+lats(i+1))/2.d0))) aq(i)=dble(2.d0*pi*(1.d0-sin(lats(i)))) enddo arcap(jg2)=dble(2.d0*pi*(1.d0-sin((lats(jg2)-pi/2.d0)/2.d0))) aq(jg2)=dble(2.d0*pi*(1.d0-sin(lats(jg2)))) aq(jg2+1)=4.d0*pi+0.001 do j=1,mg arlat(j,1)=dble((arcap(1))/real(mg)) do i=2,jg2-1 arlat(j,i)=dble((arcap(i)-arcap(i-1))/real(mg)) enddo arlat(j,jg2)=dble((4.d0*pi-arcap(jg2-1))/real(mg)) enddo do j=1,jg2 da(1+mg*(j-1):mg*j)=arlat(:,j) enddo !! Section 2: Time step over days do iday=1,tdim !to read in tracer !----------------- read(11)tra write(6,*)'day=',iday iflag(:)=0 read(15,*)atr do ijg=1,jg2 if(atr(ijg)>1e25)then iflag(ijg)=1 endif enddo ilataN=1 ilataS=2 ilast=0 mle(:)=-1.d0 ! work out band in which tracer is monotonic !------------------------------------------- do while(ilast==0) if((iflag(ilataN)==1).or.(iflag(ilataN+1)==1))then mle(ilataN)=-1.d0 ilataN=ilataN+1 if(ilataN.ge.jg2)then ilast=1 endif else ilataS=ilataN+1 if(atr(ilataN) arcap(ilataS)) then area=arcap(ilataS) endif chiint=dble(chiint+chi(index(i)-(ilataN-1)*mg)* $ da(index(i))) iskip=iskip+1 !iskip =number of bits added=1,2,3... if (area > aq(iaq)) then qa(iaq-ilataN+1)=q(index(i)-(ilataN-1)*mg) inos(iaq-ilataN+1)=iskip !inos = number of bits added=1,2,3... iskip=0 chia(iaq-ilataN+1)=chiint iaq=iaq+1 endif enddo !rest of the tracer qa(nalats)=q(index(nn)-(ilataN-1)*mg) chia(nalats)=chiint inos(nalats)=0 inos(nalats)=nn-sum(inos,1) !to get Nakamura's hat, need to differentiate !-------------------------------------------- call derivd(lats(ilataN:ilataS),chia,nalats,chihat) chihat(:)=chihat(:)/cos(lats(ilataN:ilataS)) !Nakamura's calculation !---------------------- call derivd(lats(ilataN:ilataS),qa,nalats,le) le(:)=-(dble(chihat(:)*(cos(lats(ilataN:ilataS)))**2/ $ (le(:))**2))/(2.d0*pi* $ (cos(lats(ilataN:ilataS)))**2) mle(ilataN:ilataS)=le(:) !eliminate dodgy results !----------------------- mle(ilataN:ilataN-1+nband)=0.1d0 mle(ilataS+1-nband:ilataS)=0.1d0 !deallocate arrays !----------------- deallocate(chihat, le) deallocate(qa, chia, inos) deallocate(index) deallocate(dtrdlo, sr, q, chi) deallocate(ltracer) ilataN=ilataS+1 if(ilataN.ge.jg2)then ilast=1 endif endif if(ilataN.ge.jg2)then ilast=1 endif endif if(ilataN.ge.jg2)then ilast=1 endif enddo write(13,100)mle 100 format(8e13.5) enddo close(13) end !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- subroutine deriv(x,y,ndim,z) real y(ndim,1) real*8 z(ndim,1) real*8 x(ndim) do k=2,ndim-1 z(k,1) = dble(y(k-1,1) - y(k+1,1))/(x(k-1) - x(k+1)) enddo z(1,1) = dble(y(2,1) - y(1,1))/(x(2)-x(1)) z(ndim,1) = dble(y(ndim,1)-y(ndim-1,1))/(x(ndim)-x(ndim-1)) return end !============================================================================= subroutine derivd(xs,ys,ndims,zs) real*8 ys(ndims) real*8 zs(ndims) real*8 xs(ndims) do k=2,ndims-1 zs(k) = dble(ys(k-1) - ys(k+1))/(xs(k-1) - xs(k+1)) enddo zs(1) = dble(ys(2) - ys(1))/(xs(2)-xs(1)) zs(ndims) =dble(ys(ndims)-ys(ndims-1))/(xs(ndims)-xs(ndims-1)) return end !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- subroutine sort(n,arrin,indx) integer*8 n real*8 arrin(n) integer*8 indx(n), indxt, ir, l, i, j !initial index goes 1,2,3... do j=1,n indx(j)=j enddo l=n/2+1 ir=n do if(l > 1)then l=l-1 indxt=indx(l) q=arrin(indxt) else indxt=indx(ir) q=arrin(indxt) indx(ir)=indx(1) ir=ir-1 if(ir == 1)then indx(1)=indxt return endif endif i=l j=l+l do while(j<=ir) if(j