program PWREFIND C Calculating averages integer RD,RD2,rns,runs,igit,nb,kgit,nvar PARAMETER (runs=10) PARAMETER (igit=64) PARAMETER (nb=36) parameter (kgit=56) parameter (nvar=3) parameter(ntime=1201) real uvt(igit,nb,kgit,nvar),wvel(igit,nb,kgit) integer prog CHARACTER*1 jfilein(10) c===========resid computing============================= real pi,a,H,R,cp,zoben,zunten,dz,deltf data a /6366197./ parameter(H=7000.,R=287.04,cp=1005.) parameter(zoben=157700.) real z(kgit),rgit,rb,fi,rtime real cosfi(nb),sinfi(nb),defi,omega c---------------------------------------- zonally and time averaged values real uq(nb,kgit),vq(nb,kgit),tq(nb,kgit),wq(nb,kgit) c-----------------------perturbation from zonally and time averaged values real dvq(igit,nb,kgit),dtq(igit,nb,kgit),tqdz(nb,kgit) c---------------------------------------- zonally averaged nonlinear terms real vstpsq(nb,kgit),tq3(igit,nb,kgit) real vres(nb,kgit),wres(nb,kgit) c------------------------- residual meridional circulation and derivatives real tpqzz,vstpsdz,deriv real deriw,vstpsfi,tpqdzfi pi =2.*asin(1.) dz =zoben/(float(kgit)-0.5) zunten=dz/2. rgit =float(igit) rtime =float(ntime) rb =float(nb) defi =180./rb omega =2.*pi / (24.*3600.) do k=1,kgit z(k)=zunten+(k-1)*dz enddo c=========================input file & computing area==================================== DATA jfilein/'0','1','2','3','4','5' &,'6','7','8','9'/ c---------------- geometric factors---------------------------------- do j=1,nb fi =(float(j)-18.5)*defi*pi/180. ! -87.5,...,87.5 cosfi(j)=cos(fi) sinfi(j)=sin(fi) enddo c----------------- NEUDEFINITION VON defi ==> deltf deltf=defi*pi/180. c================== ALL RUNS ================================= DO rns = 1,runs RD = 21 RD2 = 22 c================== DATA OUTPUT ================================= open(16,FILE="../../Muam_Mar_wQBO/" //jfilein(rns)// "/resid.dx", $ access='direct',status='unknown', $ recl=4*nb*kgit*2) c================= DATA INPUT ================================ OPEN(RD,FILE= *'../../Muam_Mar_wQBO/' //jfilein(rns)// '/uvt400_Mar_wQBO.dx', * form='unformatted', * access='direct', status='old',recl=nvar*nb*igit*kgit*4) OPEN(RD2,FILE= "../../Muam_Mar_wQBO/" //jfilein(rns)// "/wvel.dx", $ form='unformatted', $ access='direct', status='old',recl=nb*igit*kgit*4) c=============================================================== c================== time loop ================================= DO prog = 1,ntime print*,"run: ", rns," step: ", prog c================== reset all================================= do k=1,kgit do j=1,nb uq(j,k)=0. vq(j,k)=0. tq(j,k)=0. wq(j,k)=0. wres (j,k)=0. vres (j,k)=0. vstpsq(j,k)=0. tqdz(j,k)=0. do i=1,igit dvq(i,j,k)=0. dtq(i,j,k)=0. uvt(i,j,k,1)=0. uvt(i,j,k,2)=0. uvt(i,j,k,2)=0. wvel(i,j,k)=0. tq3(i,j,k)=0. enddo enddo enddo deriv=0;deriw=0;vstpsfi=0;tpqdzfi=0;tpqzz=0;vstpsdz=0 read(RD,rec=prog) uvt read(RD2,rec=prog) wvel c-------averaging----------------- do k=1,kgit do j=1,nb do i=1,igit c-----------teta (potential temperature) ---------- tq3(i,j,k) = uvt (i,j,k,3)*exp(R/cp/H*z(k)) c---------zonal-mean input tq(j,k)=tq(j,k)+tq3(i,j,k) /rgit uq(j,k)=uq(j,k)+uvt (i,j,k,1) /rgit vq(j,k)=vq(j,k)+uvt (i,j,k,2) /rgit wq(j,k)=wq(j,k)+wvel(i,j,k ) /rgit enddo !print*, j,k,uq(j,k),tq(j,k) enddo enddo c--------perturbation from the zonal-mean data do k=1,kgit do j=1,nb do i=1,igit dvq(i,j,k)=uvt(i,j,k,2)-vq(j,k) dtq(i,j,k)=tq3(i,j,k)-tq(j,k) vstpsq(j,k)= vstpsq(j,k)+dvq(i,j,k)*dtq(i,j,k)/rgit enddo enddo enddo c=========CALCULATING RESIDUAL:=================================== c----------d/dz(teta), ---------------------------------- c do k=1,kgit do j=1,nb if(k.ne.1.and.k.ne.kgit) then tqdz(j,k)= (tq(j,k+1)-tq(j,k-1))/2./dz ! d ttq/dz else if(k.eq.1) then tqdz(j,k)=(tq(j,k+1)-tq(j,k)) / dz else if(k.eq.kgit) then tqdz(j,k)=(tq(j,k)-tq(j,k-1)) / dz endif c enddo enddo c-------------------------------------------------------------------- do k=1,kgit do j=1,nb if(k.ne.1.and.k.ne.kgit) then c--------------d2/dz2(teta), d/dz(V*teta)----------------------------------- tpqzz =(tqdz(j,k+1)-tqdz(j,k-1))/2./dz vstpsdz =(vstpsq(j,k+1)-vstpsq(j,k-1))/2./dz else if(k.eq.1) then tpqzz =(tqdz(j,k+1)-tqdz(j,k))/dz vstpsdz =(vstpsq(j,k+1)-vstpsq(j,k))/dz else if(k.eq.kgit) then tpqzz =(tqdz(j,k)-tqdz(j,k-1))/dz vstpsdz=(vstpsq(j,k)-vstpsq(j,k-1))/dz endif c-------------------------------------------------------------------- if(j.eq.1) then c--------------d/dphi(V*teta), d2/dz*dphi(teta), ----------------------------------- vstpsfi=(vstpsq(j+1,k)-vstpsq(j,k))/deltf tpqdzfi=(tqdz(j+1,k)-tqdz(j,k))/deltf c------------------------------------------------------------------- else if(j.eq.nb) then c------------------------------------------------------------------- vstpsfi=(vstpsq(j,k)-vstpsq(j-1,k))/deltf tpqdzfi=(tqdz(j,k)-tqdz(j-1,k))/deltf c------------------------------------------------------------------- else c------------------------------------------------------------------- vstpsfi=(vstpsq(j+1,k)-vstpsq(j-1,k))/2./deltf tpqdzfi=(tqdz (j+1,k)-tqdz (j-1,k))/2./deltf c------------------------------------------------------------------- endif c-----Calculate the residual mean meridional circulation------------ c--------v*=-1/rho0(rho0*/_z)_z ----Andrews, seite 128 c c--------w*=+1/(a*cosfi)*(cosfi*/_z)_fi--------------- if(wq(j,k).ne.0.and.tqdz(j,k).ne.0.and.uq(j,k).ne.0) then deriv= 1./tqdz(j,k)*(-vstpsq(j,k)/H+vstpsdz & -vstpsq(j,k)*tpqzz/tqdz(j,k)) vres(j,k)=vq(j,k)-deriv deriw= 1./a/cosfi(j)/tqdz(j,k)*(-sinfi(j)*vstpsq(j,k) & +cosfi(j)*(vstpsfi-tpqdzfi*vstpsq(j,k)/tqdz(j,k))) wres(j,k)=wq(j,k)+deriw else vres(j,k) = 0 wres(j,k) = 0 endif enddo enddo c=========END CALCULATING RESIDUAL==================================== write(16,rec=prog) vres,wres END DO ! time loop c============close all============================== close(RD) close(RD2) close(16) ENDDO ! runs STOP END *************************************************************************