C==================================================================== PROGRAM NEWBAL PARAMETER(II=41,JJ=65,KK=3,II1=II-1,II2=II-2, . JJ1=JJ-1,JJ2=JJ-2) C C C THIS PROGRAM DETERMINES THE TIME MEAN VALUES OF THE TERMS C OF THE VORTICITY EQUATION. SECONDLY, THE HORIZONTAL ADVECTION C AND STRETCHING TERMS ARE SPLIT INTO THE COMPONENT DUE TO C THE VELOCITY OF THE MEAN FLOW AND THAT DUE TO THE EDDY C (REYNOLDS) STRESSES. C STRCH: MEAN VALUE OF THE STRETCHING TERM C ESTRC: STRETCHING OF THE EDDY FIELD C HADV: HORIZONTAL ADVECTION OF VORTICITY IN THE MEAN C EHADV: HORIZONTAL ADVECTION OF VORTICITY IN THE EDDY C BETA: HORIZONTAL ADVECTION OF PLANETARY VORTICITY C VADV: VERTICAL ADVECTION C SOLEN: THE SOLENOIDAL TERM C CURLTAU: THE CURL OF WIND FORCING C DISP: THE CURL OF VISCOUS DIFFUSION C DRAG: THE CURL OF BOTTOM DRAG C C THE VARIABLE NRCSPD GIVES THE NUMBER OF DATA RECORDS, NSTART C IS THE BEGINNING OF THE ANALYSIS PERIOD C C BASED ON A CODE WRITTEN BY BOUDRA AND CHASSIGNET, DOCUMENTED C IN THE JOURNAL OF PHYSICAL OCEANOGRAPHY, (1988) 18, 280-303. C C REWRITTEN BY M. SPALL 6/93 C MODIFICATIONS: C 1) NON-CONSERVATIVE HORIZONTAL ADVECTION OPERATOR FOR MOMENTUM C EQUATIONS ("CONSERVATION" TERMS ARE ELIMINATED). THIS IS C CONSISTENT WITH THE FORM OF ADVECTION USED IN THE MODEL OF C SMITH (1992) J. PHYS. OCEAN. 22, 686-696. C 2) PERIODIC BCS IN THE J-DIRECTION C 3) BOTTOM DRAG C 4) HORIZONTAL DIFFUSION OPERATOR MODIFIED C 5) ALLOWS FOR ZERO LAYER THICKNESS C C REAL STRCH(II,JJ),LSTR(II,JJ),ESTRC(II,JJ),MSTRC(II,JJ) REAL HADV(II,JJ),EHADV(II,JJ),MHADV(II,JJ),MADVPV(II,JJ) DIMENSION U(0:II,0:JJ,KK),V(0:II,0:JJ,KK) DIMENSION UM(0:II,0:JJ,KK),VM(0:II,0:JJ,KK) DIMENSION DPM(0:II,0:JJ,KK),DP(0:II,0:JJ,KK) DIMENSION DISP(II,JJ),US(II,JJ),UBOT(II,JJ),VBOT(II,JJ), . VS(II,JJ),PS(II,JJ) DIMENSION VADV(II,JJ),SOLEN(II,JJ) DIMENSION CORIO(II,JJ),BETA(II,JJ) DIMENSION TH(II,JJ,KK*2),P(0:II,0:JJ,KK+1),SDOT(II,JJ,KK+1) DIMENSION THETA(KK),CURLTAU(II,JJ),STRESX(II,JJ),STRESY(II,JJ) DIMENSION STRSX(II,JJ),STRSY(II,JJ),VFLU(II,JJ),VFLV(II,JJ) DIMENSION SM(II,JJ),DRAG(II,JJ),UDRAG(II,JJ),VDRAG(II,JJ) c dimension fbx(ii,jj),fby(ii,jj),dfdy(ii,jj) dimension vbxy(ii,jj),ubxy(ii,jj),dudx(ii,jj),dudy(ii,jj) dimension dvdx(ii,jj),dvdy(ii,jj),ubx(ii,jj),vby(ii,jj) dimension depths(ii,jj),diap(ii,jj),wi(ii,jj,kk+1) dimension depthv(ii,jj),depthu(ii,jj) dimension dpu(ii,jj,kk),dpv(ii,jj,kk) dimension visc(ii,jj),del2(ii,jj) dimension via(ii,jj),vib(ii,jj),s1(ii),s2(ii) dimension proxa(ii),proxb(ii),proxc(ii),proxd(ii) dimension vrt(ii,jj),vrtm(ii,jj),pvm(ii,jj),temp(ii) character*80 filein, fileout C DATA PSTRES/100.E5/ DATA THETA/.956371,.956077,.955953/ DATA SCALE/50.E5/ DATA AHU,AHV/2*8.E6/ DATA S1/-1.,II1*1./,S2/II2*1.,-1.,1./ DATA DELP0/5./ DATA BLB/20./ DATA EPS/1./ DATA G/980.6/ c data delt1/8640./ data iunit/13/ data iout/14/ c cvmgp(a,b,c)=a*(.5+sign(.5,c))+b*(.5-sign(.5,c)) uvdep(a,b)=amin1(a,b) C read(5,'(a80)') filein read(5,'(a80)') fileout read(5,*) nstart,nrcspd read(5,*) factor,vfactor read(5,*) i1,i2,j1,j2 read(5,*) nsmooth read(5,*) cdrag c open(unit=iunit,file=filein,status='old',form='unformatted') open(unit=iout,file=fileout,form='unformatted') c nout=nrcspd-nstart+1 npds=1 c X=1./SCALE read(iunit) ((corio(i,j),i=1,ii),j=1,jj) read(iunit) ((depths(i,j),i=1,ii),j=1,jj) c DO 10 I=1,II DO 10 J=1,JJ STRESX(I,J)=0. STRESY(I,J)=0. 10 continue c DO 11 K=1,KK DO 11 I=1,II DO 11 J=1,JJ TH(I,J,K)=THETA(K) 11 TH(I,J,K+KK)=THETA(K) C do 22 i=2,ii1 do 22 j=1,jj 22 DEPTHU(I,J)=UVDEP(depths(I,J),depths(I-1,J)) do 23 j=1,jj ja=mod(j-2+jj,jj)+1 do 23 i=1,ii1 23 depthv(i,j)=uvdep(depths(i,j),depths(i,ja)) C DO 100 IPD=1,NPDS C DO 110 K=1,KK DO 110 I=1,II1 DO 110 J=1,JJ1 UM(I,J,K)=0. VM(I,J,K)=0. 110 DPM(I,J,K)=0. DO 1 I=1,II1 DO 1 J=1,JJ1 BETA(I,J)=0. STRCH(I,J)=0. ESTRC(I,J)=0. MSTRC(I,J)=0. HADV(I,J)=0. EHADV(I,J)=0. MHADV(I,J)=0. MADVPV(I,J)=0. VADV(I,J)=0. LSTR(I,J)=0. SOLEN(I,J)=0. CURLTAU(I,J)=0. diap(i,j)=0. 1 DISP(I,J)=0. C DO 50 IRC=1,NRCSPD read(iunit) time0,time,nstep,sumpot,sumkin,sumxgr,eke,eape do 200 k=1,kk read(iunit) ((u(i,j,k),i=1,ii),j=1,jj) read(iunit) ((v(i,j,k),i=1,ii),j=1,jj) read(iunit) ((dp(i,j,k),i=1,ii),j=1,jj) 200 continue c if(irc.lt.nstart) go to 49 c c fill no slip boundary conditions at i=0,II call noslip(u) call noslip(v) call noslip(dp) c k=3 do 24 i=2,ii1 do 24 j=1,jj1 dpu(i,j,k)=amin1(.5*(dp(i,j,k)+dp(i-1,j,k)), .amax1(0.,depthu(i,j)-.5*(p(i,j,k)+p(i-1,j,k)))) 24 continue do 25 j=1,jj ja=mod(j-2+jj,jj)+1 do 25 i=1,ii1 25 dpv(i,j,k)=amin1(.5*(dp(i,j,k)+dp(i,ja,k)), .amax1(0.,depthv(i,j)-.5*(p(i,j,k)+p(i,ja,k)))) c NN=0 MM=0 DO 15 K=1,KK KN=K+NN DO 13 I=1,II DO 13 J=1,JJ 13 UM(I,J,K)=UM(I,J,K)+0.5*(dp(i,j,k)+dp(i-1,j,k))*U(I,J,KN)/NOUT DO 14 I=1,II DO 14 J=1,JJ 14 VM(I,J,K)=VM(I,J,K)+0.5*(dp(i,j,k)+dp(i,j-1,k))*V(I,J,KN)/NOUT DO 15 I=1,II DO 15 J=1,JJ 15 DPM(I,J,K)=DPM(I,J,K)+DP(I,J,KN)/NOUT DO 5 K=1,KK KN=K+NN DO 5 I=1,II DO 5 J=1,JJ 5 P(I,J,K+1)=P(I,J,K)+DP(I,J,KN) call noslip(p) C DO 20 K=3,3 C KN=K+NN KM=K+MM C C DO 6 I=1,II DO 6 J=1,JJ if(dp(i,j,km)+dp(i-1,j,km).gt.1) then STRSX(I,J)=STRESX(I,J)*(AMIN1(.5*(P(I,J,K+1)+ .P(I-1,J,K+1)),PSTRES)-AMIN1(.5*(P(I,J,K)+P(I-1,J,K)),PSTRES)) ./(.5*(DP(I,J,KM)+DP(I-1,J,KM))) VFLU(I,J)=((SDOT(I,J,K)+SDOT(I-1,J,K))*(U(I,J,KM)- .U(I,J,KM-1))-(SDOT(I,J,K+1)+SDOT(I-1,J,K+1))*(U(I,J,KM)- .U(I,J,KM+1)))*.5/(DELT1*(DP(I,J,KM)+DP(I-1,J,KM))) else strsx(i,j)=0.0 vflu(i,j)=0.0 endif 6 continue DO 42 J=1,JJ1 DO 42 I=1,II1 PS(I,J)=amax1(.5*(DP(I,J,KM)+DP(I-1,J,KM)),5.) 42 CONTINUE C C do 44 j=1,jj do 44 i=2,ii1 del2(i,j)=-u(i,j,k) 44 visc(i,j)=amax1(dpu(i,j,k),delp0) C do 46 j=1,jj ja=mod(j-2+jj,jj)+1 jb=mod(j ,jj)+1 do 46 i=2,ii1 us(i,j)= .-.5*(((del2(i+1,j )-del2(i,j))*(visc(i,j)+visc(i+1,j )) . +(del2(i-1,j )-del2(i,j))*(visc(i,j)+visc(i-1,j ))) . +((del2(i ,jb)-del2(i,j))*(visc(i,j)+visc(i ,jb)) . +(del2(i ,ja)-del2(i,j))*(visc(i,j)+visc(i ,ja))) ) .*ahu/amax1(dpu(i,j,k),delp0) 46 continue C DO 7 I=1,II DO 7 J=1,JJ if(dp(i,j,km)+dp(i,j-1,km).gt.1) then STRSY(I,J)=STRESY(I,J)*(AMIN1(.5*(P(I,J,K+1)+P(I,J-1,K+1)), .PSTRES)-AMIN1(.5*(P(I,J,K)+P(I,J-1,K)),PSTRES)) ./(.5*(DP(I,J,KM)+DP(I,J-1,KM))) VFLV(I,J)=((SDOT(I,J,K)+SDOT(I,J-1,K))*(V(I,J,KM)- .V(I,J,KM-1))-(SDOT(I,J,K+1)+SDOT(I,J-1,K+1))*(V(I,J,KM) .-V(I,J,KM+1)))*.5/(DELT1*(DP(I,J,KM)+DP(I,J-1,KM))) else strsy(i,j)=0.0 vflv(i,j)=0.0 endif 7 continue c do 54 j=1,jj ja=mod(j-2+jj,jj)+1 jb=mod(j ,jj)+1 do 54 i=1,ii del2(i,j)=-v(i,j,kn) 54 visc(i,j)=amax1(dpv(i,j,k),delp0) C do 93 j=1,jj do 96 i=2,ii1 96 via(i,j)=v(i-1,j,km) do 97 i=1,ii2 97 vib(i,j)=v(i+1,j,km) via(1,j)=-v(1,j,km) 93 vib(ii1,j)=-v(ii1,j,km) C do 56 j=1,jj ja=mod(j-2+jj,jj)+1 jb=mod(j ,jj)+1 do 53 i=1,ii1 proxa(i)=del2(i-1,j) proxb(i)=del2(i+1,j) proxc(i)=visc(i-1,j) 53 proxd(i)=visc(i+1,j) proxa( 1)=del2( 1,j) proxb(ii1)=del2(ii1,j) proxc( 1)=visc( 1,j) proxd(ii1)=visc(ii1,j) c do 56 i=1,ii1 vs(i,j)= .-.5*(((s2(i)*proxb(i) -del2(i,j))*(visc(i,j)+proxd(i) ) . +(s1(i)*proxa(i) -del2(i,j))*(visc(i,j)+proxc(i) )) . +((del2(i,jb )-del2(i,j))*(visc(i,j)+visc(i,jb )) . +(del2(i,ja )-del2(i,j))*(visc(i,j)+visc(i,ja )))) .*ahv/amax1(dpv(i,j,k),delp0) if(dpv(i,j,k).le.delp0) vs(i,j)=0.0 56 continue c C C --- COMPUTE BOTTOM VELOCITIES FOR QUADRATIC BOTTOM DRAG CALCULATION C DO 987 J=1,JJ DO 987 I=1,II1 UBOT(I,J)=0. 987 VBOT(I,J)=0. C DO 988 J=1,JJ JB=MOD(J,JJ)+1 DO 988 I=1,II1 BL=P(I,J,KK+1)-BLB UBOT(I,J)=UBOT(I,J)+ . (U(I,J,K+NN)+U(I+1,J,K+NN))*(AMAX1(P(I,J,K+1),BL) . -AMAX1(P(I,J,K ),BL)) 988 VBOT(I,J)=VBOT(I,J)+ . (V(I,J,K+NN)+V(I,JB,K+NN))*(AMAX1(P(I,J,K+1),BL) . -AMAX1(P(I,J,K ),BL)) C DO 989 J=1,JJ DO 989 I=1,II1 BL=P(I,J,KK+1)-BLB UBOT(I,J)=(UBOT(I,J)/(P(I,J,KK+1)-AMAX1(P(I,J,1),BL)))**2 989 VBOT(I,J)=(VBOT(I,J)/(P(I,J,KK+1)-AMAX1(P(I,J,1),BL)))**2 C DO 990 J=1,JJ DO 990 I=1,II1 990 SM(I,J)=cdrag*(.25*SQRT(UBOT(I,J)+VBOT(I,J)))*G/(BLB*1.e5) C C --- CARRY OUT DRAG CALCULATION FOR U-VELOCITY C DO 159 J=1,JJ DO 159 I=1,II1 DPUU=AMAX1(DP(I,J,K+MM)+DP(I-1,J,K+MM),EPS) UDRAG(I,J)=SM(I,J)*(1./DPUU)*AMAX1(EPS, . AMAX1(P(I,J,K+1)+P(I-1,J,K+1),P(I,J,KK+1)+P(I-1,J,KK+1)-2.*BLB) .-AMAX1(P(I,J,K )+P(I-1,J,K ),P(I,J,KK+1)+P(I-1,J,KK+1)-2.*BLB)) UDRAG(I,J)=CVMGP(UDRAG(I,J),UDRAG(I,J)*.5*DPUU/5.,DPUU*.5 .-5.) UDRAG(I,J)=UDRAG(I,J)*U(I,J,K) 159 continue C C --- CARRY OUT DRAG CALCULATION FOR V-VELOCITY C DO 160 J=1,JJ JA=MOD(J-2+JJ,JJ)+1 DO 160 I=1,II1 DPVV=AMAX1(DP(I,J,K+MM)+DP(I,JA,K+MM),EPS) VDRAG(I,J)=SM(I,J)*(1./DPVV)*AMAX1(EPS, .AMAX1(P(I,J,K+1)+P(I,JA,K+1),P(I,J,KK+1)+P(I,JA,KK+1)-2.*BLB) .-AMAX1(P(I,J,K)+P(I,JA,K),P(I,J,KK+1)+P(I,JA,KK+1)-2.*BLB)) VDRAG(I,J)=CVMGP(VDRAG(I,J),VDRAG(I,J)*.5*DPVV/5.,DPVV*.5 .-5.) VDRAG(I,J)=VDRAG(I,J)*V(I,J,K) 160 continue c c some averages and differences for use in calculating the vorticity terms do 440 j=1,jj1 ja=mod(j-2+jj,jj)+1 jb=mod(j,jj)+1 do 440 i=1,ii1 fbx(i,j)=2*corio(i,j)+corio(i-1,j)+corio(i+1,j) fby(i,j)=2*corio(i,j)+corio(i,ja)+corio(i,jb) dfdy(i,j)=corio(i,jb)-corio(i,ja) vby(i,j)=2*v(i,j,k)+v(i,ja,k)+v(i,jb,k) ubx(i,j)=2*u(i,j,k)+u(i-1,j,k)+u(i+1,j,k) vbxy(i,j)=v(i,j,k)+v(i-1,j,k)+v(i,jb,k)+v(i-1,jb,k) ubxy(i,j)=u(i,j,k)+u(i,ja,k)+u(i+1,j,k)+u(i+1,ja,k) dudx(i,j)=u(i+1,j,k)-u(i-1,j,k) dudy(i,j)=u(i,jb,k)-u(i,ja,k) dvdx(i,j)=v(i+1,j,k)-v(i-1,j,k) dvdy(i,j)=v(i,jb,k)-v(i,ja,k) 440 continue C C C COMPONENTS OF THE VORTICITY EQUATION--CALCULATED AT THE VORTICITY C POINTS. SINCE MANY OF THE AVERAGING OPERATORS REQUIRE C USE OF TWO FULL GRID INCREMENTS, WE IGNORE THE FIRST AND C LAST ROWS AND COLUMNS OF THE DOMAIN. IN FINDING THE VERTICAL C MEAN VALUE, THE VORTICITY OF EACH LAYER IS WEIGHTED FOR C ITS PRESSURE THICKNESS. C C DO 8 I=2,II1 DO 8 J=2,JJ1 C C N-L STRETCHING TERM C STRCH(I,J)=STRCH(I,J)-.0625*X*X*( . (vby (i,j)-vby (i-1,j))*(dvdy(i,j)+dvdy(i-1,j)) . -(vbxy(i,j)-vbxy(i,j-1))*(dudy(i,j)+dudy(i,j-1)) . +(ubxy(i,j)-ubxy(i-1,j))*(dvdx(i,j)+dvdx(i-1,j)) . -(ubx (i,j)-ubx (i,j-1))*(dudx(i,j)+dudx(i,j-1)) ) . /FLOAT(NOUT) C C --- LINEAR STRETCHING TERM F (DIV V) C LSTR(I,J)=LSTR(I,J)-.0625*X*( . (vbxy(i,j)-vbxy(i,j-1))*fby (i,j) . +(ubxy(i,j)-ubxy(i-1,j))*fbx (i,j) ) ./FLOAT(NOUT) c c --- further divide linear stretching term into adiabatic and diabatic c contributions c c if(fluxtype.eq.1) deltaz=wi(3) c if(fluxtype.eq.2) deltaz=wi(3)*2/(dp(i,j,kn-1)+dp(i,j,kn)) c if(fluxtype.eq.3) deltaz=wi(3)/amax1(dp(i,j,kn),100*1.e+5) c if(dp(i,j,kn)-delt1*deltaz.lt.0) deltaz=dp(i,j,kn)/delt1 c diap(i,j)=diap(i,j)-plvrt(i,j)*.25*x*(wi(k)-wi(k+1)) C C C --- RELATIVE VORTICITY ADVECTION C HADV(I,J)=HADV(I,J)-.0625*X*X*( . (ubxy(i,j)+ubxy(i-1,j))*(dvdx(i,j)-dvdx(i-1,j)) . -(ubx (i,j)+ubx (i,j-1))*(dudx(i,j)-dudx(i,j-1)) . +(vby (i,j)+vby (i-1,j))*(dvdy(i,j)-dvdy(i-1,j)) . -(vbxy(i,j)+vbxy(i,j-1))*(dudy(i,j)-dudy(i,j-1)) ) ./FLOAT(NOUT) C C --- PLANETARY VORTICITY ADVECTION (BETA) C BETA(I,J)=BETA(I,J)-.0625*X*( . (vbxy(i,j)+vbxy(i,j-1))*dfdy(i,j) ) ./FLOAT(NOUT) C C C CURL OF THE VERTICAL ADVECTION C VADV(I,J)=VADV(I,J)-(VFLV(I,J)-VFLV(I-1,J)-VFLU(I,J)+VFLU(I,J-1) .)*X/NOUT C C C THE SOLENOIDAL TERM (GRAD-P CROSS GRAD-ALPHA) C SOLEN(I,J)=SOLEN(I,J)+(-.25*(TH(I,J,KN)+TH(I-1,J,KN)) .*(P(I,J,K)+P(I,J,K+1)-P(I-1,J,K)-P(I-1,J,K+1))+ ..25*(TH(I,J-1,KN)+TH(I-1,J-1,KN))*(P(I,J-1,K)+P(I,J-1,K+1) .-P(I-1,J-1,K) .-P(I-1,J-1,K+1))+.25*(TH(I,J,KN)+TH(I,J-1,KN)) .*(P(I,J,K)+P(I,J,K+1)-P(I,J-1,K)-P(I,J-1,K+1))- ..25*(TH(I-1,J,KN)+TH(I-1,J-1,KN))*(P(I-1,J,K)+P(I-1,J,K+1) .-P(I-1,J-1,KN)-P(I-1,J-1,K+1)))*X*X/NOUT C C CURL OF THE WIND STRESS FORCING C CURLTAU(I,J)=CURLTAU(I,J)+(STRSY(I,J)-STRSY(I-1,J)-STRSX(I,J) .+STRSX(I,J-1))*X/NOUT C C CURL OF VISCOUS DIFFUSION C DISP(I,J)=DISP(I,J)+((VS(I,J)-VS(I-1,J)) . -(US(I,J)-US(I,J-1))) .*(X**3)/NOUT c c curl of bottom drag c DRAG(I,J)=DRAG(I,J)- . (VDRAG(I,J)-VDRAG(I-1,J)-(UDRAG(I,J)-UDRAG(I,J-1)))*X/NOUT c c relative vorticity vrt(i,j)=1.e+6*x*(dvdx(i,j)-dudy(i,j)) 8 continue write(1) ((vrt(i,j),i=1,ii),j=1,jj) C 20 CONTINUE C 49 continue 50 CONTINUE C DO 52 K=3,3 DO 51 I=1,II DO 51 J=1,JJ if(dpm(i,j,k)+dpm(i-1,j,k).ne.0.) then UM(I,J,K)=UM(I,J,K)/(.5*(DPM(I,J,K)+DPM(I-1,J,K))) else UM(I,J,K)=0.0 endif 51 continue C DO 52 I=1,II DO 52 J=1,JJ if(dpm(i,j,k)+dpm(i,j-1,k).ne.0.) then VM(I,J,K)=VM(I,J,K)/(.5*(DPM(I,J,K)+DPM(I,J-1,K))) else VM(I,J,K)=0.0 endif 52 continue C C CALCULATION OF VORTICITY OF THE MEAN FLOW FIELD C C C DO 80 K=3,3 C C KN=K+NN KM=K+MM C c calculate averages and differences for mean terms c do 74 j=1,jj do 74 i=1,ii vby(i,j)=0. ubx(i,j)=0. vbxy(i,j)=0. ubxy(i,j)=0. dudx(i,j)=0. dudy(i,j)=0. dvdx(i,j)=0. dvdy(i,j)=0. 74 continue do 75 j=1,jj1 ja=mod(j-2+jj,jj)+1 jb=mod(j,jj)+1 do 75 i=1,ii1 vby(i,j)=2*vm(i,j,k)+vm(i,ja,k)+vm(i,jb,k) ubx(i,j)=2*um(i,j,k)+um(i-1,j,k)+um(i+1,j,k) vbxy(i,j)=vm(i,j,k)+vm(i-1,j,k)+vm(i,jb,k)+vm(i-1,jb,k) ubxy(i,j)=um(i,j,k)+um(i,ja,k)+um(i+1,j,k)+um(i+1,ja,k) dudx(i,j)=um(i+1,j,k)-um(i-1,j,k) dudy(i,j)=um(i,jb,k)-um(i,ja,k) dvdx(i,j)=vm(i+1,j,k)-vm(i-1,j,k) dvdy(i,j)=vm(i,jb,k)-vm(i,ja,k) vrtm(i,j)=1.e+6*( (vm(i,j,k)-vm(i-1,j,k))*x . -(um(i,j,k)-um(i,j-1,k))*x ) pvm(i,j)=((vm(i,j,k)-vm(i-1,j,k)-(um(i,j,k)-um(i,j-1,k)))*x . + corio(i,j) )*4./ . (dpm(i,j,k)+dpm(i-1,j,k)+dpm(i,j-1,k)+dpm(i-1,j-1,k)) 75 continue C c k=3 jbar1=15 jbar2=jj-5 njbar=jbar2-jbar1+1 DO 80 I=2,II1 DO 80 J=2,JJ1 C C STRETCHING OF THE MEAN FLOW VORTICITY BY THE MEAN FLOW C (both linear and nonlinear terms) C MSTRC(I,J)=MSTRC(I,J) . -.0625*X*X*( . (vby (i,j)-vby (i-1,j))*(dvdy(i,j)+dvdy(i-1,j)) . -(vbxy(i,j)-vbxy(i,j-1))*(dudy(i,j)+dudy(i,j-1)) . +(ubxy(i,j)-ubxy(i-1,j))*(dvdx(i,j)+dvdx(i-1,j)) . -(ubx (i,j)-ubx (i,j-1))*(dudx(i,j)+dudx(i,j-1)) ) . -.0625*X*( . (vbxy(i,j)-vbxy(i,j-1))*fby (i,j) . +(ubxy(i,j)-ubxy(i-1,j))*fbx (i,j) ) c C --- RELATIVE VORTICITY ADVECTION C MHADV(I,J)=MHADV(I,J)-.0625*X*X*( . (ubxy(i,j)+ubxy(i-1,j))*(dvdx(i,j)-dvdx(i-1,j)) . -(ubx (i,j)+ubx (i,j-1))*(dudx(i,j)-dudx(i,j-1)) . +(vby (i,j)+vby (i-1,j))*(dvdy(i,j)-dvdy(i-1,j)) . -(vbxy(i,j)+vbxy(i,j-1))*(dudy(i,j)-dudy(i,j-1)) ) c c --- mean advection of potential vorticity (v dot grad pv) madvpv(i,j)=madvpv(i,j)+0.25*x*( . (um(i,j,k)+um(i-1,j,k)+um(i,j-1,k)+um(i-1,j-1,k))* . (pvm(i,j)-pvm(i-1,j)) . +(vm(i,j,k)+vm(i-1,j,k)+vm(i,j-1,k)+vm(i-1,j-1,k))* . (pvm(i,j)-pvm(i,j-1))) c if(j.ge.jbar1.and.j.le.jbar2) c . temp(i)=temp(i)+madvpv(i,j)/njbar 80 continue c c --- rescale terms call rescale(strch,factor) call rescale(lstr,factor) call rescale(mstrc,factor) call rescale(hadv,factor) call rescale(mhadv,factor) call rescale(vadv,factor) call rescale(solen,factor) call rescale(beta,factor) call rescale(disp,factor) call rescale(drag,factor) c call rescale(madvpv,vfactor) c C --- smooth variables do 90 ns=1,nsmooth call smooth(strch) call smooth(lstr) call smooth(mstrc) call smooth(hadv) call smooth(mhadv) call smooth(vadv) call smooth(solen) call smooth(beta) call smooth(disp) call smooth(drag) 90 continue C DO 401 I=1,II DO 401 J=1,JJ 401 SM(I,J)=STRCH(I,J)+HADV(I,J)+BETA(I,J)+DISP(I,J)+CURLTAU(I,J) .+SOLEN(I,J)+LSTR(I,J)+DRAG(I,J) C DO 86 I=1,II DO 86 J=1,JJ c if(i.eq.8) write(6,*) j,vrtm(i,j),sm(i,j) ESTRC(I,J)=STRCH(I,J)+LSTR(I,J)-MSTRC(I,J) 86 EHADV(I,J)=HADV(I,J)-MHADV(I,J) c c output terms to unit iout write(iout) ((um(i,j,3),i=i1,i2),j=j1,j2) write(iout) ((vm(i,j,3),i=i1,i2),j=j1,j2) write(iout) ((dpm(i,j,3),i=i1,i2),j=j1,j2) write(iout) ((vrtm(i,j),i=i1,i2),j=j1,j2) write(iout) ((pvm(i,j),i=i1,i2),j=j1,j2) write(iout) ((madvpv(i,j),i=i1,i2),j=j1,j2) c write(iout) ((strch(i,j),i=i1,i2),j=j1,j2) write(iout) ((lstr(i,j),i=i1,i2),j=j1,j2) write(iout) ((mstrc(i,j),i=i1,i2),j=j1,j2) write(iout) ((estrc(i,j),i=i1,i2),j=j1,j2) c write(iout) ((hadv(i,j),i=i1,i2),j=j1,j2) write(iout) ((mhadv(i,j),i=i1,i2),j=j1,j2) write(iout) ((ehadv(i,j),i=i1,i2),j=j1,j2) write(iout) ((solen(i,j),i=i1,i2),j=j1,j2) c write(iout) ((beta(i,j),i=i1,i2),j=j1,j2) write(iout) ((disp(i,j),i=i1,i2),j=j1,j2) write(iout) ((drag(i,j),i=i1,i2),j=j1,j2) write(iout) ((sm(i,j),i=i1,i2),j=j1,j2) write(6,*) 'output array dimensioned ',i2-i1+1,j2-j1+1 C c 100 CONTINUE STOP END c c subroutine smooth(ex1) parameter(ii=41,jj=65,kk=3,ii1=ii-1,ii2=ii-2, . jj1=jj-1,jj2=jj-2) dimension sm(ii,jj),ex1(ii,jj) C C --- SMOOTH EX1 C DO 10 I=1,II1 DO 10 J=1,JJ1 10 SM(I,J)=EX1(I,J) C DO 20 I=3,ii2 DO 20 J=3,jj2 20 SM(I,J)=EX1(I,J)+.125*(EX1(I+1,J)+EX1(I-1,J)+EX1(I,J+1) .+EX1(I,J-1)-4.*EX1(I,J))+.0625*(EX1(I+1,J+1)+EX1(I+1,J-1) .+EX1(I-1,J+1)+EX1(I-1,J-1)-4.*EX1(I,J)) DO 30 I=3,ii2 DO 30 J=3,jj2 30 EX1(I,J)=SM(I,J)-.3952*(SM(I+1,J)+SM(I-1,J)+SM(I,J+1) .+SM(I,J-1)-4.*SM(I,J))+.0676*(SM(I+1,J+1)+SM(I-1,J+1)+ .SM(I+1,J-1)+SM(I-1,J-1)-4.*SM(I,J)) C return end c subroutine noslip(fld) parameter(ii=41,jj=65,kk=3,ii1=ii-1,ii2=ii-2, . jj1=jj-1,jj2=jj-2) dimension fld(0:ii,0:jj) c do 10 j=1,jj fld(0,j)=-fld(1,j) 10 fld(ii,j)=-fld(ii-1,j) c return end c subroutine rescale(fld,factor) parameter(ii=41,jj=65,kk=3,ii1=ii-1,ii2=ii-2, . jj1=jj-1,jj2=jj-2) dimension fld(ii,jj) c do 10 j=1,jj do 10 i=1,ii1 fld(i,j)=fld(i,j)*factor 10 continue c return end