subroutine relax c c======================================================================= c === #ifdef sor5pt c RELAX takes as input the vorticity driving function computed === c in CLINIC (ZTD) and, using sequential over-relaxation, === c solves the Laplacian equation for the external mode of === c velocity in terms of a mass transport stream function === c (P). === #else c RELAX solves for the time derivative of the mass transport === c stream function using a pre-conditioned conjugate === c gradient algorithm with 5-point Laplacian. === c === c NOTE: abrupt topographic changes may lead to non === c convergence. In practice this may only be a === c problem for coarse resolution model with non smooth === c topography. === c === c This algorithm is described in GFDL technical report === c by J. Derber, 1987. (See also: Gill, P. E., W. Murray === c and M. H. Wright, 1981, Practical Optimizations, === c Academic Press). === #endif c === c======================================================================= c c----------------------------------------------------------------------- c Define global data c----------------------------------------------------------------------- c #include #include #include #include #include #include #include #include #include #include #include #include #ifdef bndy_rlx # include #endif c c----------------------------------------------------------------------- c Define local and equivalence data. c----------------------------------------------------------------------- c integer i,iend,istr,j,l,luptd,luptdb #ifdef shapiro * ,nn #endif #ifdef coast integer ip,ismask(imt,jmt),itest # ifdef islands integer isle,jstr,jend FLOAT * cof(imt,jmt),cofis(misle),fxd,resis # endif #endif FLOAT * fx,fxa,fxb,fxc,resmax #ifndef sor5pt FLOAT * asor,apr2,apr4,beta,oapr #endif #ifdef vel_conv FLOAT * dumm,vbar_test #endif FLOAT * ptdb(imt,jmt) equivalence (res,ptdb) c c----------------------------------------------------------------------- c Begin executable code. c----------------------------------------------------------------------- c c======================================================================= c Begin introductory section to prepare for the relaxation =========== c======================================================================= c c----------------------------------------------------------------------- c Initiate reading of relaxation solution of 2 timesteps back c for the purpose of computing an initial guess for the present pass c (input unit alternates on timesteps between 5 & 6) c----------------------------------------------------------------------- c luptdb=6-mod(itt,2) luptd =5+mod(itt,2) call ofind(kflds,nwds,(luptdb-1)*nwds+1) c c----------------------------------------------------------------------- c Initialize the work area c----------------------------------------------------------------------- c do 10 j=1,jmt do 10 i=1,imt cfn(i,j)=c0 cfs(i,j)=c0 cfe(i,j)=c0 cfw(i,j)=c0 #ifdef coast ismask(i,j)=0 #endif #ifdef islands cof(i,j)=c0 #endif 10 continue c c----------------------------------------------------------------------- c Complete reading of relaxation solution of 2 timesteps back and c initiate reading of relaxation solution of previous timestep c (input unit alternates on timesteps between 6 & 5) c----------------------------------------------------------------------- c call oget(kflds,nwds,(luptdb-1)*nwds+1,ptdb) call ofind(kflds,nwds,(luptd-1)*nwds+1) #ifdef coast c c----------------------------------------------------------------------- c Form mask array ISMASK=0 interior ocean c =1 perimeter c =2 land c----------------------------------------------------------------------- c do 20 j=2,jmt do 20 i=2,imt itest=landv(i-1,j)*landv(i,j-1)*landv(i-1,j-1)*landv(i,j) if(itest.eq.0) ismask(i,j)=1 itest=landv(i-1,j)+landv(i,j-1)+landv(i-1,j-1)+landv(i,j) if(itest.eq.0) ismask(i,j)=2 20 continue #endif c c----------------------------------------------------------------------- c Filter/smooth vorticity field. c----------------------------------------------------------------------- #ifdef shapiro c c If requested, apply Shapiro filter: c c ICNTZ = counter (for frequency of application) c NFRQZ = frequency with which filter is applied c NTIMZ = number of times filter is applied per time step c NORDZ = order of the filter c if(mixztd.eq.1) then icntz=icntz+1 if(icntz.eq.nfrqz)then icntz=0 if(nordz.gt.0)then do 30 nn=1,ntimz # ifndef cyclic call shap_lev(ztd,imtm1,jmtm1,zgrid,nordz) # else call shap_lev(ztd,imt,jmtm1,zgrid,nordz) # endif 30 continue endif endif endif #endif c c Or, apply Laplacian filter. c if(mixztd.eq.2) then call lap_filt(ztd,tgrid) endif c c Save rate change of barotropic vorticity in array ZTDB for output c purposes. c do 40 j=1,jmt do 40 i=1,imt ztdb(i,j)=ztd(i,j) 40 continue c c----------------------------------------------------------------------- c Generate arrays of coefficients for relaxation. c----------------------------------------------------------------------- c c 1st, compute coefficients of the Laplacian star. c (hold non-interior points to zero using start and end indices) c do 80 j=2,jmtm1 do 80 l=1,lseg istr=isz(j,l) if(istr.eq.0) goto 80 iend=iez(j,l) fxa=c2*cstr(j)*cstr(j) fxb=c2*cs(j )*cstr(j)*dytr(j)*dyur(j ) fxc=c2*cs(j-1)*cstr(j)*dytr(j)*dyur(j-1) do 50 i=istr,iend cfn(i,j)=fxb/(hdv(i-1,j)+hdv(i,j)) cfs(i,j)=fxc/(hdv(i-1,j-1)+hdv(i,j-1)) cfe(i,j)=fxa*dxur(i)*dxtr(i)/(hdv(i,j)+hdv(i,j-1)) cfw(i,j)=fxa*dxur(i-1)*dxtr(i)/(hdv(i-1,j)+hdv(i-1,j-1)) #if defined bndy_rlx & defined imp_bnd_rlx fx =p25*c2dtuv cfn(i,j)=cfn(i,j)*(c1+fx*(vfacbrlx(i-1,j)+vfacbrlx(i,j))) cfs(i,j)=cfs(i,j)*(c1+fx*(vfacbrlx(i-1,j-1)+vfacbrlx(i,j-1))) cfe(i,j)=cfe(i,j)*(c1+fx*(vfacbrlx(i,j)+vfacbrlx(i,j-1))) cfw(i,j)=cfw(i,j)*(c1+fx*(vfacbrlx(i-1,j)+vfacbrlx(i-1,j-1))) #endif cpf(i)=c1/(cfn(i,j)+cfs(i,j)+cfe(i,j)+cfw(i,j)) 50 continue c c 2nd, augment coefficients for implicit treatment of Coriolis term. c The update of coefficient CPF was missing in the original code c leading to an error while setting the diagonal of the elliptic c operator to 1. For clarity it is put here even though the line c in loop 60 is not required (CJL). c if(acor.ne.c0) then fx=-c2dtsf*acor*cstr(j)*dytr(j)*omega do 60 i=istr,iend cfn(i,j)=cfn(i,j)+(hv(i,j )*sine(i,j ) * -hv(i-1,j )*sine(i-1,j ))*fx*dxtr(i) cfs(i,j)=cfs(i,j)-(hv(i,j-1)*sine(i,j-1) * -hv(i-1,j-1)*sine(i-1,j-1))*fx*dxtr(i) cfe(i,j)=cfe(i,j)-(hv(i ,j)*sine(i,j)-hv(i ,j-1)* * sine(i,j-1))*fx*dxtr(i) cfw(i,j)=cfw(i,j)+(hv(i-1,j)*sine(i-1,j)-hv(i-1,j-1)* * sine(i-1,j-1))*fx*dxtr(i) cpf(i)=c1/(cfn(i,j)+cfs(i,j)+cfe(i,j)+cfw(i,j)) 60 continue endif c c 3rd, normalize coefficients and forcing term for efficiency. c do 70 i=istr,iend cfn(i,j)=cfn(i,j)*cpf(i) cfs(i,j)=cfs(i,j)*cpf(i) cfe(i,j)=cfe(i,j)*cpf(i) cfw(i,j)=cfw(i,j)*cpf(i) ztd(i,j)=ztd(i,j)*cpf(i) 70 continue 80 continue #ifdef islands c c 4th, compute coefficients on island perimeter points. c do 100 isle=1,nisle cofis(isle)=c0 istr=isis(isle) iend=ieis(isle) jstr=jsis(isle) jend=jeis(isle) do 90 j=jstr,jend fxa=c2*cstr(j)*cstr(j) fxb=c2*cs(j )*dyur(j )*dytr(j)*cstr(j) fxc=c2*cs(j-1)*dyur(j-1)*dytr(j)*cstr(j) fxd=-c2dtsf*acor*cstr(j)*dytr(j)*omega do 90 i=istr,iend if(ismask(i,j).ne.1) goto 90 if((landv(i-1,j ).ne.0).or.(landv(i ,j ).ne.0)) * cfn(i,j)=fxb/(hdv(i-1,j)+hdv(i,j)) * +fxd*dxtr(i)*(hv(i,j )*sine(i,j ) * -hv(i-1,j )*sine(i-1,j )) if((landv(i-1,j-1).ne.0).or.(landv(i ,j-1).ne.0)) * cfs(i,j)=fxc/(hdv(i-1,j-1)+hdv(i,j-1)) * -fxd*dxtr(i)*(hv(i,j-1)*sine(i,j-1) * -hv(i-1,j-1)*sine(i-1,j-1)) if((landv(i ,j ).ne.0).or.(landv(i ,j-1).ne.0)) * cfe(i,j)=fxa*dxtr(i)*dxur(i)/(hdv(i,j)+hdv(i,j-1)) * -fxd*dxtr(i)*(hv(i,j)*sine(i,j)-hv(i,j-1)* * sine(i,j-1)) if((landv(i-1,j ).ne.0).or.(landv(i-1,j-1).ne.0)) * cfw(i,j)=fxa*dxtr(i)*dxur(i-1)/(hdv(i-1,j)+hdv(i-1,j-1)) * +fxd*dxtr(i)*(hv(i-1,j)*sine(i-1,j)-hv(i-1,j-1)* * sine(i-1,j-1)) cof(i,j)=c1/(cfn(i,j)+cfs(i,j)+cfe(i,j)+cfw(i,j)) cfn(i,j)=cfn(i,j)*cof(i,j) cfs(i,j)=cfs(i,j)*cof(i,j) cfe(i,j)=cfe(i,j)*cof(i,j) cfw(i,j)=cfw(i,j)*cof(i,j) ztd(i,j)=ztd(i,j)*cof(i,j) cof(i,j)=cst(j)*dxt(i)*dyt(j)/cof(i,j) cofis(isle)=cofis(isle)+cof(i,j) 90 continue cofis(isle)=c1/cofis(isle) 100 continue #endif c c----------------------------------------------------------------------- c Compute a first guess for the relaxation by extrapolating the two c previous solutions forward in time. c----------------------------------------------------------------------- c c 1st, complete readin of relaxation solution of previous timestep. c call oget(kflds,nwds,(luptd-1)*nwds+1,ptd) c c 2nd, perform time extrapolation (accounting for mixing timestep). c fxa=c1 fxb=c2 if(mixts.or.mxpas2) fxa=p5 do 110 j=1,jmt do 110 i=1,imt ptd(i,j)=fxa*(fxb*ptd(i,j)-ptdb(i,j)) 110 continue c c Update the time change of stream function at the southern and c northern boundaries. c do 120 i=2,imtm1 ptd(i,1)=po(i,1,south)-pb(i,1) ptd(i,jmt)=po(i,1,north)-pb(i,jmt) 120 continue #ifdef cyclic c c Set Cyclic boundary conditions. c do 130 j=1,jmt ptd(1 ,j)=ptd(imtm1,j) ptd(imt,j)=ptd(2 ,j) 130 continue #else c c Update the time change of stream function at the western and c eastern boundaries. c do 140 j=1,jmt ptd(1,j)=po(j,1,west)-pb(1,j) ptd(imt,j)=po(j,1,east)-pb(imt,j) 140 continue #endif #ifdef coast c c Pass to coastal segments with end points at domain boundaries. c do 150 l=1,ncseg fx=ptd(icoast(1,l),jcoast(1,l)) do 150 ip=2,lencoast(l) ptd(icoast(ip,l),jcoast(ip,l))=fx 150 continue #endif c c Compute criterion for convergence of relaxation and set residuals to c zero. c prms=c0 do 170 j=1,jmt do 170 i=1,imt prms=prms+p(i,j)*p(i,j) 170 continue prms=sqrt(prms/FLoaT(imt*jmt)) prms=max(prms,c1) do 180 j=1,jmt do 180 i=1,imt res(i,j)=c0 180 continue c #ifdef vel_conv c----------------------------------------------------------------------- c Initialize velocity convergence routine. c----------------------------------------------------------------------- c dumm = vbar_test (ptd,.true.) c #endif c======================================================================= c End introductory section =========================================== c======================================================================= c #ifdef sor5pt c c======================================================================= c Begin section to do the relaxation ================================= c======================================================================= c c Multiply the coefficients by the over-relaxation factor. c do 190 j=1,jmt do 190 i=1,imt cfn(i,j)=cfn(i,j)*sor cfs(i,j)=cfs(i,j)*sor cfe(i,j)=cfe(i,j)*sor cfw(i,j)=cfw(i,j)*sor 190 continue c c----------------------------------------------------------------------- c Begin iteration loop using over-relaxation technique. c----------------------------------------------------------------------- c mscan=0 200 mscan=mscan+1 c c----------------------------------------------------------------------- c Compute entire field of residuals as in simultaneous relaxation. #ifndef islands c (..note.. residuals will remain zero over non-interior points) #endif c----------------------------------------------------------------------- c do 205 j=1,jmt do 205 i=1,imt res(i,j)=c0 205 continue c do 210 j=2,jmtm1 do 210 i=2,imtm1 res(i,j)=cfn(i,j)*ptd(i,j+1) * +cfs(i,j)*ptd(i,j-1) * +cfe(i,j)*ptd(i+1,j) * +cfw(i,j)*ptd(i-1,j) * -sor*(ptd(i,j)+ztd(i,j)) 210 continue #ifdef islands c c----------------------------------------------------------------------- c Reset residuals to zero over land. c----------------------------------------------------------------------- c do 260 j=2,jmtm1 do 220 i=1,imt cpf(i)=c0 220 continue do 230 l=1,lseg istr=isz(j,l) if(istr.eq.0) goto 240 iend=iez(j,l) do 230 i=istr,iend cpf(i)=res(i,j) 230 continue 240 continue do 250 i=1,imt res(i,j)=cpf(i) 250 continue 260 continue #endif #ifdef cyclic do 270 j=2,jmtm1 res(1 ,j)=res(imtm1,j) res(imt,j)=res(2 ,j) 270 continue #endif c c----------------------------------------------------------------------- c Perform correction on southern point to yield sequential relaxation c----------------------------------------------------------------------- c do 330 j=2,jmtm1 #ifndef islands do 280 i=2,imtm1 res(i,j)=res(i,j)+cfs(i,j)*res(i,j-1) 280 continue #else do 290 l=1,lseg istr=isz(j,l) if(istr.eq.0) goto 300 iend=iez(j,l) do 290 i=istr,iend res(i,j)=res(i,j)+cfs(i,j)*res(i,j-1) 290 continue 300 continue #endif c c----------------------------------------------------------------------- c Perform correction on western point to yield sequential relaxation c----------------------------------------------------------------------- c do 320 l=1,lseg istr=isz(j,l) if(istr.eq.0) goto 330 iend=iez(j,l) do 310 i=istr+1,iend res(i,j)=res(i,j)+cfw(i,j)*res(i-1,j) 310 continue 320 continue 330 continue c c----------------------------------------------------------------------- c Make a correction to ptd based on the residuals c----------------------------------------------------------------------- c do 340 j=2,jmtm1 do 340 i=2,imtm1 ptd(i,j)=ptd(i,j)+res(i,j) 340 continue c c----------------------------------------------------------------------- c Find the maximum absolute residual to determine convergence c----------------------------------------------------------------------- c resmax=c0 do 350 j=2,jmtm1 do 350 i=2,imtm1 resmax=max(abs(res(i,j)),resmax) 350 continue # ifdef islands c c----------------------------------------------------------------------- c Perform whole relaxation for each island c----------------------------------------------------------------------- c do 380 isle=1,nisle istr=isis(isle) iend=ieis(isle) jstr=jsis(isle) jend=jeis(isle) resis=c0 do 360 j=jstr,jend do 360 i=istr,iend if(ismask(i,j).eq.1) then resis=resis+(cfn(i,j)*ptd(i ,j+1) * +cfs(i,j)*ptd(i ,j-1) * +cfe(i,j)*ptd(i+1,j ) * +cfw(i,j)*ptd(i-1,j ) * -sor*(ptd(i,j)+ztd(i,j)))*cof(i,j) endif 360 continue c c Normalize the island residual and update the maximum absolute c residual of the relaxation, if necessary. c resis=resis*cofis(isle) resmax=max(abs(resis),resmax) c c Make a correction to PTD over the island and its perimeter points. c do 370 j=jstr,jend do 370 i=istr,iend if(ismask(i,j).ge.1) then ptd(i,j)=ptd(i,j)+resis endif 370 continue 380 continue # ifdef cyclic c c Set Cyclic boundary conditions. c do 390 j=1,jmt ptd(1 ,j)=ptd(imtm1,j) ptd(imt,j)=ptd(2 ,j) 390 continue # endif # endif c c----------------------------------------------------------------------- c Test maximum residual for convergence of the relaxation. c If not converged, proceed with another scan. c (..note.. if the number of scans reaches MXSCAN, leave the loop) c----------------------------------------------------------------------- c #ifndef vel_conv if((resmax.ge.crit*prms).and.(mscan.lt.mxscan)) go to 200 # else if( ((resmax.ge.crit*prms).or.(vbar_test(ptd,.false.).ge.crit)) * .and.(mscan.lt.mxscan)) go to 200 #endif c c======================================================================= c End of the relaxation ============================================== c======================================================================= c #else c c======================================================================= c Begin section to do 5-point conjugate gradient ====================== c======================================================================= c c----------------------------------------------------------------------- c Find initial residuals. c----------------------------------------------------------------------- c do 400 j=2,jmtm1 #ifdef coast do 400 l=1,lseg istr=isz(j,l) if(istr.eq.0) goto 400 iend=iez(j,l) do 395 i=istr,iend #else do 400 i=2,imtm1 #endif res(i,j)=cfn(i,j)*ptd(i,j+1)+ * cfs(i,j)*ptd(i,j-1)+ * cfe(i,j)*ptd(i+1,j)+ * cfw(i,j)*ptd(i-1,j)- * ptd(i,j)-ztd(i,j) #ifdef coast # ifdef islands res(i,j)=res(i,j)*cof(i,j) # endif 395 continue #endif 400 continue #ifdef islands c c Find initial island residual using a line integral. Update RES(i,j) c on perimeter and interior of island. c do 430 isle=1,nisle istr=isis(isle) iend=ieis(isle) jstr=jsis(isle) jend=jeis(isle) resis=c0 do 410 j=jstr,jend do 410 i=istr,iend resis=resis+res(i,j) 410 continue do 420 j=jstr,jend do 420 i=istr,iend res(i,j)=resis 420 continue 430 continue #endif #ifdef cyclic c c Set Cyclic boundary conditions. c do 440 j=2,jmtm1 res(1 ,j)=res(imtm1,j) res(imt,j)=res(2 ,j) 440 continue #endif c c Initialize BETA and GDIR. c beta=c0 do 450 j=1,jmt do 450 i=1,imt gdir(i,j)=c0 450 continue c c----------------------------------------------------------------------- c Begin iteration loop using conjugate gradient technique. c----------------------------------------------------------------------- c mscan=0 460 continue mscan=mscan+1 c c----------------------------------------------------------------------- c Compute field of residuals. c----------------------------------------------------------------------- c do 470 j=2,jmtm1 #ifdef coast do 470 l=1,lseg istr=isz(j,l) if(istr.eq.0) goto 470 iend=iez(j,l) do 465 i=istr,iend #else do 470 i=2,imtm1 #endif gdir(i,j)=cfn(i,j)*res(i,j+1)+ * cfs(i,j)*res(i,j-1)+ * cfe(i,j)*res(i+1,j)+ * cfw(i,j)*res(i-1,j)- * res(i,j) #ifdef coast # ifdef islands gdir(i,j)=gdir(i,j)*cof(i,j) # endif 465 continue #endif 470 continue #ifdef islands c c Find island residuals using a line integral. Update GDIR on c perimeter and interior of island. c do 500 isle=1,nisle istr=isis(isle) iend=ieis(isle) jstr=jsis(isle) jend=jeis(isle) resis=c0 do 480 j=jstr,jend do 480 i=istr,iend resis=resis+gdir(i,j) 480 continue do 490 j=jstr,jend do 490 i=istr,iend gdir(i,j)=resis 490 continue 500 continue #endif #ifdef cyclic c c Set Cyclic boundary conditions. c do 510 j=2,jmtm1 gdir(1 ,j)=gdir(imtm1,j) gdir(imt,j)=gdir(2 ,j) 510 continue #endif c do 530 j=2,jmtm1 asum(j)=c0 #ifdef coast do 530 l=1,lseg istr=isz(j,l) if(istr.eq.0) goto 530 iend=iez(j,l) do 520 i=istr,iend #else do 520 i=2,imtm1 #endif asum(j)=asum(j)+gdir(i,j)*res(i,j) 520 continue 530 continue c apr4=c0 do 540 j=2,jmtm1 apr4=apr4+asum(j) 540 continue c if(mscan.gt.1) then if(oapr.ne.c0) then beta=apr4/oapr else beta=c0 endif endif c c----------------------------------------------------------------------- c Calculate new values for HDIR and FDIR. c----------------------------------------------------------------------- c do 555 j=2,jmtm1 #ifdef coast do 555 l=1,lseg istr=isz(j,l) if(istr.eq.0) goto 555 iend=iez(j,l) do 550 i=istr,iend #else do 550 i=2,imtm1 #endif hdir(i,j)=res(i,j)+beta*hdir(i,j) fdir(i,j)=gdir(i,j)+beta*fdir(i,j) 550 continue 555 continue c oapr=apr4 apr2=c0 do 570 j=2,jmtm1 asum(j)=c0 #ifdef coast do 570 l=1,lseg istr=isz(j,l) if(istr.eq.0) goto 570 iend=iez(j,l) do 560 i=istr,iend #else do 560 i=2,imtm1 #endif asum(j)=asum(j)+fdir(i,j)*fdir(i,j) 560 continue 570 continue c do 580 j=2,jmtm1 apr2=apr2+asum(j) 580 continue c c----------------------------------------------------------------------- c Calculate stepsize ASOR and make a correction to PTD and RES. c----------------------------------------------------------------------- c if(apr2.ne.c0) then asor=-apr4/apr2 else asor=c0 endif c do 595 j=2,jmtm1 #ifdef coast do 595 l=1,lseg istr=isz(j,l) if(istr.eq.0) goto 595 iend=iez(j,l) do 590 i=istr,iend #else do 590 i=2,imtm1 #endif ptd(i,j)=ptd(i,j)+asor*hdir(i,j) res(i,j)=res(i,j)+asor*fdir(i,j) 590 continue 595 continue #ifdef cyclic c c Set Cyclic boundary conditions. c do 600 j=2,jmtm1 res(1 ,j)=res(imtm1,j) res(imt,j)=res(2 ,j) ptd(1 ,j)=ptd(imtm1,j) ptd(imt,j)=ptd(2 ,j) 600 continue #endif c c----------------------------------------------------------------------- c Find the maximum absolute residual to determine convergence. c----------------------------------------------------------------------- c resmax=c0 do 615 j=2,jmtm1 #ifdef coast do 615 l=1,lseg istr=isz(j,l) if(istr.eq.0) goto 615 iend=iez(j,l) do 610 i=istr,iend #else do 610 i=2,imtm1 #endif resmax=max(abs(res(i,j)),resmax) 610 continue 615 continue c c----------------------------------------------------------------------- c Take another iteration unless the maximum of residual RESMAX is c greater or equal CRIT*PRMS or the number of scans reaches MXSCAN. c----------------------------------------------------------------------- c #ifndef vel_conv if((resmax.ge.crit*prms).and.(mscan.lt.mxscan)) go to 460 # else if( ((resmax.ge.crit*prms).or.(vbar_test(ptd,.false.).ge.crit)) * .and.(mscan.lt.mxscan)) go to 460 #endif c c======================================================================= c End of conjugate gradient iteration loop ============================= c======================================================================= #endif c c Report convergence information. c if(mod(itt,ntsi).eq.0) write(stdout,900) itt,mscan,resmax/prms, * resmax c c----------------------------------------------------------------------- c Update the stream function based upon the relaxation solution c----------------------------------------------------------------------- c if(.not.mxpas2) then do 620 j=1,jmt do 620 i=1,imt ptdb(i,j)=pb(i,j)+ptd(i,j) pb(i,j)=p(i,j) p(i,j)=ptdb(i,j) 620 continue else do 630 j=1,jmt do 630 i=1,imt p(i,j)=pb(i,j)+ptd(i,j) 630 continue endif c c----------------------------------------------------------------------- c Filter/smooth transport stream function field c----------------------------------------------------------------------- #ifdef shapiro c c If requested, apply Shapiro filter: c c ICNTP = counter (for frequency of application) c NFRQP = frequency with which filter is applied c NTIMP = number of times filter is applied per time step c NORDP = order of the filter c icntp=icntp+1 if(icntp.eq.nfrqp)then icntp=0 if(nordp.gt.0)then do 640 nn=1,ntimp call shap_lev(p,imt,jmt,tgrid,nordp) 640 continue do 650 j=1,jmt do 650 i=1,imt ptdb(i,j)=p(i,j) 650 continue endif endif #endif c c----------------------------------------------------------------------- c Save PTD to compute 1st guess for relaxation next timestep c (..note.. on 1st pass of Euler backward timestep, bypass this c save, since it will be done on the 2nd pass) c (..note.. on a mixing timestep, alter ptd to be consistent with c normal, leap-frog stepping) c----------------------------------------------------------------------- c if((mix.eq.1).and.eb) return if((mxp.ne.0).or.(mix.ne.0)) then do 660 j=1,jmt do 660 i=1,imt ptd(i,j)=c2*ptd(i,j) 660 continue endif call oput(kflds,nwds,(luptdb-1)*nwds+1,ptd) c return c 900 format(' RELAX - itt,mscan,ratio,resmax:',2(1x,i6),2(1x,1pe15.7)) c end