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.,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/'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) f(j) =2.*omega*sinfi(j) enddo c----------------- NEUDEFINITION VON defi ==> deltf deltf=defi*pi/180. c================== ALL RUNS ================================= DO rns = 1,3 ! runs RD = 21 RD2 = 22 c================== DATA OUTPUT ================================= open(16,FILE="../../Muam_Mar_eQBO/" //jfilein(rns)// "/resEP.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 ================================= 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. 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 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 duq(i,j,k)=uvt(i,j,k,1)-uq(j,k) dvq(i,j,k)=uvt(i,j,k,2)-vq(j,k) dwq(i,j,k)=wvel(i,j,k )-wq(j,k) dtq(i,j,k)=tq3(i,j,k)-tq(j,k) 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 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 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 c-------------d/dfi(uquer*cosfi) und Fz und d/dfi(tpotq)---------------- if(j.eq.1) then c uqdfi (j,k)= (uq(j+1,k)-uq(j,k))/deltf uqdfi (j,k)= cosfi(j)*(uq(j+1,k)+uq(j,k))/deltf-sinfi(j)*uq(j,k) else if(j.eq.nb) then c uqdfi (j,k)= (uq(j,k)-uq(j-1,k))/deltf uqdfi (j,k)= cosfi(j)*(uq(j-1,k)+uq(j,k))/deltf-sinfi(j)*uq(j,k) else c uqdfi (j,k)=(uq(j+1,k)-uq(j-1,k))/2/deltf uqdfi (j,k)=cosfi(j)*(uq(j+1,k)-uq(j-1,k))/2./deltf- $ sinfi(j)*uq(j,k) endif 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)) 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 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 dFphi(j,k)=-2.*cosfi(j)*sinfi(j)* Fphi(j ,k)+ $ cosfi(j)*cosfi(j)*(Fphi(j+1,k)-Fphi(j,k))/deltf 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 dFphi(j,k)= -2.*cosfi(j)*sinfi(j)* Fphi(j ,k)+ $ cosfi(j)*cosfi(j)*(Fphi(j-1,k)-Fphi(j,k))/deltf vstpsfi=(vstpsq(j,k)-vstpsq(j-1,k))/deltf tpqdzfi=(tqdz(j,k)-tqdz(j-1,k))/deltf c------------------------------------------------------------------- else c------------------------------------------------------------------- c dFphi(j,k)=(Fphi(j+1,k)-Fphi(j-1,k))/2./deltf dFphi(j,k)=-2.*cosfi(j)*sinfi(j)* Fphi(j ,k)+ $ cosfi(j)*cosfi(j)*(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 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) c if(j.eq.10) c &print*,k,gradF(j,k),gradF(j,k) * exp(z(k)/10000) 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=prog) uq,vres,wres,gradF,Fphi,Fz c write(16,rec=prog) dFphi,Fphi,dFz,Fz,usvsq,gradF END DO ! time loop c============close all============================== close(RD) close(RD2) close(16) ENDDO ! runs STOP END *************************************************************************