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=373,nstart=361) 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.,dich0=1.273) real z(kgit),rho0(kgit),rgit,rb,fi,rtime,dich0 real cosfi(nb),f(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) real duq(igit,nb,kgit),dwq(igit,nb,kgit) c--------------------------------- derivatives of zonally averaged values real uqdz(nb,kgit),uqdfi(nb,kgit) c---------------------------------------- zonally averaged nonlinear terms real usvsq(nb,kgit),uswsq (nb,kgit) 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 c-----------EP-flux and divrgence real gradF(nb,kgit) real Fphi(nb,kgit),Fz(nb,kgit),dFphi(nb,kgit),dFz(nb,kgit) 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 rho0(k)=dich0*exp(-z(k)/H) enddo c=========================input file & computing area==================================== DATA jfilein/'e','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) f(j) =2.*omega*sinfi(j) c print*, j,fi,cosfi(j),sinfi(j) enddo c----------------- NEUDEFINITION VON defi ==> deltf deltf=defi*pi/180. c================== ALL RUNS ================================= DO rns = 1,1 !runs ! runs RD = 21 RD2 = 22 prog =1 open(24,file='test.dat') c================== DATA OUTPUT ================================= open(16,FILE="../../Muam_Mar_eQBO/" //jfilein(rns)// "/resEPa.dx", $ access='direct',status='unknown', $ recl=4*nb*kgit*6) c open(16,FILE="../../EL2_OGW_volf10/EP_IR/resEP.dx", c $ access='direct',status='unknown', c $ recl=4*nb*kgit*6) c================= DATA INPUT ================================ OPEN(RD,FILE= *'../../Muam_Mar_eQBO/' //jfilein(rns)// '/uvt400_Mar_eQBO.dx', * form='unformatted', * access='direct', status='old',recl=nvar*nb*igit*kgit*4) OPEN(RD2,FILE= "../../Muam_Mar_eQBO/" //jfilein(rns)// "/wvel.dx", $ form='unformatted', $ access='direct', status='old',recl=nb*igit*kgit*4) c OPEN(RD,FILE= c *'../../EL2_OGW_volf10/uvt400_Jan_EL_OGW_volf.dx', c * form='unformatted', c * access='direct', status='old',recl=nvar*nb*igit*kgit*4) c OPEN(RD2,FILE= "../../EL2_OGW_volf10/wvel.dx", c $ form='unformatted', c $ access='direct', status='old',recl=nb*igit*kgit*4) c=============================================================== c================== time loop ================================= 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. usvsq(j,k)=0; uswsq(j,k)=0; vstpsq(j,k)=0 tqdz(j,k)=0; uqdz(j,k)=0 uqdfi(j,k)=0 do i=1,igit duq(i,j,k)=0; dvq(i,j,k)=0; dwq(i,j,k)=0 uvt(i,j,k,1)=0; uvt(i,j,k,2)=0; uvt(i,j,k,3)=0. dtq(i,j,k)=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 DO prog = nstart,nstart+ntime 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/rtime uq(j,k)=uq(j,k)+uvt (i,j,k,1) /rgit/rtime vq(j,k)=vq(j,k)+uvt (i,j,k,2) /rgit/rtime wq(j,k)=wq(j,k)+wvel(i,j,k ) /rgit/rtime enddo !print*, j,k,uq(j,k),tq(j,k) enddo enddo 11 format(i5,2f10.3) c--------perturbation from the zonal-mean data do k=1,kgit do j=1,nb do i=1,igit duq(i,j,k)=duq(i,j,k)+(uvt(i,j,k,1)-uq(j,k))/rtime dvq(i,j,k)=dvq(i,j,k)+(uvt(i,j,k,2)-vq(j,k))/rtime dwq(i,j,k)=dwq(i,j,k)+(wvel(i,j,k )-wq(j,k))/rtime dtq(i,j,k)=dtq(i,j,k)+(tq3(i,j,k)-tq(j,k))/rtime enddo enddo enddo ENDDO do k=1,kgit do j=1,nb do i=1,igit usvsq(j,k)= usvsq(j,k)+duq(i,j,k)*dvq(i,j,k)/rgit uswsq(j,k)= uswsq(j,k)+duq(i,j,k)*dwq(i,j,k)/rgit vstpsq(j,k)= vstpsq(j,k)+dvq(i,j,k)*dtq(i,j,k)/rgit enddo enddo enddo print*, 0, usvsq(36,20) c=========CALCULATING RESIDUAL:=================================== c----------d/dz(teta), ---------------------------------- c do k=1,kgit do j=1,nb c----------------Correction for polar regions-------------------------- if(j.eq.1) then usvsq(j ,k) = usvsq(j+2,k) *cosfi(j )/cosfi(j+2) usvsq(j+1,k) = usvsq(j+2,k) *cosfi(j+1)/cosfi(j+2) uq(j ,k)=uq(j+2,k)*cosfi(j )/cosfi(j+2) vq(j ,k)=vq(j+2,k)*cosfi(j )/cosfi(j+2) wq(j ,k)=wq(j+2,k)*cosfi(j )/cosfi(j+2) uq(j+1,k)=uq(j+2,k)*cosfi(j+1)/cosfi(j+2) vq(j+1,k)=vq(j+2,k)*cosfi(j+1)/cosfi(j+2) wq(j+1,k)=wq(j+2,k)*cosfi(j+1)/cosfi(j+2) endif if(j.eq.nb) then usvsq(j ,k) = usvsq(j-2,k) *cosfi(j )/cosfi(j-2) usvsq(j-1,k) = usvsq(j-2,k) *cosfi(j-1)/cosfi(j-2) uq(j ,k)=uq(j-2,k)*cosfi(j )/cosfi(j-2) vq(j ,k)=vq(j-2,k)*cosfi(j )/cosfi(j-2) wq(j ,k)=wq(j-2,k)*cosfi(j )/cosfi(j-2) uq(j-1,k)=uq(j-2,k)*cosfi(j-1)/cosfi(j-2) vq(j-1,k)=vq(j-2,k)*cosfi(j-1)/cosfi(j-2) wq(j-1,k)=wq(j-2,k)*cosfi(j-1)/cosfi(j-2) endif enddo enddo 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 uqdz (j,k)=( uq(j,k+1)- uq(j,k-1))/2./dz ! d uquer/dz else if(k.eq.1) then tqdz (j,k)=( tq(j,k+1)- tq(j,k)) / dz uqdz (j,k)=( uq(j,k+1)- uq(j,k)) / dz else if(k.eq.kgit) then tqdz (j,k)=( tq(j,k)- tq(j,k-1)) / dz uqdz (j,k)=( uq(j,k)- uq(j,k-1)) / dz endif if(j.eq.36.and.k.eq.20) &print*,1,tqdz (j,k),tq(j,k+1),tq(j,k-1) if(j.eq.36.and.k.eq.20) &print*,2,uqdz (j,k),uq(j,k+1),uq(j,k-1) c-------------d/dfi(uquer*cosfi) und Fz und d/dfi(tpotq)---------------- if(j.eq.1) then uqdfi (j,k)= (uq(j+1,k)-uq(j,k))/deltf c uqdfi (j,k)= cosfi(j)*(uq(j+1,k)+uq(j,k))/deltf-sinfi(j)*uq(j,k) else if(j.eq.nb) then uqdfi (j,k)= (uq(j,k)-uq(j-1,k))/deltf c uqdfi (j,k)= cosfi(j)*(uq(j-1,k)+uq(j,k))/deltf-sinfi(j)*uq(j,k) else uqdfi (j,k)=(uq(j+1,k)-uq(j-1,k))/2/deltf c uqdfi (j,k)=cosfi(j)*(uq(j+1,k)-uq(j-1,k))/2./deltf- c $ sinfi(j)*uq(j,k) endif enddo enddo do k=1,kgit do j=1,nb Fphi(j,k)=rho0(k)*a*(uqdz(j,k)*vstpsq(j,k)/tqdz(j,k)-usvsq(j,k)) Fz(j,k)=rho0(k)*a*cosfi(j)*((f(j)-1./(a*cosfi(j))* $ uqdfi(j,k))*vstpsq(j,k)/tqdz(j,k)-uswsq(j,k)) if(j.eq.36.and.k.eq.20) &print*,3,Fphi(j,k),rho0(k),uqdz(j,k),vstpsq(j,k),tqdz(j,k), &usvsq(j,k) if(j.eq.36.and.k.eq.20) &print*,3,Fphi(j-1,k),rho0(k),uqdz(j-1,k),vstpsq(j-1,k), &tqdz(j-1,k),usvsq(j-1,k) if(j.eq.36.and.k.eq.20) &print*,4,Fz(j,k),cosfi(j),f(j),1./(a*cosfi(j)), &uqdfi(j,k) ,vstpsq(j,k),tqdz(j,k),uswsq(j,k) c enddo enddo c-------------------------------------------------------------------- do k=1,kgit do j=1,nb if(k.ne.1.and.k.ne.kgit) then dFz (j,k)=(Fz(j,k+1) - Fz(j,k-1))/2./dz if(j.eq.36.and.k.eq.20) &print*,5,dFz(j,k),Fz(j,k+1),Fz(j,k-1),dz 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 dFz (j,k)=(Fz(j,k+1)- Fz(j,k))/dz tpqzz =(tqdz(j,k+1)-tqdz(j,k))/dz vstpsdz =(vstpsq(j,k+1)-vstpsq(j,k))/dz else if(k.eq.kgit) then dFz (j,k)=(Fz(j,k)-Fz(j,k-1))/dz 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), ----------------------------------- c dFphi(j,k)=(Fphi(j+1,k)-Fphi(j,k))/deltf c----correction near pole - interpolation to 0. at the Pole dFphi(j,k)=(Fphi(j+2,k)-Fphi(j,k))/2/deltf/3 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------------------------------------------------------------------- c dFphi(j,k)=(Fphi(j,k)-Fphi(j-1,k))/deltf c----correction near pole - interpolation to 0. at the Pole dFphi(j,k)=dFphi(j-1,k)/3 vstpsfi=(vstpsq(j,k)-vstpsq(j-1,k))/deltf tpqdzfi=(tqdz(j,k)-tqdz(j-1,k))/deltf if(j.eq.36.and.k.eq.20) &print*,6,dFphi(j,k),Fphi(j,k),Fphi(j-1,k),deltf c------------------------------------------------------------------- else c------------------------------------------------------------------- dFphi(j,k)=(Fphi(j+1,k)-Fphi(j-1,k))/2./deltf 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 if(k.eq.20) &write(24,*) j,dFphi(j,k),Fphi(j,k),dFz(j,k),Fz(j,k), &vstpsq(j,k),usvsq(j,k),uswsq(j,k) enddo enddo do k=1,kgit do j=1,nb c-----EP-Flux-Divergence-------------------------------------------- gradF(j,k)=1./(a*cosfi(j)) * dFphi(j,k) + dFz(j,k) if(j.eq.36.and.k.eq.20) &print*,7,rns,dFphi(j,k),dFz(j,k),gradF(j,k)/rho0(k) gradF(j,k)= gradF(j,k) /rho0(k) Fz (j,k)=Fz (j,k) /rho0(k) Fphi (j,k)=Fphi (j,k) /rho0(k) c if(j.eq.10) c &print*,k,gradF(j,k) c gradF(j,k)=gradF(j,k)/rho0(k)*86400./a/cosfi(j) c dFz (j,k)=dFz (j,k)/rho0(k)*86400./a/ cosfi(j) c dFphi(j,k)=dFphi(j,k)/rho0(k)*86400./a/a/cosfi(j)/cosfi(j) c Fz (j,k)=Fz (j,k)/a/rho0(k)/cosfi(j) c Fphi (j,k)=Fphi (j,k)/a/rho0(k) enddo enddo c print*,'---------------------------------------' c=========END CALCULATING RESIDUAL==================================== c uqdz(j,k)*vstpsq(j,k)/tqdz(j,k)-usvsq(j,k) write(16,rec=1) uq,vres,wres,gradF,Fphi,Fz c write(16,rec=prog) dFphi,Fphi,dFz,Fz,usvsq,gradF c============close all============================== close(RD) close(RD2) close(16) ENDDO ! runs STOP END *************************************************************************