INCLUDE 'resolution.h' PARAMETER(IJ0=I0+J0,I1=I0-1,I2=I0-2,J1=J0-1,J2=J0-2, 1 K1=K0-1,K2=K0-2,NB1=NB0-1) c note on feb 2000 this has all of my 4th order modification c as well as new boundary forcing the latest of the best c most updates in 4th order corrections etc c still should compare to cray cctime5updated.f for sure about c latest on sponging bc's etc. june 7,1999 C C MAIN PROGRAM C Controls calculation, performs diagnostics, and saves graphics data. C Graphics and computer animation is done by postprocessors C C *** User-defined resolution parameters *** C Resolution (85X62X21) requires about 10 megabytes of array storage. C C Derived resolution parameters: C Horizontal resolution parameters are I0,I1,I2,I3, J0,J1,J2,J3 C Vertical resolution parameters are K0,K1,K2 C NB0,NB1 are dimension parameters for SEVP elliptic solver C C Double precision SEVP elliptic solver arrays c yes I run in double precision IMPLICIT REAL*8 (A-H,O-Z) c Lower precision logical depth and masking arrays INTEGER*2 KB,IU0,IV0,IN,IU,IV,IW C C "A" grid arrays COMMON/AGRID/ RHO(I0,J0,K1),PX(I0,J0),PY(I0,J0), 1 U1(I0,J0,K1),U2(I0,J0,K1),V1(I0,J0,K1),V2(I0,J0,K1), 2 S1(I0,J0,K1),S2(I0,J0,K1),T1(I0,J0,K1),T2(I0,J0,K1),P(I0,J0,K1), 3 ULF(I0,J0,K1),VLF(I0,J0,K1),SLF(I0,J0,K1),TLF(I0,J0,K1), 4 EV(I2,J2,K2),HV(I2,J2,K2),DAVG(K1),F(J2),TANPHI(J2),SUMIN(K1) C C "C" grid arrays COMMON/CGRID/U(I0,J0,K1),V(I0,J1,K1),W(I0,J0,K0) C User-defined scalar control parameters COMMON/CONTRO/ 1 DT,TLZ,B,G,DM0,DE0,DMZ0,RXYMX,RZMX,TAU,DRAG,FLTW,AV,DAODT,TAUDT, 2 TAUDTN,KTRM,M1,M2,M3,M4,M5,M5M,M6,M7,M8,KF,MVI,MXIT,ISAV,mstart, 3 MXSAV,LRSTRT,LOPEN,LSURF,LTURB,LWIND,LHEAT,LSOL,LSAVE,LMOVI C C Vertical grid arrays C "Z" array contains cell center AND interface depths COMMON/ZFS/Z(K0+K1),ODZ(K1),ODZW(K0) C C Logical depth and masking arrays C logical depth array "KB" is a user-defined array COMMON/TOPOG/KB(I0,J0),IU0(I0,J0),IV0(I0,J0), 1 IN(I0,J0,K1),IU(I0,J0,K1),IV(I0,J1,K1),IW(I0,J0,K0) C C SEVP elliptic solver arrays C SEVP subregion boundary array "IE" is a user-defined array COMMON/SEVP/RINV(I2,I2,NB0),RINV1(I2,I2,NB1),DUM0(I2,NB1), 1 DUM1(I2),DUM2(I2),X(I0,J0),H(I0,J0),AL(I2,J2),AB(I2,J2), 2 AC(I2,J2),AR(I2,J2),AT(I2,J2),S(I2,J2),IE(NB0), 3 ZBOTU(I1,J0),ZBOTV(I0,J1) C C Derived scalars COMMON/SCA/ODT,PRN,DAYS,ORXYMX,ORZMX,OFLTW,OM5,OM8,ITF,IT0 common/osave/ daodtsave,dtsave,odtsave,itfsave COMMON/METRICS/ 1 Y(J0),YV(J0),YVDEG(J0),YDEG(J0), 2 CS(J0),CSV(J1),OCS(J0),OCSV(J1), 3 DX(J0),DXV(J1),ODX(J0),ODXV(J1), 4 DY(J0),DYV(J1),ODY(J0),ODYV(J1) C Implicit vertical diffusion and Coriolis arrays COMMON/SURF/SSURF(I2,J2),TSURF(I2,J2) COMMON/WINDS/TAUX(I2,J2),TAUY(I2,J2) common/bdytsuv/ vnor(i0,k1), > uwst(j0,k1),vsud(i0,k1) CTS BOUNDARY CONDITION FOR NESTING COMMON/BCS/ 1 U2WST(J0,K1),V2WST(J0,K1),S2WST(J0,K1),T2WST(J0,K1), 3 U2SUD(I0,K1),V2SUD(I0,K1),S2SUD(I0,K1),T2SUD(I0,K1), 4 U2NOR(I0,K1),V2NOR(I0,K1),S2NOR(I0,K1),T2NOR(I0,K1), 5 UA2WST(J0,K1),VA2SUD(I0,K1),VA2NOR(I0,K1) COMMON/months/ > TSURFt(I2,J2,12),SSURFt(I2,J2,12), > TAUXt(I2,J2,12),TAUYt(I2,J2,12) COMMON/CALCUR/UWSTM(K1),VNORM(K1),VSUDM(K1) C common/leveldiv/ trant(k1),trann(k1),tranw(k1),trans(k1), > cor(k1) C Climate storing, scratch, and other auxillary arrays COMMON/KLIMAT/SCR(I0,J0,K1),SC1(IJ0),SC2(IJ0),TW(I0,J0,2), 1 TCLI(I0,J0,K0),UCLI(I0,J0),VCLI(I0,J0),RMSV(I0,J0), 2 XCLI(I0,J0),PCLI(I0,J0,2),HEATOP(I0,J0),SALTOP(I0,J0) C C Temporal data for diagnostics COMMON/TRNS/DTEMP(N0,K1),VKE(N0,K1),THST(N0),PXT(I0,J0) C Tide guage data (surface height anomaly at western boundary) COMMON/TIDES/ILOC(J0),WBDAYS(J2,360),WBSTPS(180) C Save CC cross-section of 18 zones nearest coast DIMENSION TBAR(18,K1),SBAR(18,K1),UBAR(18,K1),VBAR(18,K1), 1 WBAR(18,K1),SUMN(12),ISHORE(J0) C ********************************************************************** 5 FORMAT('LEVEL',I2,' DTEMP'/(8(1X,1PE9.2))) 6 FORMAT('P(I0/2,J0/2,2)'/(8(1X,1PE9.2))) 7 FORMAT('LEVEL',I2,' KE'/(8(1X,1PE9.2))) 8 FORMAT('MID-LAT STREAMFUNCTION PROFILE'/(8(1X,1PE9.2))) 9 FORMAT('MID-LAT TOP LAYER PRESSURE PROFILE'/(8(1X,1PE9.2))) 10 FORMAT(A5,I3,7E10.3/(120I1)) 12 FORMAT('*** UPWARD HEAT FLUXES (CALORIES/SEC) AT W-LEVEL ',I2 1 /'***',1X,1PE9.2,' CONVECTIVE'/'***',1X,1PE9.2,' DIFFUSIVE', 2 /'***',1X,1PE9.2,' CONTRIBUTION FROM TIME FILTERED W AND T') C C ---------------------- C OPEN OUTPUT DATA FILES C ---------------------- C Run history data OPEN(14,file='TR') C Restart data c istat OPEN(19,file='SV',form='unformatted') C Open 2-d cross-section files for DIECAST contour and vector graphics C in special subdirectories: C XY PLOTS: XYPLOTS/DATA (LOGICAL UNIT 30) C XZ PLOTS: XZPLOTS/DATA (LOGICAL UNIT 31) C YZ PLOTS: YZPLOTS/DATA (LOGICAL UNIT 32) OPEN(29,file='XYPLOTS/DEEP',form='unformatted') OPEN(30,file='XYPLOTS/DATA',form='unformatted') OPEN(31,file='XZPLOTS/DATA',form='unformatted') OPEN(32,file='YZPLOTS/DATA',form='unformatted') REWIND 14 REWIND 19 REWIND 29 REWIND 30 REWIND 31 REWIND 32 C C -------------------------------------------------- C INITIALIZE DERIVED SCALARS AND ALL NON-ZERO FIELDS C -------------------------------------------------- CALL INITFS idaodt = int(daodt) C Show initial conditions CALL FSPLTS(MXIT,1) CALL FLUSH(30) CALL FLUSH(31) CALL FLUSH(32) C Needed for change of time step on restart ITF=DAYS*DAODT C Set w-level for heat flux diagnostics (useful for GCM/climate studies) KF=2 write(6,*) ' itf,mxit ',itf,mxit IF (ITF.GE.MXIT) STOP 75 continue NSTPS=0 NDAYS=0 IT0=ITF C C ------------------------------ C MAIN TIME INTEGRATION LOOP 100 C ------------------------------ 100 IF (MOD(ITF,idaodt).NE.0) GO TO 15 CTS These data are for the end of the upcoming day CTS Read boundary values daily READ(91) TMP READ(91) ((UA2WST(J,K),J=2,J1),K=1,K1) READ(91) ((U2WST(J,K),J=2,J1),K=1,K1) READ(91) ((V2WST(J,K),J=2,J1),K=1,K1) READ(91) ((S2WST(J,K),J=2,J1),K=1,K1) READ(91) ((T2WST(J,K),J=2,J1),K=1,K1) READ(91) ((VA2SUD(I,K),I=2,I1),K=1,K1) READ(91) ((U2SUD(I,K),I=2,I1),K=1,K1) READ(91) ((V2SUD(I,K),I=2,I1),K=1,K1) READ(91) ((S2SUD(I,K),I=2,I1),K=1,K1) READ(91) ((T2SUD(I,K),I=2,I1),K=1,K1) READ(91) ((VA2NOR(I,K),I=2,I1),K=1,K1) READ(91) ((U2NOR(I,K),I=2,I1),K=1,K1) READ(91) ((V2NOR(I,K),I=2,I1),K=1,K1) READ(91) ((S2NOR(I,K),I=2,I1),K=1,K1) READ(91) ((T2NOR(I,K),I=2,I1),K=1,K1) CTS Convert to increments per time step CTS These are to be added at the end of each time step TMP=1./DAODT CTS DO 17 K=1,K1 DO 16 J=2,J1 UA2WST(J,K)=TMP*(UA2WST(J,K)-U(1,J,K)) U2WST(J,K)=TMP*(U2WST(J,K)-U1(1,J,K)) V2WST(J,K)=TMP*(V2WST(J,K)-V1(1,J,K)) S2WST(J,K)=TMP*(S2WST(J,K)-S1(1,J,K)) 16 T2WST(J,K)=TMP*(T2WST(J,K)-T1(1,J,K)) DO 17 I=2,I1 VA2SUD(I,K)=TMP*(VA2SUD(I,K)-V(I,1,K)) U2SUD(I,K)=TMP*(U2SUD(I,K)-U1(I,1,K)) V2SUD(I,K)=TMP*(V2SUD(I,K)-V1(I,1,K)) S2SUD(I,K)=TMP*(S2SUD(I,K)-S1(I,1,K)) T2SUD(I,K)=TMP*(T2SUD(I,K)-T1(I,1,K)) VA2NOR(I,K)=TMP*(VA2NOR(I,K)-V(I,J1,K)) U2NOR(I,K)=TMP*(U2NOR(I,K)-U1(I,J0,K)) V2NOR(I,K)=TMP*(V2NOR(I,K)-V1(I,J0,K)) S2NOR(I,K)=TMP*(S2NOR(I,K)-S1(I,J0,K)) 17 T2NOR(I,K)=TMP*(T2NOR(I,K)-T1(I,J0,K)) 15 ITF=ITF+1 DAYS=ITF/DAODT if(itf.eq.mstart+2) then days = 0. endif cTS now update forcing for model every day ' irfsh = 1*idaodt if(mod(itf,irfsh).eq.1) then if(itf.gt.1)call forcin(daodt,g,iday) write(6,*) ' made it back from forcing' C Calculate mean western boundary U DO 23 K=2,K1 TEMP=0. TMP1=0. UWSTM(K)=0. DO 22 J=2,J1 TMP2=IN(2,J,K)/ODY(J) TMP1=TMP1+TMP2 22 TEMP=TEMP+U1(1,J,K)*TMP2 IF (TMP1.NE.0.) TMP1=TEMP/TMP1 23 UWSTM(K)=TMP1 C Calculate mean southern boundary V DO 27 K=2,K1 TEMP=0. TMP1=0. VSUDM(K)=0. DO 25 I=2,I1 TMP2=IN(I,2,K)/ODXV(1) TMP1=TMP1+TMP2 25 TEMP=TEMP+V1(I,1,K)*TMP2 IF (TMP1.NE.0.) TMP1=TEMP/TMP1 27 VSUDM(K)=TMP1 c Calculate mean northern boundary v do 29 k=2,k1 TEMP=0. TMP1=0. VNORM(K)=0. DO 26 I=2,I1 TMP2=IN(I,J1,K)/ODXV(J1) TMP1=TMP1+TMP2 26 TEMP=TEMP+V1(I,J0,K)*TMP2 IF (TMP1.NE.0.) TMP1=TEMP/TMP1 VNORM(K)=TMP1 29 CONTINUE write(6,*) ' ** uwstm ',uwstm(1),uwstm(2) write(6,*) ' ** taux tauy ',taux(30,30),tauy(30,30) 7112 continue endif C C FS marches one time step each time it is called C =========== CALL FS CALL SETR(SCR,I0*J0*K1,0.) DO 800 K=1,K1 TMP=1./ODZ(K) DO 800 J=2,J1 DO 800 I=2,I1 SCR(I,J,2)=SCR(I,J,2)+TMP*IN(I,J,K)*U2(I,J,K) 800 SCR(I,J,3)=SCR(I,J,3)+TMP*IN(I,J,K)*V2(I,J,K) DO 805 J=2,J1 DO 805 I=2,I1 TMP=1./Z(2*KB(I,J)+1) SCR(I,J,2)=TMP*SCR(I,J,2) 805 SCR(I,J,3)=TMP*SCR(I,J,3) cTS save data every day to iyr idsav= 2*idaodt ldsav = idaodt CTS CTS if(mod(ITF,idaodt*10).eq.0)then CTS write(221) (((u2(i,j,k),i=1,i0),j=1,j0), CTS 1 k=1,k1) CTS write(222) (((v2(i,j,k),i=1,i0),j=1,j0), CTS 1 k=1,k1) CTS endif CTS C C Monitor solution progress in detail and check for possible instability C NOTE: DIECAST has been found to ALWAYS BE STABLE unless DT is too large C or user-specified i.c.s or b.c.s are UNPHYSICAL) CALL RANGE(V2,SC1,SC2,IN,I0,2,2,I1,J1,ILO,JLO,IHI,JHI,VLO,VHI) IF (ITF-IT0.GT.50.AND.MOD(ITF,idaodt).NE.0) GO TO 103 TMP=0. DO 101 J=2,J2 DO 101 I=2,I2 101 TMP=MAX(TMP,(1.-IU0(I,J))*ABS(U(I,J,1)), 1 (1.-IV0(I,J))*ABS(V(I,J,1))) TEMP=DAYS/360. WRITE(*,102) ITF,TEMP,ILO,JLO,VLO,IHI,JHI,VHI,TMP WRITE(14,102) ITF,TEMP,ILO,JLO,VLO,IHI,JHI,VHI,TMP 102 FORMAT('AT ITF=',I6,', YEAR',F7.2,/ 1 '(ILO,JLO,VLO,IHI,JHI,VHI)=(',2I4,1X,F9.4,2I4,1X,F9.4,')'/ 2 'max vel on land=',F8.4) CALL FLUSH(14) 103 IF ((VHI-VLO).LT.1000.) GO TO 105 WRITE(14,102) ITF,TEMP,ILO,JLO,VLO,IHI,JHI,VHI,TMP WRITE(14,104) ITF WRITE(6,104) ITF 104 FORMAT('STOP. UNIX-PECTEDLY LARGE VELOCITY AT ITF=',I6) CALL FSPLTS(MXIT,0) STOP C ------------------------ C COMPUTER-GENERAGED MOVIE C ------------------------ 105 IF (MOD(ITF,MVI).NE.0.OR.LMOVI.NE.1) GO TO 106 DO 1002 J=2,J1 DO 1002 I=2,I1 C SCR(I,J,1)=ODX(J)*.5*(U2(I,J,1)+U2(I-1,J,1)) C SCR(I,J,2)=ODY(J)*.5*(V2(I,J,1)+V2(I,J-1,1)) C SCR(I,J,3)=ODX(J)*.5*(U2(I,J,KTRM)+U2(I-1,J,KTRM)) C1002 SCR(I,J,4)=ODY(J)*.5*(V2(I,J,KTRM)+V2(I,J-1,KTRM)) C MOVIE module reads CELL CENTERED velocity !!! SCR(I,J,1)=ODX(J)*U2(I,J,1) SCR(I,J,2)=ODY(J)*V2(I,J,1) SCR(I,J,3)=ODX(J)*U2(I,J,KTRM) 1002 SCR(I,J,4)=ODY(J)*V2(I,J,KTRM) C to animate mean density above thermocline as well as top- and C thermocline- level pressure and geostrophic velocity and vorticity WRITE(*,1001) ITF 1001 FORMAT(/'===> Movie output at time step',I7) CALL COMXY(P,DAYS,IN) CALL COMXY(SCR,DAYS,IN) CALL COMXY(SCR(1,1,2),DAYS,IN) CALL COMXY(P(1,1,KTRM),DAYS,IN(1,1,KTRM)) CALL COMXY(SCR(1,1,3),DAYS,IN(1,1,KTRM)) CALL COMXY(SCR(1,1,4),DAYS,IN(1,1,KTRM)) C ---------------- C SPECIAL GRAPHICS C ---------------- 106 CONTINUE C106 DO 108 J=2,J1 C DO 108 I=2,I1 C TW(I,J,1)=OM8*(M8M*TW(I,J,1)+.5*(T2(I,J,KF-1)+T2(I,J,KF))) C108 TW(I,J,2)=OM8*(M8M*TW(I,J,2)+W(I,J,KF)) IF (MOD(ITF,M6).NE.0) GO TO 117 C Use DAVG as scratch for mean T DO 110 K=1,K1 DAVG(K)=0. DO 109 J=2,J1 DO 109 I=2,I1 109 DAVG(K)=DAVG(K)+T2(I,J,K)*IN(I,J,K) 110 DAVG(K)=SUMIN(K)*DAVG(K) WRITE(14,111) DAVG 111 FORMAT(' HORIZONTAL MEAN TEMPERATURE'/(6(1X,F9.4))) IF (ITF.NE.MXIT) GO TO 117 C Depth of specified isotherm ("TEMP" degrees) TEMP=15. TMP=1.E10 DO 114 J=2,J1 DO 114 I=2,I1 SCR(I,J,1)=0. DO 113 K=2,K1 IF (IN(I,J,K).EQ.0.OR.T2(I,J,K-1).LE.TEMP) GO TO 114 IF (T2(I,J,K).GT.TEMP) GO TO 113 SCR(I,J,1)=.01*(Z(2*K-2)+(Z(2*K)-Z(2*K-2))*(T2(I,J,K-1)-TEMP)/ 1 (T2(I,J,K-1)-T2(I,J,K))) TMP=MIN(TMP,SCR(I,J,1)) 113 CONTINUE 114 CONTINUE DO 115 J=2,J1 DO 115 I=2,I1 115 SCR(I,J,1)=MAX(SCR(I,J,1),TMP) CALL XYPLOT('D','depth of 15 deg. isotherm (m.)',SCR,1,1) THF=0. HFLT=0. HCND=0. DO 116 J=2,J1 DO 116 I=2,I1 SCR(I,J,4)=-W(I,J,KF)*(T2(I,J,KF-1)+T2(I,J,KF)) C SCR(I,J,5)=-TW(I,J,1)*TW(I,J,2) HCND=HCND+IW(I,J,KF)*DX(J)*DY(J)*(T2(I,J,KF)-T2(I,J,KF-1)) 1 *HV(I-1,J-1,KF-1) THF=THF+SCR(I,J,4)*DX(J)*DY(J) 116 HFLT=HFLT+SCR(I,J,5)*DX(J)*DY(J) THF=.5*THF HFLT=HFLT WRITE(14,12) KF,THF,HCND,HFLT IF (ITF.NE.MXIT) GO TO 117 C This provides proper 3rd subscript for TW C CALL XYPLOT('T','running time average T(deg. c)',TW,KF,KF) C CALL XYPLOT('W','running time average w(cm/sec)',TW(1,1,2),KF,KF) C CALL XYPLOT('h','upward heat flux ',SCR(1,1,4),KF,KF) C CALL XYPLOT('H','time average upward heat flux ',SCR(1,1,5),KF,KF) C 117 IF (MOD(ITF,M3).NE.0) GO TO 130 C Save zonal profile of latitudinally averaged top level pressure J=ITF/M3+1 JHST=MIN(J,J1) IF (JHST.EQ.J) GO TO 119 DO 118 J=2,J2 DO 118 I=2,I1 118 PXT(I,J)=PXT(I,J+1) 119 DO 121 I=2,I1 J=J0/2 JLO=MAX(2,J-5) JHI=MIN(J1,J+5) TMP=0. DO 120 J=JLO,JHI 120 TMP=TMP+P(I,J,1) 121 PXT(I,JHST)=TMP/(JHI-JLO+1.) 130 IF (MOD(ITF,M2).NE.0) GO TO 140 NHST=MOD(ITF/M2,N0)+1 THST(NHST)=P(I2/2,J2/2,1) DO 135 K=1,K1 VKE(NHST,K)=0. DTEMP(NHST,K)=0. DO 134 J=2,J1 DO 134 I=2,I1 VKE(NHST,K)=VKE(NHST,K)+.25*(U2(I-1,J,K)**2+U2(I,J,K)**2+ 1 V2(I,J-1,K)**2+V2(I,J,K)**2) 134 DTEMP(NHST,K)=DTEMP(NHST,K)+ABS(T2(I,J,K)-T1(I,J,K)) VKE(NHST,K)=VKE(NHST,K)*SUMIN(K) 135 DTEMP(NHST,K)=DTEMP(NHST,K)*SUMIN(K)/(TLZ*ODZ(K)) IF (NHST.NE.N0.AND.ITF.NE.MXIT) GO TO 140 DO 136 K=1,K1 WRITE(14,5) K,(DTEMP(I,K),I=1,N0) 136 WRITE(14,7) K,(VKE(I,K),I=1,N0) DO 138 I=1,NHST VKE(I,1)=SQRT(VKE(I,1))/ODZ(1) DO 137 K=2,K1 137 VKE(I,1)=VKE(I,1)+SQRT(VKE(I,K))/ODZ(K) 138 VKE(I,1)=SQRT(2.)*VKE(I,1)/TLZ WRITE(14,139) (VKE(I,1),I=1,NHST) 139 FORMAT(' RMS VELOCITY'/(6(1X,1PE9.2))) WRITE(14,6) THST 140 IF (MOD(ITF,M4).NE.0) GO TO 199 C C Integrate latitudinal velocity vertically and longitudinally to get C barotropic streamfunction. Evaluate at natural psi locations C (grid line intersects) CALL SETR(SCR,2*I0*J0,0.) DO 150 K=1,K1 TMP=1.E-12/ODZ(K) C SVR=TMP*DXV(1) C DO 148 I=2,I1 C148 SCR(I,1,1)=SCR(I-1,1,1)+V(I,1,K)*SVR DO 150 J=1,J1 SVR=TMP*DXV(J) DO 150 I=2,I1 150 SCR(I,J,2)=SCR(I,J,2)+V(I,J,K)*SVR C DO 152 I=I2,1,-1 DO 152 J=1,J1 152 SCR(I,J,1)=SCR(I+1,J,1)-SCR(I+1,J,2) C C ------------------------------- C MEAN (TIME AVERAGE) DIAGNOSTICS C ------------------------------- C STREAMFUNCTION (XCLI(I,J)) AV=AV+1. WT=1./AV OWT=1.-WT DO 158 J=1,J1 DO 158 I=1,I1 C Long term mean (climate) 158 XCLI(I,J)=OWT*XCLI(I,J)+WT*SCR(I,J,1) C Moving window mean C158 XCLI(I,J)=M5M*OM5*XCLI(I,J)+OM5*SCR(I,J,1) K=0 159 K=K+1 IF (Z(2*K).LT.7.E4) GO TO 159 KL1=K-1 KL2=K WT2=(7.E4-Z(2*KL1))/(Z(2*KL2)-Z(2*KL1)) WT1=1.-WT2 DO 160 J=2,J1 DO 160 I=2,I1 TMP= WT1*U2(I,J,KL1)+WT2*U2(I,J,KL2) TEMP=WT1*V2(I,J,KL1)+WT2*V2(I,J,KL2) C Climate mean U and V at depth 700 meters UCLI(I,J)=OWT*UCLI(I,J)+WT*TMP VCLI(I,J)=OWT*VCLI(I,J)+WT*TEMP RMSV(I,J)=OWT*RMSV(I,J)+WT*((UCLI(I,J)-TMP)**2+(VCLI(I,J)-TMP)**2) C Mean surface pressure PCLI(I,J,1)=OWT*PCLI(I,J,1)+WT*P(I,J,1) C Mean squared deviation of surface pressure from mean 160 PCLI(I,J,2)=OWT*PCLI(I,J,2)+WT*(P(I,J,1)-PCLI(I,J,1))**2 C Mean temperature DO 162 K=1,K1 DO 162 J=2,J1 DO 162 I=2,I1 C Mean temperature (long term mean) 162 TCLI(I,J,K)=OWT*TCLI(I,J,K)+WT*T2(I,J,K) 199 IF (MOD(ITF,M6).NE.0) GO TO 220 C Interpolate streamfunction to "P" cell centers for plotting DO 200 J=1,J2 DO 200 I=1,I2 SCR(I+1,J+1,2)= 1 .25*(SCR(I,J,1)+SCR(I+1,J,1)+SCR(I,J+1,1)+SCR(I+1,J+1,1)) 200 SCR(I+1,J+1,3)= 1 .25*(XCLI(I,J)+XCLI(I+1,J)+XCLI(I,J+1)+XCLI(I+1,J+1)) CALL XYPLOT('X','streamfunction (Sverdrups) ',SCR(1,1,2),1,1) CALL XYPLOT('X',' long term avg xi (Sverdrups) ',SCR(1,1,3),1,1) UMAX=0. VMAX=0. DO 210 J=2,J1 DO 210 I=2,I1 UMAX=MAX(UMAX,ABS(SCR(I,J,1)-SCR(I,J-1,1)))*ODY(J) 210 VMAX=MAX(VMAX,ABS(SCR(I,J,1)-SCR(I-1,J,1)))*ODXV(J) UMAX=1.E12*UMAX/TLZ VMAX=1.E12*VMAX/TLZ WRITE(14,211) UMAX,VMAX 211 FORMAT(' BAROTROPIC (UMAX,VMAX)= ','(',1PE9.2,',',1PE9.2,')') C Calculate top layer vorticity at natural (grid line intersect) points C and divergence at P-points CALL SETR(SCR,3*I0*J0,0.) TMP=7.*24.*3600. DO 217 J=2,J2 DO 217 I=2,I2 217 SCR(I,J,2)= 1 (V(I+1,J,1)-V(I,J,1))*ODX(J)-(U(I,J+1,1)-U(I,J,1))*ODYV(J) DO 218 J=1,J2 DO 218 I=1,I2 C Interpolate vorticity to P-points 218 SCR(I+1,J+1,1)=.25*TMP* 1 (SCR(I,J,2)+SCR(I+1,J,2)+SCR(I,J+1,2)+SCR(I+1,J+1,2)) CALL XYPLOT('s','vorticity (per week) ',SCR,1,1) 220 IF (ITF.NE.MXIT) GO TO 250 TMP=7.*24.*3600. C C -------------------------- C SPECIAL LONG TERM AVERAGES C -------------------------- C Conversion to deg C per year and g/kg/yr TEMP=52.*TMP/DT C Conversion to watts/m-m WATTS=TEMP*0.00134/ODZ(1) DO 225 J=2,J1 DO 225 I=2,I1 SCR(I,J,1)=PXT(I,J)/980. SCR(I,J,2)=WATTS*HEATOP(I,J) SCR(I,J,3)=TEMP*SALTOP(I,J) SCR(I,J,4)=SQRT(RMSV(I,J)) 225 SCR(I,J,5)=TMP*((U(I,J,2)-U(I,J,2))*ODX(J) 1 +(V(I,J,2)-V(I,J,2))*ODY(J)) CALL XYPLOT('d','horizontal divergence(per wk.)',SCR(1,1,5),1,1) C CALL XYPLOT('p','lat-avg. surf. p (eq. fsa, cm)',SCR,1,1) CALL XYPLOT('q','mean Tsurf restoring(watts/m2)',SCR(1,1,2),1,1) CALL XYPLOT('n','mean Ssurf restoring (g/kg/yr)',SCR(1,1,3),1,1) CALL XYPLOT('c','rms v anomaly (cm/sec) at 700m',SCR(1,1,4),1,1) CALL XYPLOT('U','mean u (cm/sec) at depth 700m ',UCLI,1,1) CALL XYPLOT('V','mean v (cm/sec) at depth 700m ',VCLI,1,1) C Calculate mean zonal heating rates SUMLAT=0. SCR(3,1,1)=0. DO 240 J=2,J1 SUM=0. INSUM=0 DO 228 I=2,I1 INSUM=INSUM+IN(I,J,1) 228 SUM=SUM+HEATOP(I,J)*IN(I,J,1) C Zonal mean sensible heating rate (surface restoring) in deg C per year SCR(1,J,1)=TEMP*SUM/INSUM C Convert to zonal mean watts/m2 SCR(2,J,1)=0.00134*SCR(1,J,1)/ODZ(1) C C Sum total SURFACE sensible heat inputs over latitude in terawatts C When added to the latitudinal sponge layer restorings, this matches the C model's equilibrium mean advective latitudinal heat transport C C Area must be in m-m, because SCR(2,J,1) is watts/m-m AREA=1.E-4*INSUM/(ODX(J)*ODY(J)) SUMLAT=SUMLAT+1.E-12*AREA*SCR(2,J,1) 240 SCR(3,J,1)=SUMLAT WRITE(14,241) (SCR(1,J,1),J=2,J1) 241 FORMAT(/'Lat. profile of zonal mean sensible heating rate (surface 1 restoring deg. C/year)'/(10F7.2)) WRITE(14,242) (SCR(2,J,1),J=2,J1) 242 FORMAT(/'Lat. profile of zonal mean sensible heating rate (surface 1 restoring watts/m-m)'/(10F7.2)) WRITE(14,243) (SCR(3,J,1),J=2,J1) 243 FORMAT(/'Lat. profile of latitudinally integrated sensible heating 1 rate (terawatts)',1P,/(10(1X,E8.1))) 250 IF (ITF.NE.1.AND.MOD(ITF,80).NE.0) GO TO 300 goto 300 CALL SETR(SBAR,18*K1,0.) CALL SETR(TBAR,18*K1,0.) CALL SETR(UBAR,18*K1,0.) CALL SETR(VBAR,18*K1,0.) CALL SETR(WBAR,18*K1,0.) DO 260 J=2,J1 N=1 DO 259 I=2,I1 259 N=N+IN(I,J,1) C Last wet point is I=ISHORE(J) 260 ISHORE(J)=N DO 280 K=1,K1 CALL SETR(SUMN,18,0.) DO 270 J=31,j1 IHI=ISHORE(J) write(6,*) ' ihi ',ihi N=0 DO 270 I=IHI-17,IHI N=N+1 SUMN(N)=SUMN(N)+IN(I,J,K) SBAR(N,K)=SBAR(N,K)+IN(I,J,K)*S2(I,J,K) TBAR(N,K)=TBAR(N,K)+IN(I,J,K)*T2(I,J,K) UBAR(N,K)=UBAR(N,K)+IN(I,J,K)*U2(I,J,K) VBAR(N,K)=VBAR(N,K)+IN(I,J,K)*V2(I,J,K) 270 WBAR(N,K)=WBAR(N,K)+IN(I,J,K)*.5*(W(I,J,K)+W(I,J,K+1)) DO 280 N=1,18 IF (SUMN(N).EQ.0) GO TO 280 SBAR(N,K)=SBAR(N,K)/SUMN(N) TBAR(N,K)=TBAR(N,K)/SUMN(N) UBAR(N,K)=UBAR(N,K)/SUMN(N) VBAR(N,K)=VBAR(N,K)/SUMN(N) WBAR(N,K)=WBAR(N,K)/SUMN(N) 280 CONTINUE WRITE(15,289) DAYS WRITE(*,289) DAYS 289 FORMAT('day',F7.2) WRITE(15,290) (.01*Z(2*K),(SBAR(N,K),N=1,18),K=1,K1) WRITE(*,290) (.01*Z(2*K),(SBAR(N,K),N=1,18),K=1,K1) 290 FORMAT( 1 /'latitudinal mean S cross-section to 18 points from coast' 2 /(F6.0,18F6.2)) WRITE(15,291) (.01*Z(2*K),(TBAR(N,K),N=1,18),K=1,K1) WRITE(*,291) (.01*Z(2*K),(TBAR(N,K),N=1,18),K=1,K1) 291 FORMAT( 1 /'latitudinal mean T cross-section to 18 points from coast' 2 /(F6.0,18F6.2)) WRITE(15,292) (.01*Z(2*K),(UBAR(N,K),N=1,18),K=1,K1) WRITE(*,292) (.01*Z(2*K),(UBAR(N,K),N=1,18),K=1,K1) 292 FORMAT( 1 /'latitudinal mean U cross-section to 18 points from coast' 2 /(F6.0,18F6.2)) WRITE(15,293) (.01*Z(2*K),(VBAR(N,K),N=1,18),K=1,K1) WRITE(*,293) (.01*Z(2*K),(VBAR(N,K),N=1,18),K=1,K1) 293 FORMAT( 1 /'latitudinal mean V cross-section to 18 points from coast' 2 /(F6.0,18F6.2)) WRITE(15,294) (.01*Z(2*K),(100.*WBAR(N,K),N=1,18),K=1,K1) WRITE(*,294) (.01*Z(2*K),(100.*WBAR(N,K),N=1,18),K=1,K1) 294 FORMAT( 1 /'100 * latitudinal mean W cross-section to 18 points from coast' 2 /(F6.0,18F6.2)) 300 IF (MOD(ITF,ISAV).NE.0.AND.ITF.NE.MXIT) GO TO 410 ISAV=MIN0(2*ISAV,MXSAV) C Flag to do plots at all levels only at every 24th save C Otherwise, top level only I=MOD(ITF/ISAV,24) C C ---------------------------------------------------------------------- C Contour and vector plots: C save data for post-process screen viewing and optional hardcopy C ---------------------------------------------------------------------- c CALL FSPLTS(MXIT,I) call fsplts(mxit,1) C C Save 3-d T data for empirical orthogonal function analysis IF (LSAVE.NE.1) GO TO 410 C C Save restart data 403 REWIND 19 WRITE(19) ITF,DAYS,AV,X,U1,U2,V1,V2,S1,S2,T1,T2,P,ULF,VLF,SLF,TLF, 1 U,V,W,EV,HV,TW,TCLI,UCLI,VCLI,RMSV,HEATOP,SALTOP,XCLI,PCLI,PXT, 2 VNORM,UWSTM,VSUDM NXTSAV=ITF+M2 WRITE(14,404) ITF WRITE(*,404) ITF 404 FORMAT(' SAVE AT ITF = ',I6) 410 IF (MOD(ITF,10).NE.0) GO TO 420 416 FORMAT(I6) 420 IF (MOD(ITF,M2).EQ.0) WRITE(14,425) ITF,NXTSAV 425 FORMAT(' ITF=',I6,', next save at ITF=',I6) IF (I0.NE.0) GO TO 450 DLAPS=(ITF-IT0)/DAODT IF (DLAPS.LT.360) GO TO 450 NSTPS=NSTPS+1 I=ILOC(J0/2) WBSTPS(NSTPS)=(P(I,J0/2,1)-PCLI(I,J0/2,1))/980. C IF (NSTPS.EQ.180) WRITE(16,430) WBSTPS 430 FORMAT(10F8.3) IF (NSTPS.EQ.180) NSTPS=0 I=DAODT IF (MOD(ITF,I).NE.0) GO TO 450 NDAYS=NDAYS+1 DO 440 J=2,J1 I=ILOC(J) 440 WBDAYS(J-1,NDAYS)=(P(I,J,1)-PCLI(I,J,1))/980. C IF (NDAYS.EQ.360) WRITE(15,430) WBDAYS IF (NDAYS.EQ.360) NDAYS=0 450 IF (ITF.LT.MXIT) GO TO 100 END ************************************************** C MAIN DYNAMICS PROGRAM AND ASSOCIATED SUBROUTINES ************************************************** SUBROUTINE FS C ---------------------------------------------------------------------- C FS is a wind stress and buoyancy driven free stream model, with C rigid top, variable DZ, and bottom topography. C FS resolves bottom features using a staircase approximation. C EXACT conservation is satisfied. C Bottom fluxes are specified as SOURCES when coupled to TS (thin-shell C submodel). W is positive downward. W(I,J,1) is at free surface. C W(I,J,KB(I,J)+1), where 0.LE.KB.LE.K1, is at bottom C ---------------------------------------------------------------------- INCLUDE 'resolution.h' PARAMETER(IJ0=I0+J0,I1=I0-1,I2=I0-2,I3=I0-3,J1=J0-1, 1 J2=J0-2,J3=J0-3,K1=K0-1,K2=K0-2,NB1=NB0-1) IMPLICIT REAL*8 (A-H,O-Z) INTEGER*2 KB,IU0,IV0,IN,IU,IV,IW,ileft,irite COMMON/AGRID/ RHO(I0,J0,K1),PX(I0,J0),PY(I0,J0), 1 U1(I0,J0,K1),U2(I0,J0,K1),V1(I0,J0,K1),V2(I0,J0,K1), 2 S1(I0,J0,K1),S2(I0,J0,K1),T1(I0,J0,K1),T2(I0,J0,K1),P(I0,J0,K1), 3 ULF(I0,J0,K1),VLF(I0,J0,K1),SLF(I0,J0,K1),TLF(I0,J0,K1), 4 EV(I2,J2,K2),HV(I2,J2,K2),DAVG(K1),F(J2),TANPHI(J2),SUMIN(K1) COMMON/CGRID/U(I0,J0,K1),V(I0,J1,K1),W(I0,J0,K0) COMMON/OPENB/AROUT COMMON/CONTRO/ 1 DT,TLZ,B,G,DM0,DE0,DMZ0,RXYMX,RZMX,TAU,DRAG,FLTW,AV,DAODT,TAUDT, 2 TAUDTN,KTRM,M1,M2,M3,M4,M5,M5M,M6,M7,M8,KF,MVI,MXIT,ISAV,mstart, 3 MXSAV,LRSTRT,LOPEN,LSURF,LTURB,LWIND,LHEAT,LSOL,LSAVE,LMOVI COMMON/SPONGE/DM0X(I0,j0) common/bloks/ileft(j0,k1),irite(j0,k1) COMMON/FLX/ 1 UX(I1,J2),UY(I2,J1),UZ(I2,J2,2),VX(I1,J2),VY(I2,J1),VZ(I2,J2,2), 2 SX(I1,J2),SY(I2,J1),SZ(I2,J2,2),TX(I1,J2),TY(I2,J1),TZ(I2,J2,2) COMMON/TRNS/DTEMP(N0,K1),VKE(N0,K1),THST(N0),PXT(I0,J0) COMMON/SURF/SSURF(I2,J2),TSURF(I2,J2) COMMON/KLIMAT/SCR(I0,J0,K1),SC1(IJ0),SC2(IJ0),TW(I0,J0,2), 1 TCLI(I0,J0,K0),UCLI(I0,J0),VCLI(I0,J0),RMSV(I0,J0), 2 XCLI(I0,J0),PCLI(I0,J0,2),HEATOP(I0,J0),SALTOP(I0,J0) COMMON/PUNIX/UBT(I1,J0),VBT(I0,J1) C Northern and southern sponge layer Levitus Climatology COMMON/SCA/ODT,PRN,DAYS,ORXYMX,ORZMX,OFLTW,OM5,OM8,ITF,IT0 CTS NESTING BOUNDARY COMMON/BCS/ 1 U2WST(J0,K1),V2WST(J0,K1),S2WST(J0,K1),T2WST(J0,K1), 3 U2SUD(I0,K1),V2SUD(I0,K1),S2SUD(I0,K1),T2SUD(I0,K1), 4 U2NOR(I0,K1),V2NOR(I0,K1),S2NOR(I0,K1),T2NOR(I0,K1), 5 UA2WST(J0,K1),VA2SUD(I0,K1),VA2NOR(I0,K1) COMMON/METRICS/ 1 Y(J0),YV(J0),YVDEG(J0),YDEG(J0), 2 CS(J0),CSV(J1),OCS(J0),OCSV(J1), 3 DX(J0),DXV(J1),ODX(J0),ODXV(J1), 4 DY(J0),DYV(J1),ODY(J0),ODYV(J1) COMMON/SEVP/RINV(I2,I2,NB0),RINV1(I2,I2,NB1),DUM0(I2,NB1), 1 DUM1(I2),DUM2(I2),X(I0,J0),H(I0,J0),AL(I2,J2),AB(I2,J2), 2 AC(I2,J2),AR(I2,J2),AT(I2,J2),S(I2,J2),IE(NB0), 3 ZBOTU(I1,J0),ZBOTV(I0,J1) COMMON/TOPOG/KB(I0,J0),IU0(I0,J0),IV0(I0,J0), 1 IN(I0,J0,K1),IU(I0,J0,K1),IV(I0,J1,K1),IW(I0,J0,K0) C Z IS DISTANCE FROM TOP OCEAN BOUNDARY, INCREASING WITH INCREASING K COMMON/ZFS/Z(K0+K1),ODZ(K1),ODZW(K0) common/bdytsuv/ vnor(i0,k1), > uwst(j0,k1),vsud(i0,k1) COMMON/CALCUR/UWSTM(K1),VNORM(K1),VSUDM(K1) O12 = 1./12. O24 = 1./24. C ====================================================== C Hydrostatic equation with static balance subracted out C ====================================================== C For better rho diagnostics TEMP=1.E5 DO 96 K=1,K1 C For minimum roundoff error subtract each layer mean C TEMP=1.E5 DO 96 J=2,J1 ILO = ILEFT(J,K) IHI = IRITE(J,K) DO 96 I=ILO,IHI CTS DO 96 I=2,I1 TMP=T2(I,J,K) SLT=S2(I,J,K) c dietrich orginial eos RLAM=1727.0849+12.2662*TMP-.045253*TMP**2-(3.6305+.036811*TMP)*SLT P01= 5832.294+42.8068*TMP-.31542*TMP**2+3.41777*SLT 1 +.00041517*TMP**3-SLT*TMP*.088353 RHO(I,J,K)=P01/(RLAM+.703845*P01) 96 CONTINUE DO 99 K=1,K2 L=K+1 TMP=.5*G/ODZW(L) DO 99 J=2,J1 DO 99 I=2,I1 99 P(I,J,L)=P(I,J,K)+TMP*(RHO(I,J,K)+RHO(I,J,L)) C Variable horizontal diffusivity arrays DMX and DMY are inactivated C to reduce storage in this version C Set moving window indices and top level arrays 202 LB=1 LT=2 DO 250 J=1,J2 DO 250 I=1,I2 UZ(I,J,1)=0. VZ(I,J,1)=0. SZ(I,J,1)=0. 250 TZ(I,J,1)=0. C ================================= C CALCULATE ALL INTERNAL FLUXES and C SUBSTITUTE INTO CONSERVATION LAWS C ================================= C Level-by-level moving window reduces storage. C Note: all BOUNDARY fluxes are masked out in this loop, C but are included later. This is the main computation loop. CTS Remove the parallel setup: Some wrong index problems CTS and it is not very efficient DO 500 K=1,K1 L=K+1 C --------------------------- C 4th order pressure gradient C --------------------------- C First pressure differences DO 300 J=2,J1 DO 300 I=2,I1 300 UX(I,J-1)=IU(I,J,K)*(P(I+1,J,K)-P(I,J,K)) DO 301 J=2,J2 DO 301 I=2,I1 301 VY(I-1,J)=IV(I,J,K)*(P(I,J+1,K)-P(I,J,K)) C Zero p.g. over land IF (K.NE.1) GO TO 304 DO 302 J=2,J1 DO 302 I=2,I1 302 UX(I,J-1)=IU0(I,J)*UX(I,J-1) DO 303 J=2,J2 DO 303 I=2,I1 303 VY(I-1,J)=IV0(I,J)*VY(I-1,J) 304 DO 306 J=2,J1 C Basic second order DO 305 I=2,I1 PX(I,J)=6.*(UX(I,J-1)+UX(I-1,J-1)) 305 PY(I,J)=6.*(VY(I-1,J)+VY(I-1,J-1)) C Fourth order correction DO 306 I=3,I2 C Dan Wright's suggestion: less noisy 306 PX(I,J)=PX(I,J)+UX(I,J-1)+UX(I-1,J-1)-UX(I+1,J-1)-UX(I-2,J-1) DO 307 J=3,J2 DO 307 I=2,I1 307 PY(I,J)=PY(I,J)+VY(I-1,J)+VY(I-1,J-1)-VY(I-1,J+1)-VY(I-1,J-2) DO 308 J=2,J1 DO 308 I=2,I1 PX(I,J)=O12*PX(I,J) 308 PY(I,J)=O12*PY(I,J) C --------------- C Vertical fluxes C --------------- IF (K.NE.K1) GO TO 317 C Set bottom fluxes to zero (drag law applied after label 500) DO 310 J=1,J2 do 310 i=1,I2 UZ(I,J,LT)=0. VZ(I,J,LT)=0. SZ(I,J,LT)=0. 310 TZ(I,J,LT)=0. GO TO 359 c alternate 4th order vertical fluxes 317 DO 341 J=2,J1 do 341 i=2,I1 SCR(I,J,1)=6.*(U2(I,J,K)+U2(I,J,L)) SCR(I,J,2)=6.*(V2(I,J,K)+V2(I,J,L)) SCR(I,J,3)=6.*(S2(I,J,K)+S2(I,J,L)) 341 SCR(I,J,4)=6.*(T2(I,J,K)+T2(I,J,L)) C Already skipped if K=K1 IF (K.EQ.1.OR.K.EQ.K2) GO TO 343 DO 342 J=2,J1 do 342 i=2,I1 TMP=IW(I,J,K)*IW(I,J,K+2) SCR(I,J,1)=SCR(I,J,1)+TMP* 1 (-ULF(I,J,K-1)+U2(I,J,K)+U2(I,J,L)-U2(I,J,K+2)) SCR(I,J,2)=SCR(I,J,2)+TMP* 1 (-VLF(I,J,K-1)+V2(I,J,K)+V2(I,J,L)-V2(I,J,K+2)) SCR(I,J,3)=SCR(I,J,3)+TMP* 1 (-SLF(I,J,K-1)+S2(I,J,K)+S2(I,J,L)-S2(I,J,K+2)) 342 SCR(I,J,4)=SCR(I,J,4)+TMP* 1 (-TLF(I,J,K-1)+T2(I,J,K)+T2(I,J,L)-T2(I,J,K+2)) 343 DO 344 N=1,4 DO 344 J=2,J1 do 344 i=1,I1 344 SCR(I,J,N)=O12*SCR(I,J,N) DO 349 J=2,J1 do 349 i=2,I1 UZ(I-1,J-1,LT)=W(I,J,L)*SCR(I,J,1) 1 -EV(I-1,J-1,K)*(U1(I,J,L)-U1(I,J,K))*IW(I,J,L) VZ(I-1,J-1,LT)=W(I,J,L)*SCR(I,J,2) 1 -EV(I-1,J-1,K)*(V1(I,J,L)-V1(I,J,K))*IW(I,J,L) SZ(I-1,J-1,LT)=W(I,J,L)*SCR(I,J,3) 1 -HV(I-1,J-1,K)*(S1(I,J,L)-S1(I,J,K))*IW(I,J,L) 349 TZ(I-1,J-1,LT)=W(I,J,L)*SCR(I,J,4) 1 -HV(I-1,J-1,K)*(T1(I,J,L)-T1(I,J,K))*IW(I,J,L) C ------------------- C Longitudinal fluxes C ------------------- 359 DO 370 J=2,J1 C Higher order interpolations DO 361 I=1,I1 SCR(I,J,1)=6.*(U2(I,J,K)+U2(I+1,J,K)) SCR(I,J,2)=6.*(V2(I,J,K)+V2(I+1,J,K)) SCR(I,J,3)=6.*(S2(I,J,K)+S2(I+1,J,K)) 361 SCR(I,J,4)=6.*(T2(I,J,K)+T2(I+1,J,K)) choose either old scheme or my scheme careful with comments!! check in Latidunal loop also c goto 1199 do 1362 i=3,i2 c do 1362 i=2,i2 scr(i,j,1) = scr(i,j,1) + iu(i-1,j,k)*(u2(i,j,k)- > u2(i-1,j,k))+iu(i+1,j,k)*(u2(i+1,j,k)-u2(i+2,j,k)) > +(1-iu(i-1,j,k))*u2(i,j,k)+1.*(1-iu(i+1,j,k))*u2(i+1,j,k) c can make above coast condition symetric or zeroi, 2 or 1. scr(i,j,2) = scr(i,j,2) + iu(i-1,j,k)*(v2(i,j,k)- > v2(i-1,j,k))+iu(i+1,j,k)*(v2(i+1,j,k)-v2(i+2,j,k)) scr(i,j,3) = scr(i,j,3) + iu(i-1,j,k)*(s2(i,j,k)- > s2(i-1,j,k))+iu(i+1,j,k)*(s2(i+1,j,k)-s2(i+2,j,k)) 1362 scr(i,j,4) = scr(i,j,4) + iu(i-1,j,k)*(t2(i,j,k)- > t2(i-1,j,k))+iu(i+1,j,k)*(t2(i+1,j,k)-t2(i+2,j,k)) goto 363 1199 continue C Skip loop 362 for original 2nd order scheme DO 362 I=3,I3 TMP=IN(I-1,J,K)*IN(I,J,K)*IN(I+1,J,K)*IN(I+2,J,K) SCR(I,J,1)=SCR(I,J,1)+TMP* 1 (-U2(I-1,J,K)+U2(I,J,K)+U2(I+1,J,K)-U2(I+2,J,K)) SCR(I,J,2)=SCR(I,J,2)+TMP* 1 (-V2(I-1,J,K)+V2(I,J,K)+V2(I+1,J,K)-V2(I+2,J,K)) SCR(I,J,3)=SCR(I,J,3)+TMP* 1 (-S2(I-1,J,K)+S2(I,J,K)+S2(I+1,J,K)-S2(I+2,J,K)) 362 SCR(I,J,4)=SCR(I,J,4)+TMP* 1 (-T2(I-1,J,K)+T2(I,J,K)+T2(I+1,J,K)-T2(I+2,J,K)) 363 DO 364 N=1,4 DO 364 I=1,I1 c DO 364 i=1,I2 364 SCR(I,J,N)=O12*SCR(I,J,N) DO 370 I=1,I1 c DO 370 i=1,I2 AM=DM0X(I,j)*ODX(J) AH=AM C UX could be simply U(I,J,K)**2 + ...., but this does not conserve C energy and is unstable UX(I,J-1)=(U(I,J,K)*SCR(I,J,1) 1 -AM*(U1(I+1,J,K)-U1(I,J,K)))*IU(I,J,K) VX(I,J-1)=(U(I,J,K)*SCR(I,J,2) 1 -AM*(V1(I+1,J,K)-V1(I,J,K)))*IU(I,J,K) SX(I,J-1)=(U(I,J,K)*SCR(I,J,3) 1 -AH*(S1(I+1,J,K)-S1(I,J,K)))*IU(I,J,K) 370 TX(I,J-1)=(U(I,J,K)*SCR(I,J,4) 1 -AH*(T1(I+1,J,K)-T1(I,J,K)))*IU(I,J,K) C ------------------ C Latitudinal fluxes C ------------------ DO 380 J=1,J1 DO 371 I=2,I1 SCR(I,J,1)=6.*(U2(I,J,K)+U2(I,J+1,K)) SCR(I,J,2)=6.*(V2(I,J,K)+V2(I,J+1,K)) SCR(I,J,3)=6.*(S2(I,J,K)+S2(I,J+1,K)) 371 SCR(I,J,4)=6.*(T2(I,J,K)+T2(I,J+1,K)) C IF (I0.NE.0) GO TO 373 Cdietrich Skip loop 372 for original 2nd order scheme IF (J.EQ.1.OR.J.EQ.J1) GO TO 373 c comment for revised 4th order, else uses old second order c goto 1399 do 1372 i=2,i1 scr(i,j,1) = scr(i,j,1) + iv(i,j-1,k)*(u2(i,j,k)- > u2(i,j-1,k))+iv(i,j+1,k)*(u2(i,j+1,k)-u2(i,j+2,k)) scr(i,j,2) = scr(i,j,2) + iv(i,j-1,k)*(v2(i,j,k)- > v2(i,j-1,k))+iv(i,j+1,k)*(v2(i,j+1,k)-v2(i,j+2,k)) scr(i,j,3) = scr(i,j,3) + iv(i,j-1,k)*(s2(i,j,k)- > s2(i,j-1,k))+iv(i,j+1,k)*(s2(i,j+1,k)-s2(i,j+2,k)) 1372 scr(i,j,4) = scr(i,j,4) + iv(i,j-1,k)*(t2(i,j,k)- > t2(i,j-1,k))+iv(i,j+1,k)*(t2(i,j+1,k)-t2(i,j+2,k)) goto 373 1399 continue DO 372 I=2,I1 TMP=IN(I,J-1,K)*IN(I,J,K)*IN(I,J+1,K)*IN(I,J+2,K) SCR(I,J,1)=SCR(I,J,1)+TMP* 1 (-U2(I,J-1,K)+U2(I,J,K)+U2(I,J+1,K)-U2(I,J+2,K)) SCR(I,J,2)=SCR(I,J,2)+TMP* 1 (-V2(I,J-1,K)+V2(I,J,K)+V2(I,J+1,K)-V2(I,J+2,K)) SCR(I,J,3)=SCR(I,J,3)+TMP* 1 (-S2(I,J-1,K)+S2(I,J,K)+S2(I,J+1,K)-S2(I,J+2,K)) 372 SCR(I,J,4)=SCR(I,J,4)+TMP* 1 (-T2(I,J-1,K)+T2(I,J,K)+T2(I,J+1,K)-T2(I,J+2,K)) 373 DO 374 N=1,4 DO 374 I=2,I1 C AM=DM0*ODYV(J) C AH=DE0*ODYV(J) 374 SCR(I,J,N)=O12*SCR(I,J,N) DO 380 I=2,I1 AM=DM0X(I,j)*ODX(J) AH=AM UY(I-1,J)=CSV(J)*(V(I,J,K)*SCR(I,J,1) 1 -AM*(U1(I,J+1,K)-U1(I,J,K)))*IV(I,J,K) VY(I-1,J)=CSV(J)*(V(I,J,K)*SCR(I,J,2) 1 -AM*(V1(I,J+1,K)-V1(I,J,K)))*IV(I,J,K) SY(I-1,J)=CSV(J)*(V(I,J,K)*SCR(I,J,3) 1 -AH*(S1(I,J+1,K)-S1(I,J,K)))*IV(I,J,K) 380 TY(I-1,J)=CSV(J)*(V(I,J,K)*SCR(I,J,4) 1 -AH*(T1(I,J+1,K)-T1(I,J,K)))*IV(I,J,K) C ====================== C CONSERVATION EQUATIONS C ====================== C Note: These are only partial updates. Pressure and boundary forcings C are added after statement label 500. The below time difference fields C are interpolated to the "C" grid where their interpolated values are C added to the "old" time levels. This avoids interpolations of large C old time level fields and thus reduces numerical dispersion. DO 440 J=2,J1 TEMP=OCS(J)*ODY(J) ILO=ILEFT(J,K) IHI=IRITE(J,K) DO 440 I=ILO,IHI CTS DO 440 I=2,I1 DTIN=DT*IN(I,J,K) C Longitudinal momentum U2(I,J,K)=U1(I,J,K)-DTIN* 1 ((UX(I,J-1)-UX(I-1,J-1)+PX(I,J))*ODX(J) 2 +(UY(I-1,J)-UY(I-1,J-1))*TEMP 3 +(UZ(I-1,J-1,LT)-UZ(I-1,J-1,LB))*ODZ(K)) C Latitudinal momentum V2(I,J,K)=V1(I,J,K)-DTIN* 1 ((VX(I,J-1)-VX(I-1,J-1))*ODX(J) 1 +(VY(I-1,J)-VY(I-1,J-1))*TEMP+PY(I,J)*ODY(J) 2 +(VZ(I-1,J-1,LT)-VZ(I-1,J-1,LB))*ODZ(K)) C Salinity S2(I,J,K)=S1(I,J,K)-DTIN*((SX(I,J-1)-SX(I-1,J-1))*ODX(J) 1 +(SY(I-1,J)-SY(I-1,J-1))*TEMP 2 +(SZ(I-1,J-1,LT)-SZ(I-1,J-1,LB))*ODZ(K)) C Temperature 440 T2(I,J,K)=T1(I,J,K)-DTIN*((TX(I,J-1)-TX(I-1,J-1))*ODX(J) 1 +(TY(I-1,J)-TY(I-1,J-1))*TEMP 2 +(TZ(I-1,J-1,LT)-TZ(I-1,J-1,LB))*ODZ(K)) I=LB LB=LT 500 LT=I C C =============== C Surface effects C =============== C C Wind stress IF (LWIND.NE.0) CALL WIND(U2,V2,DT,ODZ) C Heat exchanges with atmosphere c IF (LHEAT.EQ.1) CALL QSURF(T2,DT,ODZ) C Radiative heating c IF (LSOL.EQ.1) CALL SUN(T2,DT,ODZ,QPROF,KB) IF (LOPEN.EQ.0) GO TO 510 CTS CTS ======================== CTS Open boundary conditions CTS ======================== CTS THE FOLLOWING IS REMOVED BY NESTING WITH CALCUR model CTS These are all determined by "known" normal boundary veloctiy (NBV) CTS i.e. boundary normal flux is UPWINDED for both inflow and outflow. CTS For inflow, exterior "ghost zone" fields are fluxed inward. The ghost CTS zone fields can be set to climatology. CTS Loop 506 below does SCALAR fluxes on boundaries CTS Loop 632 below does MOMENTUM fluxes CTS Loop 644 below determines NBV CTS UPWIND METHOD based on NBV. CTS Temporarily use ABS() function. CTS Later use boundary flag arrays if a vectorized function exists which CTS sets an integer array to the sign bit (i.e., 0 OR 1) CTS SKIP on first time step because advection velocity divergence in CTS boundary zones has not been cleared out yet. CTS cTS revised 8/5/99 to be nondamping, note that tangential momentum is cTS still not advected in ghostzone u2,v2 set after nbv update cTS before the pressure solver. TMP=.5*DT DO 506 K=1,K1 DO 504 J=2,J1 C West TMPIN=IN(2,J,K)*TMP*ODX(J) TEMP=ABS(U(1,J,K)) TEMP1=TMPIN*(TEMP+U(1,J,K)) TEMP2=TMPIN*(TEMP-U(1,J,K)) U2(2,J,K)=U2(2,J,K)-TEMP2*ULF(2,J,K)+TEMP1*U2(1,J,K) V2(2,J,K)=V2(2,J,K)-TEMP2*VLF(2,J,K)+TEMP1*V2(1,J,K) S2(2,J,K)=S2(2,J,K)-TEMP2*SLF(2,J,K)+TEMP1*S2(1,J,K) 504 T2(2,J,K)=T2(2,J,K)-TEMP2*TLF(2,J,K)+TEMP1*T2(1,J,K) C East boundary is land C Latitudinal boundaries DO 506 I=2,I1 C North TMPIN=IN(I,J1,K)*TMP*ODY(J1)*CSV(J1)*OCS(J1) TEMP=ABS(V(I,J1,K)) TEMP1=TMPIN*(TEMP+V(I,J1,K)) TEMP2=TMPIN*(TEMP-V(I,J1,K)) U2(I,J1,K)=U2(I,J1,K)-TEMP1*ULF(I,J1,K)+TEMP2*U2(I,J0,K) V2(I,J1,K)=V2(I,J1,K)-TEMP1*VLF(I,J1,K)+TEMP2*V2(I,J0,K) S2(I,J1,K)=S2(I,J1,K)-TEMP1*SLF(I,J1,K)+TEMP2*S2(I,J0,K) T2(I,J1,K)=T2(I,J1,K)-TEMP1*TLF(I,J1,K)+TEMP2*T2(I,J0,K) C South TMPIN=IN(I,2,K)*TMP*ODY(2)*CSV(1)*OCS(2) TEMP=ABS(V(I,1,K)) TEMP1=TMPIN*(TEMP+V(I,1,K)) TEMP2=TMPIN*(TEMP-V(I,1,K)) U2(I,2,K)=U2(I,2,K)+TEMP1*U2(I,1,K)-TEMP2*ULF(I,2,K) V2(I,2,K)=V2(I,2,K)+TEMP1*V2(I,1,K)-TEMP2*VLF(I,2,K) S2(I,2,K)=S2(I,2,K)+TEMP1*S2(I,1,K)-TEMP2*SLF(I,2,K) 506 T2(I,2,K)=T2(I,2,K)+TEMP1*T2(I,1,K)-TEMP2*TLF(I,2,K) C C ================================ C Climatological surface restoring C ================================ 510 TMP=1.-TAUDT C 180 day running time mean for surface heat and salinity fluxes TMP1=1./(180.*DAODT) TMP2=1.-TMP1 DO 580 J=2,J1 DO 580 I=2,I1 SAVE=S2(I,J,1) C First add dynamically adjusted running time average restoring rate S2(I,J,1)=SAVE+SALTOP(I,J) C Then nudge toward Levitus climatology S2(I,J,1)=TMP*S2(I,J,1)+TAUDT*SSURF(I-1,J-1) C Accumulate changes due to combined restoring rate C and nudging (slowly modifies long term restoring rate) SALTOP(I,J)=TMP2*SALTOP(I,J)+TMP1*(S2(I,J,1)-SAVE) SAVE=T2(I,J,1) C Add long term mean heating rate effect T2(I,J,1)=SAVE+HEATOP(I,J) C Accelerate surface temperature restoring for large difference between C climatological model and Levitus T T2(I,J,1)=TMP*T2(I,J,1)+TAUDT*TSURF(I-1,J-1) 580 HEATOP(I,J)=TMP2*HEATOP(I,J)+TMP1*(T2(I,J,1)-SAVE) C CTS REMOVE THE SPONGE LAYER C Sponge layer must be applied more slowly for stabililty !!! C ===================================== C Climatological sponge layer restoring C 10 boundary latitudes only C ===================================== CTS C ============= C Bottom stress C ============= TMP=DT*DRAG DO 600 J=2,J1 DO 600 I=2,I1 K=KB(I,J) DRG=IN(I,J,K)*TMP*ODZ(K)*SQRT(U1(I,J,K)**2+V1(I,J,K)**2) U2(I,J,K)=U2(I,J,K)-DRG*U1(I,J,K) 600 V2(I,J,K)=V2(I,J,K)-DRG*V1(I,J,K) C Trapezoidal Coriolis DO 602 K=1,K1 DO 602 J=2,J1 DO 602 I=2,I1 TMP=.5*(F(J-1)+ULF(I,J,K)*TANPHI(J-1))*DT TEMP=1./(1.+TMP**2) QU=U2(I,J,K)+TMP*V1(I,J,K) QV=V2(I,J,K)-TMP*U1(I,J,K) U2(I,J,K)=TEMP*(QU+TMP*QV) 602 V2(I,J,K)=TEMP*(QV-TMP*QU) C C ===================================================== C Interpolate "A" grid quantities to "C" grid locations C ===================================================== C Unnecessary masking has been removed here c 4th order interpolations do 1635 K=1,K1 DO 1625 J=2,J1 DO 1621 I=2,i2 1621 SCR(I,J,1)=6.*(U2(I,J,K)+U2(I+1,J,K)) C GO TO 623 C Skip loop 622 for 2nd order scheme ILO=MAX(3,ILO) IHI=MIN(I3,IHI) CTS do 1622 i=3,I3 do 1622 i=ilo,ihi 1622 SCR(I,J,1)=SCR(I,J,1) 1 -U2(I-1,J,K)+U2(I,J,K)+U2(I+1,J,K)-U2(I+2,J,K) c this allows for modification of coast effect above treats c land normal velocity as 0 rather than negative symmetric 1623 ILO=ILEFT(J,K) IHI=IRITE(J,K) DO 1625 I=2,i2 1625 U(I,J,K)=O12*SCR(I,J,1)*IU(I,J,K) DO 1631 J=2,J2 ILO=MAX(ILEFT(J,K),ILEFT(J+1,K)) IHI=MIN(IRITE(J,K),IRITE(J+1,K)) DO 1631 I=2,i1 1631 SCR(I,J,1)=6.*(V2(I,J,K)+V2(I,J+1,K)) C GO TO 633 C Skip loop 632 for 2nd order scheme DO 1632 J=3,J3 ILO=MAX(ILEFT(J,K),ILEFT(J+1,K)) IHI=MIN(IRITE(J,K),IRITE(J+1,K)) DO 1632 I=2,i1 1632 SCR(I,J,1)=SCR(I,J,1) 1 -V2(I,J-1,K)+V2(I,J,K)+V2(I,J+1,K)-V2(I,J+2,K) 1633 DO 1635 J=2,J2 ILO=MAX(ILEFT(J,K),ILEFT(J+1,K)) IHI=MIN(IRITE(J,K),IRITE(J+1,K)) DO 1635 I=2,i1 1635 V(I,J,K)=O12*SCR(I,J,1)*IV(I,J,K) C Nudge ghost zone climatology to accommodate westward and southward C surface Ekman layer drift (ignore Ekman effects for K > 1) C This ignores Ekman drift variation along boundary (LATER we may adjust C for that by keeping time average of normal flow interior to boundary) c goto 4048 TMP=0. DO 630 J=2,J1 630 TMP=TMP+(U(2,J,1)-U(1,J,1))/ODY(J) TMP=.01*TMP/(YV(J1)-YV(1)) c tmp is average change in U in y direction divided by 100 c so nudge velocities by adjusting mean and individual velocities UWSTM(1)=UWSTM(1)+TMP DO 631 J=2,J1 631 U(1,J,1)=U(1,J,1)+TMP TMP=0. TEMP=0. N=0 DO 632 I=2,I1 N=N+IN(I,2,1) TMP=TMP+IV(I,2,1)*V(I,2,1) C Yes, "IN" rather than "IV" 632 TEMP=TEMP+IN(I,2,1)*V(I,1,1) C Total flux (per unit depth) difference TMP=TMP*DXV(2)-TEMP*DXV(1) C Nudge velocity to decrease total flux difference TMP=.01*TMP/(N*DXV(1)) VSUDM(1)=VSUDM(1)+TMP DO 633 I=2,I1 633 V(I,1,1)=V(I,1,1)+IN(I,2,1)*TMP TMP=0. TEMP=0. N=0 DO 634 I=2,I1 N=N+IN(I,J1,1) TMP=TMP+IV(I,J2,1)*V(I,J2,1) 634 TEMP=TEMP+IN(I,J1,1)*V(I,J1,1) TMP=TMP*DXV(J2)-TEMP*DXV(J1) TMP=.01*TMP/(N*DXV(J1)) VNORM(1)=VNORM(1)+TMP DO 635 I=2,I1 635 V(I,J1,1)=V(I,J1,1)+IN(I,J1,1)*TMP IF(LOPEN.EQ.0) GO TO 651 C ================================================= C Specify time dependent open b.c.'s from CALCUR model C ================================================= C NBV update must be after loop 500 (to keep non-divergent bt mode) C Open boundary is totally determined by full CALCUR model DO 640 K=1,K1 DO 639 J=2,J1 C west U(1,J,K)=U(1,J,K)+UA2WST(J,K) U1(1,J,K)=U1(1,J,K)+U2WST(J,K) V1(1,J,K)=V1(1,J,K)+V2WST(J,K) S1(1,J,K)=S1(1,J,K)+S2WST(J,K) T1(1,J,K)=T1(1,J,K)+T2WST(J,K) U2(1,J,K)=U1(1,J,K) V2(1,J,K)=V1(1,J,K) S2(1,J,K)=S1(1,J,K) 639 T2(1,J,K)=T1(1,J,K) DO 640 I=2,I1 C south V(I,1,K)=V(I,1,K)+VA2SUD(I,K) U1(I,1,K)=U1(I,1,K)+U2SUD(I,K) V1(I,1,K)=V1(I,1,K)+V2SUD(I,K) S1(I,1,K)=S1(I,1,K)+S2SUD(I,K) T1(I,1,K)=T1(I,1,K)+T2SUD(I,K) U2(I,1,K)=U1(I,1,K) V2(I,1,K)=V1(I,1,K) S2(I,1,K)=S1(I,1,K) T2(I,1,K)=T1(I,1,K) C north V(I,J1,K)=V(I,J1,K)+VA2NOR(I,K) U1(I,J0,K)=U1(I,J0,K)+U2NOR(I,K) V1(I,J0,K)=V1(I,J0,K)+V2NOR(I,K) S1(I,J0,K)=S1(I,J0,K)+S2NOR(I,K) T1(I,J0,K)=T1(I,J0,K)+T2NOR(I,K) U2(I,J0,K)=U1(I,J0,K) V2(I,J0,K)=V1(I,J0,K) S2(I,J0,K)=S1(I,J0,K) 640 T2(I,J0,K)=T1(I,J0,K) C C ------------------------------ C Subtract mean boundary outflow as required by incompressibility C ------------------------------ C TMP is the volume flux exiting region C Check inflow calculation consistency TMP=0. DO 646 K=1,K1 DO 645 J=2,J1 645 TMP=TMP-U(1,J,K)*IN(2,J,K)/(ODY(J)*ODZ(K)) TMP1=1./(ODZ(K)*ODXV(1)) TMP2=1./(ODZ(K)*ODXV(J1)) DO 646 I=2,I1 646 TMP=TMP-IN(I,2,K)*TMP1*V(I,1,K) 1 +IN(I,J1,K)*TMP2*V(I,J1,K) c write(6,*) ' imbalance ',itf,tmp TMP=TMP/AROUT WRITE(6,647) TMP 647 FORMAT(1PE8.1,' cm/sec outflow vel adjust (+)=net div') C Adjust outflow to get zero net inflow c adjust on all three sides. DO 648 K=1,K1 DO 648 I=2,I1 V(I,1,K)=V(I,1,K)+TMP*IN(I,2,K) v(i,j1,k)=v(i,j1,k)-tmp*in(i,j1,k) c V1(I,1,K)=V(I,1,K) c v1(i,j0,k) = v(i,j1,k) c v2(i,j0,k) = v(i,j1,k) c V2(I,1,K)=V(I,1,K) c vlf(i,j0,k) = v(i,j1,k) c 648 VLF(I,1,K)=V(I,1,K) 648 continue do 1749 k=1,k1 do 1749 j=2,j1 u(1,j,k) = u(1,j,k)+tmp*in(2,j,k) c u1(1,j,k) = u(1,j,k) c u2(1,j,k) = u(1,j,k) c1749 ulf(1,j,k) = u(1,j,k) 1749 continue c TMP=0. c DO 1646 K=1,K1 c DO 1645 J=2,J1 c 1645 TMP=TMP-U(1,J,K)*IN(2,J,K)/(ODY(J)*ODZ(K)) c TMP1=1./(ODZ(K)*ODXV(1)) c TMP2=1./(ODZ(K)*ODXV(J1)) c DO 1646 I=2,I1 c 1646 TMP=TMP-IN(I,2,K)*TMP1*V(I,1,K) c 1 +IN(I,J1,K)*TMP2*V(I,J1,K) c TMP=TMP/AROUT c WRITE(6,*) ' after compressibbility correction ',TMP c c Iterate solver to get zero normal flow at shore c 651 MXMASK=3 C C Use SCR and SCR to store changes for incompressibility DO 652 K=1,4 DO 652 J=1,J0 DO 652 I=1,I0 652 SCR(I,J,K)=0. CALL SETR(UBT,I1*J0+I0*J1,0.) C C --------------------------------------------------- C Iterate EVP solver to get zero normal flow at shore C and exactly non-divergent bt mode C --------------------------------------------------- C MXMASK is maximum number of iterations allowed MXMASK=5 C C ================================ C Dan Wright's bt mode formulation C ================================ C C Calculate volume fluxes (depth weighted bt velocity) DO 656 K=1,K1 TMP=1./ODZ(K) DO 654 J=2,J1 UBT(1,J)=UBT(1,J)+TMP*U(1,J,K)*IN(2,J,K) UBT(I1,J)=UBT(I1,J)+TMP*U(I1,J,K)*IN(I1,J,K) DO 654 I=2,I2 654 UBT(I,J)=UBT(I,J)+TMP*U(I,J,K) DO 655 J=2,J2 DO 655 I=2,I1 655 VBT(I,J)=VBT(I,J)+TMP*V(I,J,K) DO 656 I=2,I1 VBT(I,1)=VBT(I,1)+TMP*V(I,1,K)*IN(I,2,K) 656 VBT(I,J1)=VBT(I,J1)+TMP*V(I,J1,K)*IN(I,J1,K) ITMASK=0 C 659 ITMASK=ITMASK+1 C C Zero out advection velocity over land TMP=0. DO 660 J=2,J1 DO 660 I=2,I2 U(I,J,1)=IU0(I,J)*U(I,J,1) 660 UBT(I,J)=IU0(I,J)*UBT(I,J) DO 661 J=2,J2 DO 661 I=2,I1 V(I,J,1)=IV0(I,J)*V(I,J,1) 661 VBT(I,J)=IV0(I,J)*VBT(I,J) DO 662 J=2,J1 TEMP=OCS(J)*ODY(J) DO 662 I=2,I1 662 S(I-1,J-1)=((UBT(I,J)-UBT(I-1,J))*ODX(J) 1 +(CSV(J)*VBT(I,J)-CSV(J-1)*VBT(I,J-1))*TEMP) C For rigid top get top level pressure adjustment from EVP solver and C barotropically adjust U, V, W CALL REP(AL,AB,AC,AR,AT,RINV,RINV1,DUM0,DUM1,DUM2,S,H,X,IE,I0,I2, 1 NB0) C Rigid-lid pressure adjustment DO 664 J=2,J1 DO 664 I=2,I1 664 P(I,J,1)=P(I,J,1)+ODT*X(I,J) C Barotropic mode adjustment (from pressure gradient change) DO 668 J=2,J1 DO 668 I=2,I2 668 SCR(I,J,1)=(X(I+1,J)-X(I,J))*ODX(J) DO 669 J=2,J2 DO 669 I=2,I1 669 SCR(I,J,2)=(X(I,J+1)-X(I,J))*ODYV(J) C C Accumulate top layer pressure gradient and adjust top layer bt mode C Remaining layers are handled by loop 684 below DO 674 J=2,J1 DO 674 I=2,I2 SCR(I,J,3)=SCR(I,J,3)-SCR(I,J,1) U(I,J,1)=U(I,J,1)-SCR(I,J,1) 674 UBT(I,J)=UBT(I,J)-ZBOTU(I,J)*SCR(I,J,1) DO 675 J=2,J2 DO 675 I=2,I1 SCR(I,J,4)=SCR(I,J,4)-SCR(I,J,2) V(I,J,1)=V(I,J,1)-SCR(I,J,2) 675 VBT(I,J)=VBT(I,J)-ZBOTV(I,J)*SCR(I,J,2) C Check convergence to zero flow over land TMP=0. DO 677 J=2,J1 DO 677 I=2,I2 677 TMP=MAX(TMP,(1.-IU0(I,J))*ABS(UBT(I,J))) DO 678 J=2,J2 DO 678 I=2,I1 678 TMP=MAX(TMP,(1.-IV0(I,J))*ABS(VBT(I,J))) TMP=TMP/Z(3) IF (ITF-IT0.LT.20.OR.MOD(ITF,36).EQ.0) WRITE(14,679) ITMASK,TMP IF (ITF-IT0.LT.99.OR.MOD(ITF,36).EQ.0) WRITE(*,679) ITMASK,TMP 679 FORMAT 1 (9X,'solver iteration #',I2,', max vel on land=',F9.5,' cm/sec') IF (TMP.GT..2.AND.ITMASK.LT.MXMASK.OR.ITMASK.LT.2) GO TO 659 c IF (TMP.GT..2.AND.ITMASK.LT.MXMASK) GO TO 659 C C Add bt mode changes to all layers except top layer (see loop 674) DO 684 K=2,K1 DO 683 J=2,J1 DO 683 I=2,I2 683 U(I,J,K)=U(I,J,K)+SCR(I,J,3)*IU(I,J,K) DO 684 J=2,J2 DO 684 I=2,I1 684 V(I,J,K)=V(I,J,K)+SCR(I,J,4)*IV(I,J,K) C This completes the advanced time level advection velocity C having exactly zero barotropic divergence and satisfying all C kinematic and specified open boundary conditions DO 685 K=1,K1 TMP=1./ODZ(K) DO 685 J=2,J1 TEMP=OCS(J)*ODY(J) DO 685 I=2,I1 685 W(I,J,K+1)=IW(I,J,K+1)*(W(I,J,K)-((U(I,J,K)-U(I-1,J,K))*ODX(J) 1 +(CSV(J)*V(I,J,K)-CSV(J-1)*V(I,J-1,K))*TEMP)*TMP) C C C Interpolate incompressibility induced changes back to cell centers C High order C NOTE: top layer has free-flow over land; C changes are independent of depth elsewhere DO 687 J=2,J1 DO 686 I=2,I1 686 SCR(I,J,1)=12.*(SCR(I-1,J,3)+SCR(I,J,3)) DO 687 I=3,I2 687 SCR(I,J,1)=SCR(I,J,1) 1 -SCR(I-2,J,3)+SCR(I-1,J,3)+SCR(I,J,3)-SCR(I+1,J,3) DO 689 K=1,K1 DO 689 J=2,J1 c new 8.5/99 ilo = max(3,ileft(j,k)) ihi = min(i2,irite(j,k)) do 688 i=ilo,ihi CTS do 688 i=2,I1 688 U2(I,J,K)=IN(I,J,K)*(U2(I,J,K)+O24*SCR(I,J,1)) u2(2,j,k) = 0.5*(u(1,j,k)+u(2,j,k)) u2(i1,j,k) = 0.5*(u(i2,j,k)+u(i1,j,k)) 689 CONTINUE DO 696 J=2,J1 DO 696 I=2,I1 696 SCR(I,J,1)=12.*(SCR(I,J-1,4)+SCR(I,J,4)) DO 697 J=3,J2 DO 697 I=2,I1 697 SCR(I,J,1)=SCR(I,J,1) 1 -SCR(I,J-2,4)+SCR(I,J-1,4)+SCR(I,J,4)-SCR(I,J+1,4) DO 700 K=1,K1 DO 698 J=2,J1 ilo = ileft(j,k) ihi = irite(j,k) DO 698 I=ilo,ihi CTS DO 698 I=2,J1 698 V2(I,J,K)=IN(I,J,K)*(V2(I,J,K)+O24*SCR(I,J,1)) c new 8/5/99 do 699 i=2,i1 v2(i,2,k) = 0.5*(v(i,1,k)+v(i,2,k)) v2(i,j1,k) = 0.5*(v(i,j2,k)+v(i,j1,k)) 699 continue 700 continue C IF (ITF-IT0.GT.50.AND.MOD(ITF,50).NE.0) GO TO 712 C ======================= C Check incompressibility C ======================= TMP=0. ERR=0. DO 710 K=1,K1 DO 710 J=2,J1 DO 710 I=2,I1 TMP1=(U(I,J,K)-U(I-1,J,K))*ODX(J) TMP2=(CSV(J)*V(I,J,K)-CSV(J-1)*V(I,J-1,K))*OCS(J)*ODY(J) TMP3=(W(I,J,K+1)-W(I,J,K))*ODZ(K) TMP=TMP+MAX(ABS(TMP1),ABS(TMP2),ABS(TMP3))*IN(I,J,K) 710 ERR=ERR+ABS(TMP1+TMP2+TMP3)*IN(I,J,K) ERR=ERR/TMP WRITE(*,711) ERR WRITE(14,711) ERR 711 FORMAT(' *** NORMALIZED mean incompressibility error = ',1PE9.2) 712 IF (MOD(ITF,M6).LT.3) WRITE(14,720) P(I0/2,J0/2,1),X(I0/2,J0/2) 720 FORMAT(' P,X ',2(1X,1PE10.3)) PSM=0. DO 724 J=2,J1 DO 724 I=2,I1 724 PSM=PSM+P(I,J,1) PSM=PSM/(I2*J2) DO 725 J=2,J1 DO 725 I=2,I1 725 P(I,J,1)=P(I,J,1)-PSM C C ======================== C Update using FLTW method C ======================== DO 745 K=1,K1 DO 745 J=2,J1 ilo = ileft(j,k) ihi = irite(j,k) DO 745 I=ilo,ihi S1(I,J,K)=OFLTW*SLF(I,J,K)+FLTW*(S1(I,J,K)+S2(I,J,K)) T1(I,J,K)=OFLTW*TLF(I,J,K)+FLTW*(T1(I,J,K)+T2(I,J,K)) U1(I,J,K)=OFLTW*ULF(I,J,K)+FLTW*(U1(I,J,K)+U2(I,J,K)) V1(I,J,K)=OFLTW*VLF(I,J,K)+FLTW*(V1(I,J,K)+V2(I,J,K)) SLF(I,J,K)=S2(I,J,K) TLF(I,J,K)=T2(I,J,K) ULF(I,J,K)=U2(I,J,K) 745 VLF(I,J,K)=V2(I,J,K) c new code 9/27/99 c Relax swamp layer and masked water to interior to reduce c seepage affects. note that swamp layer does advect and c use quantities over land do 777 j=2,j1 do 777 i=2,i1 t2(i,j,1) = in(i,j,1)*t2(i,j,1) + (1.-in(i,j,1))* > 0.25*(t2(i-1,j,1)+t2(i+1,j,1)+t2(i,j+1,1)+t2(i,j-1,1)) s2(i,j,1) = in(i,j,1)*s2(i,j,1) + (1.-in(i,j,1))* > 0.25*(s2(i-1,j,1)+s2(i+1,j,1)+s2(i,j+1,1)+s2(i,j-1,1)) 777 continue CTS REMOVED THE COMMENTS FOR INCREASING TOPLAYER MIXING C ========================== C Set vertical diffusivities C ========================== C Set background vertical diffusivities to DMZ0 and add numerical C contribution according to directed cell Reynolds number. c latest dietrich vertical diffusivities code c i modify initial wind induced surface forcing part by c using gaussian shape in vertical with 40 m depth scale c this reduces contribution to .006 at 400 meters c I wonder if i shouldn't just use some factor times tau instead c of a constant do 747 k=1,k2 tump = 20.*exp(-0.5*(z(2*k))/4000)*odzw(k+1) do 747 j=1,j2 do 747 i=1,i2 EV(i,j,k) = 0.0 747 HV(i,j,k) = 0.0 DEZ0 = 0.2*DMZ0 Do 760 k=1,k2 l = k+1 do 760 j=1,j2 do 760 i=1,i2 temp = orzmx*abs(w(i+1,j+1,l)) c pacanowski and philander(1981) evisc = a*(1/1+b*ri))**N c evisc = turbulent eddy viscosity diffuse=turbulent diffusivity ri=max(0.,980.*(rho(i,j,l)-rho(i,j,k))*odzw(l)/(odzw(l)**2* > (0.001+(u2(i,j,l)-u2(i,j,k))**2+(v2(i,j,l)-v2(i,j,k))**2))) tmp = 1./(1.+5.*ri) c reduce cross-pycnocline irreversible diffusion effects c kosher because internal waves dominate w at pycnocline, but c result in no net transports temp = temp*tmp evisc = 100.*tmp**2 evisc = min(evisc,10.) diffuse = (evisc*tmp+DEZ0)*odzw(l) ev(i,j,k) = (ev(i,j,k)+(evisc+dmz0)*odzw(l)+temp) > *iw(i+1,j+1,l) hv(i,j,k) = (hv(i,j,k)+diffuse+temp)*iw(i+1,j+1,l) 760 continue if(mod(itf,144).eq.0) then hmin = 999999. hmax = -999999. emin = 99999. emax = -99999. do k=1,k2 do j= 2,j2 do i=2,i2 hmin = min(hmin,hv(30,30,k)) hmax = max(hmax,hv(30,30,k)) emin = min(emin,ev(30,30,k)) emax = max(emax,ev(30,30,k)) enddo enddo enddo write(6,*) ' ** hminmax ',hmin,hmax write(6,*) ' ** eminmax ',emin,emax endif 793 continue CTS ================= CTS NO FILTER APPLIED CTS ================= C Damp cross-flow near all open boundaries (reduces reflections) c goto 877 c .999 originally N=J2/10 N= 10 DO 950 K=1,K1 C Damp LATITUDINAL velocity near open LONGITUDINAL boundaries C Skip closed eastern boundary C Western boundary TMP=1. N= 6 DO 930 I=N,2,-1 TMP=.9995*TMP DO 930 J=2,J1 930 V1(I,J,K)=TMP*V1(I,J,K) c goto 877 C Damp LONGITUDINAL velocity near open LATITUDINAL boundaries C Northern boundary TMP=1. N= 18 DO 940 J=J1-N,J1 TMP=.9995*TMP DO 940 I=2,I1 940 U1(I,J,K)=TMP*U1(I,J,K) C Southern boundary TMP=1. N = 8 DO 950 J=N,2,-1 TMP=.9995*TMP DO 950 I=2,I1 950 U1(I,J,K)=TMP*U1(I,J,K) 877 continue C C Add turbulence closure scheme contributions to vertical diffusivities. C IF (LTURB.EQ.1) CALL DIFUSE(HV,EV,UC,VC,RHO,ODZW,SCR) C CALL DIFUSE(HV,EV,U2,V2,RHO,ODZW,SCR) END C ---------------------------------------------------------------------- SUBROUTINE DIFUSE(HV,EV,U,V,RHO,ODZW,SCR) C ---------------------------------------------------------------------- INCLUDE 'resolution.h' PARAMETER(I2=I0-2,J1=J0-1,J2=J0-2,K1=K0-1,K2=K0-2) IMPLICIT REAL*8 (A-H,O-Z) INTEGER*2 KB,IU0,IV0,IN,IU,IV,IW COMMON/WINDS/TAUX(I2,J2),TAUY(I2,J2) DIMENSION EV(I2,J2,K2),HV(I2,J2,K2),RHO(I0,J0,K1), 1 U(I0,J0,K1),V(I0,J0,K1),ODZW(K0),SCR(I0,J0,K1) COMMON/TOPOG/KB(I0,J0),IU0(I0,J0),IV0(I0,J0), 1 IN(I0,J0,K1),IU(I0,J0,K1),IV(I0,J1,K1),IW(I0,J0,K0) C C This subroutine adds turbulence contributions to diffusivities HV,EV. C This is not a transportive turbulence closure scheme. C The velocity arrays U,V and temperature, T, have perimeter "ghost zones". C There are no vertical "ghost zones"; K1=total number of layers. C The diffusivity arrays HV,EV have no "ghost zones". C HV and EV are stored at the vertical velocity (layer interface) locations. C They are not defined or stored at the top surface or bottom interface of the C deepest layer, where diffusive flux effects are specified by other means. C C Add any desired values HVAD,EVAD to HV and EV here. C Take SQRT(wind stress) once-and-for-all at a give horizontal location DO 100 J=1,J2 DO 100 I=1,I2 100 SCR(I,J,1)=5.*SQRT(TAUX(I,J)*TAUX(I,J)+TAUY(I,J)*TAUY(I,J)) C Add to HV only for top five layers DO 200 K=1,5 DO 200 J=1,J2 DO 200 I=1,I2 DVDZ=ODZW(K)*SQRT((U(I+1,J+1,K+1)-U(I+1,J+1,K))**2 1 +(V(I+1,J+1,K+1)-V(I+1,J+1,K))**2) C DZ increases downward in DRHODZ= d(rh0)/dz DRHODZ=(RHO(I+1,J+1,K+1)-RHO(I+1,J+1,K))*ODZW(K) C Make RI=0 for negative stratification DRHODZ=MAX(DRHODZ,0.) RI=980.*DRHODZ/(0.0001+DVDZ**2) EVAD=SCR(I,J,1)/((1.0+10.0*RI)**1.5)*ODZW(K+1)*IW(I+1,J+1,K+1) HVAD=0.5*EVAD c--- increase temp mixing for cold water over warm c if(t(i,j,k).lt.t(i,j,k+1)) hvad=10.0 c--- HV(I,J,K)=HV(I,J,K)+HVAD 200 EV(I,J,K)=EV(I,J,K)+EVAD END C ---------------------------------------------------------------------- SUBROUTINE WIND(U,V,DT,ODZ) C ---------------------------------------------------------------------- INCLUDE 'resolution.h' PARAMETER(I1=I0-1,J1=J0-1,I2=I0-2,J2=J0-2) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION U(I0,J0),V(I0,J0),ODZ(*) COMMON/WINDS/TAUX(I2,J2),TAUY(I2,J2) C C The velocity arrays U,V have a perimeter "ghost zone". C Thus, we set wind forcing only in the interior zones. C TAUX,TAUY are surface wind stress components. C TAUX,TAUY units are force per unit area (i.e., energy per unit volume). C C TMP=ODZ(1)/RHO C All units are cgs, so we use RHO=1. TMP=DT*ODZ(1) DO 100 J=2,J1 DO 100 I=2,I1 U(I,J)=U(I,J)+TMP*TAUX(I-1,J-1) 100 V(I,J)=V(I,J)+TMP*TAUY(I-1,J-1) END C ---------------------------------------------------------------------- SUBROUTINE QSURF(T,DT,ODZ) C ---------------------------------------------------------------------- INCLUDE 'resolution.h' PARAMETER(I2=I0-2,J2=J0-2) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION T(I0,J0),ODZ(*) C C The temperature array T has a permiter "ghost zone". C Thus, we set atmospheric heat exchange only in the interior zones. C C Typical latitudinal range: -50 to 50 Watts per square meter. C 1 Watt per square meter = 1 erg per second per equare cm. C C ------- C EXAMPLE: surface heating QDOT = 50 Watts per square meter C ------- QDOT=50. C TMP is conversion factor from Watts per square meter C to deg C per time step in top model layer. C TMP=DT*ODZ(1)/RHO/CP C All units are cgs, so we use rho=1. C CP for water is 2.5E4 ergs/gram/deg C. CHECK THIS CP VALUE!!! TMP=DT*ODZ(1)/2.5E4 DTEMP=TMP*QDOT DO 100 J=1,J2 DO 100 I=1,I2 100 T(I+1,J+1)=T(I+1,J+1)+DTEMP C ------- C EXAMPLE C ------- C END C ---------------------------------------------------------------------- C ****************** C GRAPHICS PACKAGE * C ****************** C ---------------------------------------------------------------------- SUBROUTINE FSPLTS(MXIT,NFLG) C SAVES PLOT DATA C ---------------------------------------------------------------------- INCLUDE 'resolution.h' PARAMETER(IJ0=I0+J0,I1=I0-1,I2=I0-2,J1=J0-1,J2=J0-2, 1 K1=K0-1,K2=K0-2) IMPLICIT REAL*8 (A-H,O-Z) INTEGER*2 KB,IU0,IV0,IN,IU,IV,IW COMMON/AGRID/ RHO(I0,J0,K1),PX(I0,J0),PY(I0,J0), 1 U1(I0,J0,K1),U2(I0,J0,K1),V1(I0,J0,K1),V2(I0,J0,K1), 2 S1(I0,J0,K1),S2(I0,J0,K1),T1(I0,J0,K1),T2(I0,J0,K1),P(I0,J0,K1), 3 ULF(I0,J0,K1),VLF(I0,J0,K1),SLF(I0,J0,K1),TLF(I0,J0,K1), 4 EV(I2,J2,K2),HV(I2,J2,K2),DAVG(K1),F(J2),TANPHI(J2),SUMIN(K1) COMMON/SURF/SSURF(I2,J2),TSURF(I2,J2) COMMON/CGRID/U(I0,J0,K1),V(I0,J1,K1),W(I0,J0,K0) COMMON/KLIMAT/SCR(I0,J0,K1),SC1(IJ0),SC2(IJ0),TW(I0,J0,2), 1 TCLI(I0,J0,K0),UCLI(I0,J0),VCLI(I0,J0),RMSV(I0,J0), 2 XCLI(I0,J0),PCLI(I0,J0,2),HEATOP(I0,J0),SALTOP(I0,J0) COMMON/SCA/ODT,PRN,DAYS,ORXYMX,ORZMX,OFLTW,OM5,OM8,ITF,IT0 COMMON/METRICS/ 1 Y(J0),YV(J0),YVDEG(J0),YDEG(J0), 2 CS(J0),CSV(J1),OCS(J0),OCSV(J1), 3 DX(J0),DXV(J1),ODX(J0),ODXV(J1), 4 DY(J0),DYV(J1),ODY(J0),ODYV(J1) CTS NESTING BOUNDARY COMMON/BCS/ 1 U2WST(J0,K1),V2WST(J0,K1),S2WST(J0,K1),T2WST(J0,K1), 3 U2SUD(I0,K1),V2SUD(I0,K1),S2SUD(I0,K1),T2SUD(I0,K1), 4 U2NOR(I0,K1),V2NOR(I0,K1),S2NOR(I0,K1),T2NOR(I0,K1), 5 UA2WST(J0,K1),VA2SUD(I0,K1),VA2NOR(I0,K1) COMMON/TOPOG/KB(I0,J0),IU0(I0,J0),IV0(I0,J0), 1 IN(I0,J0,K1),IU(I0,J0,K1),IV(I0,J1,K1),IW(I0,J0,K0) COMMON/ZFS/Z(K0+K1),ODZ(K1),ODZW(K0) C C --------- C X-Y PLOTS C --------- KSTOP=K1 IF (NFLG.NE.0.AND.ITF.NE.MXIT) KSTOP=1 K=1 10 DO 15 J=2,J1 DO 15 I=2,I1 C Scale up density to kg/m3 SCR(I,J,5)=1000.*RHO(I,J,K) C Convert to equivalent free surface anomaly SCR(I,J,1)=P(I,J,K)/980. 15 SCR(I,J,4)=T2(I,J,K)-TCLI(I,J,K) CALL XYPLOT('P','pressure-equivalent f.s.a.(cm)',SCR,K,K) CALL XYPLOT('u','longitudinal velocity (cm/sec)',U2(1,1,K),K,K) CALL XYPLOT('v','latitudinal velocity (cm/sec) ',V2(1,1,K),K,K) C CALL XYPLOT('R','rho - (hor. avg. rho) ',SCR(1,1,5),K,K) CALL XYPLOT('S','salinity (parts per thousand) ',S2(1,1,K),K,K) CALL XYPLOT('T','temperature (degrees c) ',T2(1,1,K),K,K) c CALL XYPLOT('t','T-(time average T) (degrees C)',SCR(1,1,4),K,K) C IF (ITF.NE.MXIT) GO TO 30 DO 25 J=2,J1 DO 25 I=2,I1 SCR(I,J,1)=TMP*(W(I,J,K)+W(I,J,K+1)) IF (I.NE.2.AND.J.NE.2.AND.I.NE.I1.AND.J.NE.J1) GO TO 24 N=IN(I-1,J-1,K)*IN(I,J-1,K)*IN(I+1,J-1,K) 1 *IN(I-1, J,K)*IN(I, J,K)*IN(I+1, J,K) 2 *IN(I-1,J+1,K)*IN(I,J+1,K)*IN(I+1,J+1,K) GO TO 25 24 N=IN(I-2,J-2,K)*IN(I-1,J-2,K)*IN(I,J-2,K)*IN(I+1,J-2,K) 1 *IN(I+2,J-2,K)*IN(I-2,J-1,K)*IN(I-1,J-1,K)*IN(I,J-1,K) 2 *IN(I+1,J-1,K)*IN(I+2,J-1,K)*IN(I-2,J,K)*IN(I-1,J,K) 3 *IN(I,J,K)*IN(I+1,J,K)*IN(I+2,J,K)*IN(I-2,J+1,K)*IN(I-1,J+1,K) 4 *IN(I,J+1,K)*IN(I+1,J+1,K)*IN(I+2,J+1,K)*IN(I-2,J+2,K) 5 *IN(I-1,J+2,K)*IN(I,J+2,K)*IN(I+1,J+2,K)*IN(I+2,J+2,K) 25 SCR(I,J,2)=N*SCR(I,J,1) c CALL XYPLOT('w','vertical velocity (meters/day)',SCR,K,K) CALL XYPLOT('W','interior vertical vel. (m/day)',SCR(1,1,2),K,K) 30 IF (ITF.NE.MXIT) K=MIN(2*K,K1) IF (ITF.EQ.MXIT) K=K+1 IF (K.LT.KSTOP) GO TO 10 C IF (KSTOP.EQ.1) GO TO 100 C IF (KSTOP.NE.K1) RETURN IF (KSTOP.NE.K1.OR.ITF.NE.MXIT) RETURN C --------------------- C CLIMATOLOGICAL FIELDS C --------------------- DO 35 J=2,J1 DO 35 I=2,I1 SCR(I,J,1)=PCLI(I,J,1)/980. SCR(I,J,2)=SQRT(PCLI(I,J,2))/980. SCR(I,J,3)=SQRT(U2(I,J,1)**2+V2(I,J,1)**2) 35 SCR(I,J,4)=HV(I-1,J-1,1)/ODZW(2) C35 SCR(I,J,5)=EV(I-1,J-1,1)/ODZW(2) CALL XYPLOT('x','time average P (eq. fsa, cm) ',SCR,1,1) CALL XYPLOT('e','eddy rms surface elevation(cm)',SCR(1,1,2),1,1) CALL XYPLOT('m','horizontal velocity magnitude ',SCR(1,1,3),1,1) CALL XYPLOT('m','longitudinal velocity (cm/sec)',U2,1,1) CALL XYPLOT('m','latititudinal velocity(cm/sec)',V2,1,1) C Show unit interval isotherms for top layer only (to compare with Levitus) CALL XYPLOT('T','time average T (degrees C) ',TCLI,1,1) CALL XYPLOT('C','time average T (degrees C) ',TCLI,1,K1) C CALL XYPLOT('z','z-heat diffusivity (cm-cm/sec)',SCR(1,1,4),1,1) C CALL XYPLOT('z','z-mom. diffusivity (cm-cm/sec)',SCR(1,1,5),1,1) C --------- C X-Z PLOTS C --------- 100 DO 130 J=2,J1,12 TMP=7.*24.*36. TEMP=25.*TMP DO 110 K=1,K1 DO 110 I=2,I1 SCR(I,1,K)=P(I,J,K)/980. SCR(I,2,K)=T2(I,J,K)-TCLI(I,J,K) SCR(I,4,K)=.5*TMP*(W(I,J,K)+W(I,J,K+1)) 110 SCR(I,5,K)=TEMP*( 1 (V(I,J-1,K)-V(I-1,J-1,K)+V(I+1,J-1,K)-V(I,J-1,K))*ODXV(J-1) 2 +(V(I+1,J,K)-V(I,J,K)+V2(I,J,K)-V2(I-1,J,K))*ODXV(J) 3 -(U(I-1,J,K)-U(I-1,J-1,K)+U(I,J,K)-U(I,J-1,K))*ODYV(J-1) 4 -(U(I,J+1,K)-U(I,J,K)+U2(I-1,J+1,K)-U2(I-1,J,K))*ODYV(J)) CALL XZPLOT('P','pressure-equivalent f.s.a.(cm)',SCR,J,J) CALL XZPLOT('S','salinity (parts per thousand) ',S2(1,J,1),J,J) CALL XZPLOT('T','temperature (degrees C) ',T2(1,J,1),J,J) CALL XZPLOT('w','vertical velocity(meters/week)',SCR(1,4,1),J,J) CALL XZPLOT('s','vorticity (per week) ',SCR(1,5,1),J,J) 130 CALL XZPLOT('t','T-(time average T) (degrees C)',SCR(1,2,1),J,J) C --------- C Y-Z PLOTS C --------- C Hardwire for 1/2 deg. North Atlantic deep water and GIN Sea DO 230 I=2,I1,12 DO 210 K=1,K1 DO 210 J=2,J1 C Scale up density SCR(3,J,K)=1000.*RHO(I,J,K) SCR(1,J,K)=P(I,J,K)/980. 210 SCR(2,J,K)=TMP*.5*(W(I,J,K)+W(I,J,K+1)) CALL YZPLOT('P','pressure-equivalent f.s.a.(cm)',SCR,I,I) CALL YZPLOT('R','(RHO-RHOmin) (kg/m3) ',SCR(3,1,1),I,I) CALL YZPLOT('S','salinity (parts per thousand) ',S2(I,1,1),I,I) CALL YZPLOT('T','temperature (degrees c) ',T2(I,1,1),I,I) 230 CALL YZPLOT('w','vertical velocity(meters/week)',SCR(2,1,1),I,I) END C ---------------------------------------------------------------------- SUBROUTINE COMXY(FLD,DAYS,IN) C SAVES MOVIE DATA C ---------------------------------------------------------------------- INCLUDE 'resolution.h' PARAMETER(IJ0=I0+J0,I1=I0-1,J1=J0-1,K1=K0-1) IMPLICIT REAL*8 (A-H,O-Z) INTEGER*2 NFLD,IN COMMON/KLIMAT/SCR(I0,J0,K1),SC1(IJ0),SC2(IJ0),TW(I0,J0,2), 1 TCLI(I0,J0,K0),UCLI(I0,J0),VCLI(I0,J0),RMSV(I0,J0), 2 XCLI(I0,J0),PCLI(I0,J0,2),HEATOP(I0,J0),SALTOP(I0,J0) COMMON/FLOAT/NFLD(I0,J0) DIMENSION IN(I0,1),FLD(I0,1) CALL RANGE(FLD,SC1,SC2,IN,I0,2,2,I1,J1,ILO,JLO,IHI,JHI, 1 FMIN,FMAX) C WRITE(*,5) ILO,IHI,JLO,JHI,FMIN,FMAX C5 FORMAT('ILO,IHI,JLO,JHI,FMIN,FMAX',4I4,2(1X,1PE8.1)) FRNG=FMAX-FMIN RMIN=FMIN+1.E-5*FRNG RMAX=FMAX-1.E-5*FRNG IF ((RMAX-RMIN).EQ.0.) RETURN 2 RF=9999./FRNG C NFLD=9999.*(FLD-RMIN)/FRNG DO 11 J=2,J1 DO 11 I=2,I1 11 NFLD(I,J)=MIN(9999.,RF*(FLD(I,J)-RMIN)) WRITE(21) 1 DAYS,ILO,JLO,IHI,JHI,RMIN,RMAX,((NFLD(I,J),I=2,I1),J=2,J1) END C ---------------------------------------------------------------------- SUBROUTINE XYPLOT(FN,FNAME,FLD,KL,KH) C ---------------------------------------------------------------------- INCLUDE 'resolution.h' PARAMETER(IJ0=I0+J0,I1=I0-1,J1=J0-1,K1=K0-1) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER FN*1,FNAME*30 INTEGER*2 KB,IU0,IV0,IN,IU,IV,IW,NFLD COMMON/KLIMAT/SCR(I0,J0,K1),SC1(IJ0),SC2(IJ0),TW(I0,J0,2), 1 TCLI(I0,J0,K0),UCLI(I0,J0),VCLI(I0,J0),RMSV(I0,J0), 2 XCLI(I0,J0),PCLI(I0,J0,2),HEATOP(I0,J0),SALTOP(I0,J0) COMMON/FLOAT/NFLD(I0,J0) COMMON/TOPOG/KB(I0,J0),IU0(I0,J0),IV0(I0,J0), 1 IN(I0,J0,K1),IU(I0,J0,K1),IV(I0,J1,K1),IW(I0,J0,K0) COMMON/SCA/ODT,PRN,DAYS,ORXYMX,ORZMX,OFLTW,OM5,OM8,ITF,IT0 DIMENSION FLD(I0,J0,1) DO 100 K=KL,KH KK=K-KL+1 CALL RANGE(FLD(1,1,KK),SC1,SC2,IN(1,1,K),I0,2,2,I1,J1, 1 ILO,JLO,IHI,JHI,FMIN,FMAX) C WRITE(*,5) FN,K,FMIN,ILO,JLO,FMAX,IHI,JHI C IF (ITF.LT.100) WRITE(14,5) FN,K,FMIN,ILO,JLO,FMAX,IHI,JHI 5 FORMAT(A1,' AT LEVEL ',I2,': MIN=',1PE9.2,' AT (',I3,',',I3, 1 '); MAX=',1PE9.2,' AT (',I3,',',I3,')') FRNG=FMAX-FMIN RMIN=FMIN+1.E-5*FRNG RMAX=FMAX-1.E-5*FRNG IF ((RMAX-RMIN).EQ.0.) RETURN 2 RF=9999./FRNG C NFLD=9999.*(FLD-RMIN)/FRNG DO 11 J=2,J1 DO 11 I=2,I1 11 NFLD(I,J)=MIN(9999.,RF*(FLD(I,J,KK)-RMIN)) IF (DAYS.NE.0.) GO TO 99 WRITE(28) FN,FNAME,DAYS,K,ILO,JLO,IHI,JHI,RMIN,RMAX, 1 ((NFLD(I,J),I=2,I1),J=2,J1) GO TO 100 99 IF (K.EQ.1.OR.(K.EQ.4.AND.(FN.EQ.'W'.OR.FN.EQ.'w'))) 1 WRITE(30) FN,FNAME,DAYS,K,ILO,JLO,IHI,JHI,RMIN,RMAX, 2 ((NFLD(I,J),I=2,I1),J=2,J1) IF (K.NE.1.OR.FN.EQ.'X') 1 WRITE(29) FN,FNAME,DAYS,K,ILO,JLO,IHI,JHI,RMIN,RMAX, 2 ((NFLD(I,J),I=2,I1),J=2,J1) 100 CONTINUE END C ---------------------------------------------------------------------- SUBROUTINE XZPLOT(FN,FNAME,FLD,JL,JH) C ---------------------------------------------------------------------- INCLUDE 'resolution.h' PARAMETER(IJ0=I0+J0,I1=I0-1,J1=J0-1,K1=K0-1) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER FN*1,FNAME*30 INTEGER*2 KB,IU0,IV0,IN,INXZ,IU,IV,IW,NFLD COMMON/KLIMAT/SCR(I0,J0,K1),SC1(IJ0),SC2(IJ0),TW(I0,J0,2), 1 TCLI(I0,J0,K0),UCLI(I0,J0),VCLI(I0,J0),RMSV(I0,J0), 2 XCLI(I0,J0),PCLI(I0,J0,2),HEATOP(I0,J0),SALTOP(I0,J0) COMMON/TOPOG/KB(I0,J0),IU0(I0,J0),IV0(I0,J0), 1 IN(I0,J0,K1),IU(I0,J0,K1),IV(I0,J1,K1),IW(I0,J0,K0) COMMON/SCA/ODT,PRN,DAYS,ORXYMX,ORZMX,OFLTW,OM5,OM8,ITF,IT0 DIMENSION FLD(I0,J0,*),NFLD(I0,K1),XZ(I0,K1),INXZ(I0,K1) DO 100 J=JL,JH JJ=J-JL+1 DO 10 K=1,K1 DO 10 I=2,I1 XZ(I,K)=FLD(I,JJ,K) 10 INXZ(I,K)=IN(I,J,K) CALL RANGE(XZ,SC1,SC2,INXZ,I0,2,1,I1,K1,ILO,KLO,IHI,KHI, 1 FMIN,FMAX) C IF (ITF.LT.100) WRITE(14,5) FN,J,FMIN,FMAX 5 FORMAT(A1,': MIN,MAX AT LATITUDE ',I2,' = ',2(1X,1PE9.2)) FRNG=FMAX-FMIN RMIN=FMIN+1.E-5*FRNG RMAX=FMAX-1.E-5*FRNG IF ((RMAX-RMIN).EQ.0.) RETURN 2 RF=9999./FRNG C NFLD=9999.*(FLD-RMIN)/FRNG DO 11 K=1,K1 DO 11 I=2,I1 11 NFLD(I,K)=MIN(9999.,RF*(XZ(I,K)-RMIN)*INXZ(I,K)) 100 WRITE(31) FN,FNAME,DAYS,J,ILO,KLO,IHI,KHI,RMIN,RMAX, 1 ((NFLD(I,K),I=2,I1),K=1,K1) END C ---------------------------------------------------------------------- SUBROUTINE YZPLOT(FN,FNAME,FLD,IL,IH) C ---------------------------------------------------------------------- INCLUDE 'resolution.h' PARAMETER(IJ0=I0+J0,I1=I0-1,J1=J0-1,K1=K0-1) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER FN*1,FNAME*30 INTEGER*2 KB,IU0,IV0,IN,INYZ,IU,IV,IW,NFLD COMMON/KLIMAT/SCR(I0,J0,K1),SC1(IJ0),SC2(IJ0),TW(I0,J0,2), 1 TCLI(I0,J0,K0),UCLI(I0,J0),VCLI(I0,J0),RMSV(I0,J0), 2 XCLI(I0,J0),PCLI(I0,J0,2),HEATOP(I0,J0),SALTOP(I0,J0) COMMON/TOPOG/KB(I0,J0),IU0(I0,J0),IV0(I0,J0), 1 IN(I0,J0,K1),IU(I0,J0,K1),IV(I0,J1,K1),IW(I0,J0,K0) COMMON/SCA/ODT,PRN,DAYS,ORXYMX,ORZMX,OFLTW,OM5,OM8,ITF,IT0 DIMENSION FLD(I0,J0,*),NFLD(J0,K1),YZ(J0,K1),INYZ(J0,K1) DO 100 I=IL,IH II=I-IL+1 DO 10 K=1,K1 DO 10 J=2,J1 YZ(J,K)=FLD(II,J,K) 10 INYZ(J,K)=IN(I,J,K) C C For plotting zonal averages (indicated by I=1) IF (I.NE.1) GO TO 20 DO 15 K=1,K1 DO 15 J=2,J1 N=0 DO 14 III=2,I1 14 N=N+IN(III,J,K) INYZ(J,K)=0 15 INYZ(J,K)=MIN(1,N) C 20 CALL RANGE(YZ,SC1,SC2,INYZ,J0,2,1,J1,K1,JLO,KLO,JHI,KHI, 1 FMIN,FMAX) C IF (ITF.LT.100) WRITE(14,21) FN,I,FMIN,FMAX 21 FORMAT(A1,': MIN,MAX AT LONGITUDE ',I3,' = ',2(1X,1PE9.2)) FRNG=FMAX-FMIN RMIN=FMIN+1.E-5*FRNG RMAX=FMAX-1.E-5*FRNG IF ((RMAX-RMIN).EQ.0.) RETURN RF=9999./FRNG C NFLD=9999.*(FLD-RMIN)/FRNG DO 30 K=1,K1 DO 30 J=2,J1 30 NFLD(J,K)=MIN(9999.,RF*(YZ(J,K)-RMIN)*INYZ(J,K)) 100 WRITE(32) FN,FNAME,DAYS,I,JLO,KLO,JHI,KHI,RMIN,RMAX, 1 ((NFLD(J,K),J=2,J1),K=1,K1) END C ---------------------------------------------------------------------- SUBROUTINE RANGE 1 (FLD,SEGMN,SEGMX,IN,IR,IL,JL,IH,JH,ILO,JLO,IHI,JHI,FMIN,FMAX) C VECTORIZED SEARCH FOR MAX'S AND MIN'S C ---------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) INTEGER*2 IN DIMENSION FLD(IR,*),SEGMN(*),SEGMX(*),IN(IR,*) LBLK=1 LAB=1 IHM=IH-1 C ---------------- C SEARCH FOR MAXES C ---------------- C C Initialize column maxes DO 2 I=IL,IH 2 SEGMX(I)=-1.E30 IF (LBLK.EQ.0) GO TO 9 C C Search for column maxes C with islands DO 5 J=JL,JH DO 5 I=IL,IH 5 SEGMX(I)=(1.-IN(I,J))*SEGMX(I)+IN(I,J)*MAX(SEGMX(I),FLD(I,J)) GO TO 15 C without islands 9 DO 10 J=JL,JH DO 10 I=IL,IH 10 SEGMX(I)=MAX(SEGMX(I),FLD(I,J)) C C Search for max of column maxes 15 FMAX=SEGMX(IL) DO 20 I=IL,IHM 20 FMAX=MAX(FMAX,SEGMX(I+1)) IF (LAB.EQ.0) GO TO 50 C C Search for location of max DO 29 I=IL,IH 29 IF (SEGMX(I).EQ.FMAX) IHI=I DO 30 J=JL,JH 30 IF (FLD(IHI,J).EQ.FMAX) JHI=J C --------------- C SEARCH FOR MINS C --------------- 50 DO 52 I=IL,IH 52 SEGMN(I)=1.E30 IF (LBLK.EQ.0) GO TO 59 DO 55 J=JL,JH DO 55 I=IL,IH 55 SEGMN(I)=(1.-IN(I,J))*SEGMN(I)+IN(I,J)*MIN(SEGMN(I),FLD(I,J)) GO TO 65 59 DO 60 J=JL,JH DO 60 I=IL,IH 60 SEGMN(I)=MIN(SEGMN(I),FLD(I,J)) 65 FMIN=SEGMN(IL) DO 70 I=IL,IHM 70 FMIN=MIN(FMIN,SEGMN(I+1)) IF (LAB.EQ.0) RETURN DO 79 I=IL,IH 79 IF (SEGMN(I).EQ.FMIN) ILO=I DO 80 J=JL,JH 80 IF (FLD(ILO,J).EQ.FMIN) JLO=J END C ---------------------------------------------------------------------- SUBROUTINE SAV3D(FN,FNAME,FLD,KL,KH,DAYS) C ---------------------------------------------------------------------- INCLUDE 'resolution.h' PARAMETER(IJ0=I0+J0,I1=I0-1,J1=J0-1,K1=K0-1) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER FN*1,FNAME*30 INTEGER*2 KB,IU0,IV0,IN,IU,IV,IW,NFLD COMMON/KLIMAT/SCR(I0,J0,K1),SC1(IJ0),SC2(IJ0),TW(I0,J0,2), 1 TCLI(I0,J0,K0),UCLI(I0,J0),VCLI(I0,J0),RMSV(I0,J0), 2 XCLI(I0,J0),PCLI(I0,J0,2),HEATOP(I0,J0),SALTOP(I0,J0) COMMON/FLOAT3/NFLD(I0,J0,K0) COMMON/TOPOG/KB(I0,J0),IU0(I0,J0),IV0(I0,J0), 1 IN(I0,J0,K1),IU(I0,J0,K1),IV(I0,J1,K1),IW(I0,J0,K0) DIMENSION FLD(I0,J0,1) RMAX=-1.E20 RMIN=1.E20 DO 9 K=KL,KH CALL RANGE(FLD(1,1,K),SC1,SC2,IN(1,1,K),I0,2,2,I1,J1,IL,JL,IH, 1 JH,FMIN,FMAX) RMAX=MAX(RMAX,FMAX) RMIN=MIN(RMIN,FMIN) IF (RMAX.NE.FMAX) GO TO 5 IHI=IH JHI=JH KHI=K GO TO 9 5 IF (RMIN.NE.FMIN) GO TO 9 ILO=IL JLO=JL KLO=K 9 CONTINUE FRNG=RMAX-RMIN RMIN=RMIN+1.E-5*FRNG RMAX=RMAX-1.E-5*FRNG IF ((RMAX-RMIN).EQ.0.) RETURN RF=9999./FRNG C NFLD=9999.*(FLD-RMIN)/FRNG DO 11 K=KL,KH DO 11 J=2,J1 DO 11 I=2,I1 11 NFLD(I,J,K)=MIN(9999.,RF*(FLD(I,J,K)-RMIN)) DO 12 K=2,K1 DO 12 J=2,J1 DO 12 I=2,I1 12 NFLD(I,J,K)=NFLD(I,J,K)*IN(I,J,K) C This binary output was used for original SAV3D.4YR C WRITE(40) FN,FNAME,DAYS,KL,KH,ILO,JLO,KLO,IHI,JHI,KHI,RMIN,RMAX, C 1 (((NFLD(I,J,K),I=2,I1),J=2,J1),K=KL,KH) C WRITE(40,20) FN,FNAME,DAYS,KL,KH,ILO,JLO,KLO,IHI,JHI,KHI,RMIN, C 1 RMAX,(((NFLD(I,J,K),I=2,I1),J=2,J1),K=KL,KH) C20 FORMAT(A1,A30,E10.3,8I2,2E12.5/(20I4)) END C ---------------------------------------------------------------------- SUBROUTINE SETR(A,N,R) C ---------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) real r DIMENSION A(1) DO 10 I=1,N 10 A(I)=R END C ***************** C ELLIPTIC SOLVER * C ***************** C ---------------------------------------------------------------------- SUBROUTINE REP(AX,AY,BB,CX,CY,RINV,RINV1,DUM0,DUM1,DUM2,F,H,X,IE, 1 I0,I2,NBLK) IMPLICIT REAL*8 (A-H,O-Z) C ---------------------------------------------------------------------- c REAL*8 RINV,RINV1,DUM0,DUM1,DUM2,X,H DIMENSION AX(I2,1),AY(I2,1),BB(I2,1),CX(I2,1),CY(I2,1), 1 RINV(I2,I2,1),RINV1(I2,I2,1),H(I0,1),IE(1),DUM0(I2,1),DUM1(1), 2 DUM2(1),F(I2,1),X(I0,1) JS=1 DO 150 NB=1,NBLK JF=IE(NB)-2 DO 105 J=JS,JF DO 105 I=1,I2 105 X(I+1,J+2)=(F(I,J)-AX(I,J)*X(I,J+1)-AY(I,J)*X(I+1,J)-BB(I,J)* 1 X(I+1,J+1)-CX(I,J)*X(I+2,J+1))/CY(I,J) IF (NB.EQ.NBLK) GO TO 150 J=IE(NB)-1 DO 115 I=1,I2 115 DUM1(I)=F(I,J)-AX(I,J)*X(I,J+1)-AY(I,J)*X(I+1,J)-BB(I,J)* 1 X(I+1,J+1)-CX(I,J)*X(I+2,J+1)-CY(I,J)*X(I+1,J+2) J=IE(NB) DO 120 N=1,I2 DUM2(N)=0. DO 118 M=1,I2 118 DUM2(N)=DUM2(N)+DUM1(M)*RINV1(M,N,NB) DUM0(N,NB)=X(N+1,J) 120 X(N+1,J)=X(N+1,J)-DUM2(N) 150 JS=IE(NB) DO 300 NBS=1,NBLK NB=NBLK-NBS+1 JS=1 IF (NB.NE.1) JS=IE(NB-1) JF=IE(NB)-2 IF (NB.EQ.NBLK) GO TO 201 J=IE(NB) DO 200 N=1,I2 200 X(N+1,J)=DUM0(N,NB) 201 N=IE(NB) DO 202 J=JS,N DO 202 I=1,I0 202 H(I,J)=0. J=IE(NB)-1 DO 210 I=1,I2 210 DUM1(I)=F(I,J)-AX(I,J)*X(I,J+1)-AY(I,J)*X(I+1,J)-BB(I,J)* 1 X(I+1,J+1)-CX(I,J)*X(I+2,J+1)-CY(I,J)*X(I+1,J+2) DO 220 N=1,I2 DUM2(N)=0. DO 218 M=1,I2 218 DUM2(N)=DUM2(N)+DUM1(M)*RINV(M,N,NB) H(N+1,JS+1)=DUM2(N) 220 X(N+1,JS+1)=X(N+1,JS+1)+DUM2(N) IF (NB.EQ.1) GO TO 250 DO 230 M=1,I2 230 DUM1(M)=H(M+1,JS+1)*CY(M,JS-1) J=IE(NB-1) DO 240 N=1,I2 DUM2(N)=0. DO 238 M=1,I2 238 DUM2(N)=DUM2(N)+DUM1(M)*RINV1(M,N,NB-1) 240 H(N+1,J)=DUM2(N) 250 DO 300 J=JS,JF DO 299 I=1,I2 H(I+1,J+2)=(-AX(I,J)*H(I,J+1)-AY(I,J)*H(I+1,J)-BB(I,J)* 1 H(I+1,J+1)-CX(I,J)*H(I+2,J+1))/CY(I,J) 299 X(I+1,J+2)=X(I+1,J+2)+H(I+1,J+2) 300 CONTINUE END C **************** C INITIALIZATION * C **************** C ---------------------------------------------------------------------- SUBROUTINE INITFS C INITIALIZE ALL CONTROLLING ARRAYS AND DERIVED SCALARS FOR FS C (SCALAR CONTROL PARAMETERS ARE INITIALIZED IN BLOCK DATA PROGRAM AT END C OF THIS FILE) C ---------------------------------------------------------------------- INCLUDE 'resolution.h' PARAMETER(IJ0=I0+J0,I1=I0-1,I2=I0-2,J1=J0-1,J2=J0-2, 1 K1=K0-1,K2=K0-2,NB1=NB0-1) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER CASE*5,DSCRIB*63 INTEGER*2 KB,IU0,IV0,IN,IU,IV,IW,ileft,irite COMMON/AGRID/ RHO(I0,J0,K1),PX(I0,J0),PY(I0,J0), 1 U1(I0,J0,K1),U2(I0,J0,K1),V1(I0,J0,K1),V2(I0,J0,K1), 2 S1(I0,J0,K1),S2(I0,J0,K1),T1(I0,J0,K1),T2(I0,J0,K1),P(I0,J0,K1), 3 ULF(I0,J0,K1),VLF(I0,J0,K1),SLF(I0,J0,K1),TLF(I0,J0,K1), 4 EV(I2,J2,K2),HV(I2,J2,K2),DAVG(K1),F(J2),TANPHI(J2),SUMIN(K1) COMMON/CGRID/U(I0,J0,K1),V(I0,J1,K1),W(I0,J0,K0) COMMON/SURF/SSURF(I2,J2),TSURF(I2,J2) COMMON/WINDS/TAUX(I2,J2),TAUY(I2,J2) COMMON/OPENB/AROUT common/bloks/ileft(j0,k1),irite(j0,k1) COMMON/CONTRO/ 1 DT,TLZ,B,G,DM0,DE0,DMZ0,RXYMX,RZMX,TAU,DRAG,FLTW,AV,DAODT,TAUDT, 2 TAUDTN,KTRM,M1,M2,M3,M4,M5,M5M,M6,M7,M8,KF,MVI,MXIT,ISAV,mstart, 3 MXSAV,LRSTRT,LOPEN,LSURF,LTURB,LWIND,LHEAT,LSOL,LSAVE,LMOVI COMMON/SPONGE/DM0X(I0,J0) COMMON/KLIMAT/SCR(I0,J0,K1),SC1(IJ0),SC2(IJ0),TW(I0,J0,2), 1 TCLI(I0,J0,K0),UCLI(I0,J0),VCLI(I0,J0),RMSV(I0,J0), 2 XCLI(I0,J0),PCLI(I0,J0,2),HEATOP(I0,J0),SALTOP(I0,J0) CTS NESTING BOUNDARY COMMON/BCS/ 1 U2WST(J0,K1),V2WST(J0,K1),S2WST(J0,K1),T2WST(J0,K1), 3 U2SUD(I0,K1),V2SUD(I0,K1),S2SUD(I0,K1),T2SUD(I0,K1), 4 U2NOR(I0,K1),V2NOR(I0,K1),S2NOR(I0,K1),T2NOR(I0,K1), 5 UA2WST(J0,K1),VA2SUD(I0,K1),VA2NOR(I0,K1) COMMON/SCA/ODT,PRN,DAYS,ORXYMX,ORZMX,OFLTW,OM5,OM8,ITF,IT0 common/osave/ daodtsave,dtsave,odtsave,itfsave COMMON/METRICS/ 1 Y(J0),YV(J0),YVDEG(J0),YDEG(J0), 2 CS(J0),CSV(J1),OCS(J0),OCSV(J1), 3 DX(J0),DXV(J1),ODX(J0),ODXV(J1), 4 DY(J0),DYV(J1),ODY(J0),ODYV(J1) COMMON/SEVP/RINV(I2,I2,NB0),RINV1(I2,I2,NB1),DUM0(I2,NB1), 1 DUM1(I2),DUM2(I2),X(I0,J0),H(I0,J0),AL(I2,J2),AB(I2,J2), 2 AC(I2,J2),AR(I2,J2),AT(I2,J2),S(I2,J2),IE(NB0), 3 ZBOTU(I1,J0),ZBOTV(I0,J1) COMMON/TOPOG/KB(I0,J0),IU0(I0,J0),IV0(I0,J0), 1 IN(I0,J0,K1),IU(I0,J0,K1),IV(I0,J1,K1),IW(I0,J0,K0) COMMON/ZFS/Z(K0+K1),ODZ(K1),ODZW(K0) COMMON/TITLE/CASE,DSCRIB COMMON/TIDES/ILOC(J0),WBDAYS(J2,360),WBSTPS(180) COMMON/TRNS/DTEMP(N0,K1),VKE(N0,K1),THST(N0),PXT(I0,J0) common/bdytsuv/ vnor(i0,k1), > uwst(j0,k1),vsud(i0,k1) COMMON/CALCUR/UWSTM(K1),VNORM(K1),VSUDM(K1) C 1 FORMAT('CASE ',A5,' (',A63,')'/'----------'/ 1 'SW corner of scalar cell at ( 2, 2) location is (X0DEG,Y0DEG)'/ 2 'NW corner of scalar cell at (I1,J1) location is (X1DEG,Y1DEG)'/ 3 ' LRSTRT MXIT I0 J0 K0 KTRM LSAVE'/2I7,5I5// 4 8X,'DT',6X,'FLTW',7X,'TLZ'/1P,3(1X,E9.2)//5X, 5 'X0DEG',5X,'X1DEG',5X,'Y0DEG',5X,'Y1DEG',5X,'DXDEG'/ 6 0P,5F10.3//7X,'TAU',9X,'G',7X,'DM0',7X,'DE0',6X,'DMZ0',5X, 7 'RXYMX',6X,'RZMX',6X,'DRAG'/1P,8(1X,E9.2)) 2 FORMAT( 1 /'((KB(I,J),I=2,I1),J=J1,2,-1) (INTERIOR LOGICAL DEPTH ARRAY)' 2 /(20I4)) 6 FORMAT(2I4/4F13.7,A10) 8 FORMAT('PSI(SV)'/(10(1X,F7.3))) 9 FORMAT(5E15.7) CTS ************************************** CTS *** OPEN AND READ INPUT DATA FILES *** CTS ************************************** CTS CTS ------------------------------------ CTS NESTING BOUNDARY DATA FROM MED MODEL CTs ------------------------------------ if(LRSTRT.eq.0)THEN OPEN(81,file='PREP/DATA/mon_ics.dat',form='unformatted') REWIND 81 ENDIF OPEN(91,file='PREP/DATA/mon_bcs.dat',form='unformatted') REWIND 91 CTS KB boundary values from Med model CTS This should agree with full array read below if(LRSTRT.eq.0)THEN READ(81) (KB(2,J),J=2,J1) READ(81) (KB(I1,J),J=2,J1) READ(81) (KB(I,2),I=2,I1) READ(81) (KB(I,J1),I=2,I1) C Read initial conditions C The data is stored as real instead of real*8 DO K=1,K1 READ(81) ((U1(I,J,K),I=2,I1),J=2,J1), 1 ((U2(I,J,K),I=2,I1),J=2,J1), 2 ((ULF(I,J,K),I=2,I1),J=2,J1) READ(81) ((V1(I,J,K),I=2,I1),J=2,J1), 1 ((V2(I,J,K),I=2,I1),J=2,J1), 2 ((VLF(I,J,K),I=2,I1),J=2,J1) READ(81) ((S1(I,J,K),I=2,I1),J=2,J1), 1 ((S2(I,J,K),I=2,I1),J=2,J1), 2 ((SLF(I,J,K),I=2,I1),J=2,J1) READ(81) ((T1(I,J,K),I=2,I1),J=2,J1), 1 ((T2(I,J,K),I=2,I1),J=2,J1), 2 ((TLF(I,J,K),I=2,I1),J=2,J1) ENDDO CLOSE(81) ENDIF C ----------------------------------- C SCALAR DATA DEFINING RUN TO BE MADE C ----------------------------------- c cycle reads in monthly fields for use in daily spline interpolation call cycle OPEN(9,file='PREP/DATA/RUNDATA',form='unformatted') REWIND 9 READ(9) CASE,DSCRIB,DT,TLZ,B,G,DM0,DE0,DMZ0,EV,HV,F, 1 TANPHI,RXYMX,RZMX,TAU,TAUN,DRAG,FLTW,DAODT,X0DEG,Y0DEG,DXDEG, 2 X1DEG,Y1DEG,ODT,PRN,ORXYMX,ORZMX,OFLTW,Y,YV,YVDEG,YDEG, 3 CS,CSV,OCS,OCSV,DX,DXV,ODX,ODXV,DY,DYV,ODY,ODYV,OM5,OM8, 4 KTRM,LRSTRT,istart,iend,ISAV,MXSAV,M1,M2,M3,M4,M5,M5M,M6,M7, 5 M8,M8M,LOPEN,LSURF,LTURB,LWIND,LHEAT,LSOL,LSAVE,LMOVI,MVI read(9) ileft,irite idaodt = int(daodt) c iend = 155 mxit = idaodt*iend mstart = idaodt*istart itf = mstart days = float(itf)/daodt write(6,*) ' istart,iend,mstart,mxit ',istart,iend,mstart,mxit write(6,*) 'itf days ',itf,days write(6,*) ' dt,daodt,odt,fltw,ofltw ',dt,daodt,odt,fltw,ofltw C Open boundary data IF(LRSTRT.EQ.0)OPEN(28,file='XYPLOTS/DATA_0',form='unformatted') do 11 j=1,j0 do 11 i=1,i0 11 dm0x(i,j) = dm0 c optional sponge this increases diffusion to about 1.5 timest c 1.03 = 2 times CTS REMOVE THE SPONGE LAYERS CTS C ------------------------------------------ C SURFACE SALINITY AND TEMPERATURE RESTORING C AND WIND STRESS DATA C ------------------------------------------ C When LWIND=0, we still have TAUX,TAUY written as place holders C SSURF,TSURF are needed to avoid overwrites of new restoring data C when reading restart OPEN(60,file='PREP/DATA/ZKB',form='unformatted') REWIND 60 READ(60) TLZ,Z,ODZ,ODZW,KB,IN,IU,IV,IW,IU0,IV0 c compute arout ARout=0. DO 414 K=1,K1 DO 412 I=2,I1 TMPs=IN(I,2,K)/(ODXV(1)*ODZ(K)) tmpn = IN(i,j1,k)/(ODXV(j1)*odz(k)) ARout=ARout+TMPs+tmpn 412 continue do 413 j=2,j1 tmpw = IN(2,j,k)/(ody(j)*odz(k)) arout = arout + tmpw 413 continue 414 continue write(6,*) ' ARout correction area ',arout c c call forcing to define boundary forcing at initial time c note forcing will define taux,tauy,ssurf,tsurf,t2,s2,u1,v1 c in ghost zones call forcin(daodt,g,iday) DO 5551 J=2,J1 DO 5551 I=2,I1 SSURF(I-1,J-1)=S1(I,J,1) 5551 TSURF(I-1,J-1)=T1(I,J,1) c for fresh start ensure consistent variables do k=1,k1 do j=1,j0 do i=1,i0 slf(i,j,k) = s2(i,j,k) s1(i,j,k) = s2(i,j,k) tlf(i,j,k) = t2(i,j,k) t1(i,j,k) = t2(i,j,k) enddo enddo enddo OPEN(99,file='PREP/DATA/SEVP',form='unformatted') REWIND 99 READ(99) ZBOTU,ZBOTV,AL,AR,AB,AT,AC,RINV,RINV1,IE CLOSE(99) C ************************** C *** NO MORE INPUT DATA *** C ************************** C 16 IF (LMOVI.EQ.0) GO TO 19 C --------------------------------- C OPTIONAL COMPUTER-GENERATED MOVIE C --------------------------------- C XY animation: MOVIE/MV0 and MOVIE/DATA (LOGICAL UNITS 21 AND 23) OPEN(21,file='MOVIE/DATA',form='unformatted') OPEN(23,file='MOVIE/MV0') REWIND 21 REWIND 23 TMP=MVI*DT WRITE(23,17) CASE,I0,J0,K0,TMP,Z 17 FORMAT(A5,3I3,E11.4/(7E11.4)) WRITE(23,18) ((IN(I,J,1),I=2,I1),J=2,J1) WRITE(23,18) ((IN(I,J,KTRM),I=2,I1),J=2,J1) 18 FORMAT(/(90I1)) CLOSE(23) C C Write formatted run description and parameters 19 WRITE(14,1) CASE,DSCRIB,LRSTRT,MXIT,I0,J0,K0,KTRM,LSAVE,DT,FLTW, 1 TLZ,X0DEG,X1DEG,Y0DEG,Y1DEG,DXDEG,TAU,G,DM0,DE0, 2 DMZ0,RXYMX,RZMX,DRAG DO 21 K=1,K1 SUMIN(K)=0. DO 20 J=2,J1 DO 20 I=2,I1 20 SUMIN(K)=SUMIN(K)+IN(I,J,K) 21 SUMIN(K)=1./SUMIN(K) FLTW=.5*FLTW TAUDT=1./(TAU*DAODT) TAUDTN=1./(TAUN*DAODT) DO 32 K=1,K1 DO 32 J=1,J0 DO 32 I=1,I0 SLF(I,J,K)=S1(I,J,K) TLF(I,J,K)=T1(I,J,K) S2(I,J,K)=S1(I,J,K) 32 T2(I,J,K)=T1(I,J,K) IF (LRSTRT.EQ.0) 1WRITE(28) CASE,DSCRIB,LRSTRT,MXIT,I0,J0,K0,IN,ODX,ODY,ODZ,Z WRITE(30) CASE,DSCRIB,LRSTRT,MXIT,I0,J0,K0,IN,ODX,ODY,ODZ,Z WRITE(31) CASE,DSCRIB,LRSTRT,MXIT,I0,J0,K0,IN,ODX,ODY,ODZ,Z WRITE(32) CASE,DSCRIB,LRSTRT,MXIT,I0,J0,K0,IN,ODX,ODY,ODZ,Z C Depth READ(60) ((SCR(I,J,1),I=1,I0),J=1,J0) CLOSE(60) DO 45 J=2,J1 DO 45 I=2,I1 C Convert to kilometers 45 SCR(I,J,1)=1.E-5*SCR(I,J,1) CALL XYPLOT('Z','depth (kilometers) ',SCR,1,1) C C Calculate zonally averaged Levitus Climatology DO 50 K=1,K1 DO 50 J=2,J1 TMP=0. TEMP=0. INSUM=0. DO 49 I=2,I1 INSUM=INSUM+IN(I,J,K) TMP=TMP+IN(I,J,K)*T1(I,J,K) 49 TEMP=TEMP+IN(I,J,K)*S1(I,J,K) SCR(1,J,K)=0. SCR(2,J,K)=0. IF (INSUM.EQ.0) GO TO 50 SCR(1,J,K)=TMP/INSUM SCR(2,J,K)=TEMP/INSUM 50 CONTINUE CTS INITIALIZE FOR NESTING READ(91) DAY0 WRITE(*,90) DAY0 90 FORMAT('Nested model INITIAL time is day',F9.1,' of Calcur model') IF (LRSTRT.EQ.0) GO TO 170 C ----------------- C READ RESTART DATA C ----------------- READ(19) ITF,DAYS,AV,X,U1,U2,V1,V2,S1,S2,T1,T2,P,ULF,VLF,SLF,TLF, 1 U,V,W,EV,HV,TW,TCLI,UCLI,VCLI,RMSV,HEATOP,SALTOP,XCLI,PCLI,PXT, 2 VNORM,UWSTM,VSUDM WRITE(*,92) ITF,DAYS write(6,*) ' dt,daodt ',dt,daodt WRITE(14,92) ITF,DAYS 92 FORMAT('RESTART DATA READ FOR TIME STEP ',I6,', DAY ',F8.2) I=DAYS*DAODT+.5 IF (I.EQ.ITF) GO TO 99 c time step has been changed itfsave = itf daodtsave = float(itf)/days write(6,*) ' old new daodt ',daodtsave,daodt c now i need to interpolate time levels u1,u2 to reflect new time c step write(6,*) ' restart u,v,s,t ' write(6,*) u1(70,70,1),v1(70,70,1),s1(70,70,1),t1(70,70,1) write(6,*) ulf(70,70,1),vlf(70,70,1),slf(70,70,1),tlf(70,70,1) write(6,*) u2(70,70,1),v2(70,70,1),s2(70,70,1),t2(70,70,1) odaodtsave = 1./daodtsave odaodt = 1./daodt c goto 7161 do k=1,k1 do j=1,j0 do i=1,i0 u1(i,j,k)=u1(i,j,k)+(1.-daodtsave/daodt)*(ulf(i,j,k)-u1(i,j,k)) v1(i,j,k)=v1(i,j,k)+(1.-daodtsave/daodt)*(vlf(i,j,k)-v1(i,j,k)) s1(i,j,k)=s1(i,j,k)+(1.-daodtsave/daodt)*(slf(i,j,k)-s1(i,j,k)) t1(i,j,k)=t1(i,j,k)+(1.-daodtsave/daodt)*(tlf(i,j,k)-t1(i,j,k)) enddo enddo enddo c note that ulf and u2 are at same time level so no c adjustment is produced by time interpolation do k=1,k1 do j=1,j0 do i=1,i0 u2(i,j,k) = ulf(i,j,k) + odaodt/odaodtsave*(u2(i,j,k)-ulf(i,j,k)) v2(i,j,k) = vlf(i,j,k) + odaodt/odaodtsave*(v2(i,j,k)-vlf(i,j,k)) s2(i,j,k) = slf(i,j,k) + odaodt/odaodtsave*(s2(i,j,k)-slf(i,j,k)) t2(i,j,k) = tlf(i,j,k) + odaodt/odaodtsave*(t2(i,j,k)-tlf(i,j,k)) enddo enddo enddo 7161 continue write(6,*) ' after adjustment u,v,s,t u1,ulf,u2' write(6,*) u1(70,70,1),v1(70,70,1),s1(70,70,1),t1(70,70,1) write(6,*) ulf(70,70,1),vlf(70,70,1),slf(70,70,1),tlf(70,70,1) write(6,*) u2(70,70,1),v2(70,70,1),s2(70,70,1),t2(70,70,1) write(6,*) 'dt,daodt,odt,itf,days',dt,daodt,odt,itf,days 99 CONTINUE C C Calculate model zonally averaged climatology C and differences from Levitus climatology RMST=0. RMSS=0. DO 150 K=1,K1 VOLSUM=0. DO 150 J=2,J1 TMP=0. TEMP=0. INSUM=0. DO 148 I=2,I1 INSUM=INSUM+IN(I,J,K) TMP=TMP+IN(I,J,K)*T2(I,J,K) 148 TEMP=TEMP+IN(I,J,K)*S2(I,J,K) VOL=INSUM/(ODX(J)*ODZ(K)) VOLSUM=VOLSUM+VOL SCR(3,J,K)=0. SCR(4,J,K)=0. IF (INSUM.EQ.0) GO TO 149 SCR(3,J,K)=TMP/INSUM SCR(4,J,K)=TEMP/INSUM 149 SCR(5,J,K)=SCR(3,J,K)-SCR(1,J,K) SCR(6,J,K)=SCR(4,J,K)-SCR(2,J,K) RMST=RMST+VOL*SCR(5,J,K)**2 RMSS=RMSS+VOL*SCR(6,J,K)**2 150 CONTINUE RMSS=SQRT(RMSS/VOLSUM) RMST=SQRT(RMST/VOLSUM) WRITE(14,151) RMST,RMSS 151 FORMAT('Volume weighted full basin rms difference between'/ 1' zonally averaged model T'/ 2'and zonally averaged Levitus annual mean T=',F6.3,' deg. C.', 3//'Volume weighted full basin rms difference between'/ 4' zonally averaged model S'/ 5'and zonally averaged Levitus annual mean S=',F6.3,' ppt.') CALL YZPLOT('S','mean zonal Levitus salt (ppt) ',SCR(2,1,1),1,1) CALL YZPLOT('T','mean zonal Levitus T (deg. C) ',SCR,1,1) CALL YZPLOT('S','mean zonal S-S(Levitus) (ppt) ',SCR(6,1,1),1,1) CALL YZPLOT('T','mean zonal T-T(Levitus)(deg.C)',SCR(5,1,1),1,1) WRITE(14,155) ((SCR(5,J,K),K=1,K1),J=2,J1) 155 FORMAT('((T-Tlev (K,J),K=1,K1),J=2,J1)'/(10(1X,F3.0))) C RETURN CTS ADD BOUNDARY CTS Read boundary values at initial time CTS Nested model DAYS is relative to first day of Cal model data saves 170 NREC=DAYS+1 print*,'NREC',NREC DO 175 N=1,NREC print*,'in',N READ(91) ((U(1,J,K),J=2,J1),K=1,K1) READ(91) ((U1(1,J,K),J=2,J1),K=1,K1) READ(91) ((V1(1,J,K),J=2,J1),K=1,K1) READ(91) ((S1(1,J,K),J=2,J1),K=1,K1) READ(91) ((T1(1,J,K),J=2,J1),K=1,K1) READ(91) ((V(I,1,K),I=2,I1),K=1,K1) READ(91) ((U1(I,1,K),I=2,I1),K=1,K1) READ(91) ((V1(I,1,K),I=2,I1),K=1,K1) READ(91) ((S1(I,1,K),I=2,I1),K=1,K1) READ(91) ((T1(I,1,K),I=2,I1),K=1,K1) READ(91) ((V(I,J1,K),I=2,I1),K=1,K1) READ(91) ((U1(I,J0,K),I=2,I1),K=1,K1) READ(91) ((V1(I,J0,K),I=2,I1),K=1,K1) READ(91) ((S1(I,J0,K),I=2,I1),K=1,K1) READ(91) ((T1(I,J0,K),I=2,I1),K=1,K1) print*,'out',N 175 IF (N.NE.NREC) READ(91) TMP WRITE(*,176) TMP 176 FORMAT('Nested model RESTART time is day',F9.1,' of Calcu model.'/ 1'INIT time is time nested model is initialized by Cal model data') DO 23 K=1,K1 TEMP=0. TMP1=0. UWSTM(K)=0. DO 22 J=2,J1 TMP2=IN(2,J,K)/ODY(J) TMP1=TMP1+TMP2 22 TEMP=TEMP+U(1,J,K)*TMP2 IF (TMP1.NE.0.) TMP1=TEMP/TMP1 23 UWSTM(K)=TMP1 CTS Calculate mean southern boundary V DO 27 K=1,K1 TEMP=0. TMP1=0. VSUDM(K)=0. DO 25 I=2,I1 TMP2=IN(I,2,K)/ODXV(1) TMP1=TMP1+TMP2 25 TEMP=TEMP+V(I,1,K)*TMP2 IF (TMP1.NE.0.) TMP1=TEMP/TMP1 VSUDM(K)=TMP1 TEMP=0. TMP1=0. VNORM(K)=0. DO 26 I=2,I1 TMP2=IN(I,J1,K)/ODXV(J1) TMP1=TMP1+TMP2 26 TEMP=TEMP+V(I,J1,K)*TMP2 IF (TMP1.NE.0.) TMP1=TEMP/TMP1 VNORM(K)=TMP1 27 CONTINUE WRITE(14,24) UWSTM WRITE(*,24) UWSTM 24 FORMAT('vertical profile of latitudinal avg western U'/(10F7.3)) WRITE(14,28) VNORM WRITE(*,28) VNORM 28 FORMAT('vertical profile of longitudinal avg northern V'/(10F7.3)) WRITE(14,29) VSUDM WRITE(*,29) VSUDM 29 FORMAT('vertical profile of longitudinal avg southern V'/(10F7.3)) C C ----------------------------------------- C SHOW INITIAL FIELDS WHEN STARTING NEW RUN C ----------------------------------------- C DO 440 J=1,J2 DO 440 I=1,I2 SCR(I+1,J+1,1)=SSURF(I,J) 440 SCR(I+1,J+1,2)=TSURF(I,J) CALL XYPLOT('S','surface restoring salinty(ppt)',SCR,1,1) CALL XYPLOT('T','surface restoring temperature ',SCR(1,1,2),1,1) K=2 IF (LRSTRT.NE.1) K=K1 CALL XYPLOT('S','initial salinity (ppt) ',S2,2,K) CALL XYPLOT('T','initial temperature(degrees c)',T2,2,K) C Wind stress DO 450 J=2,J1 DO 450 I=2,I1 C Needed for vector plots C IF (TAUX(I-1,J-1).EQ.0.) TAUX(I-1,J-1)=1.E-4 C IF (TAUY(I-1,J-1).EQ.0.) TAUY(I-1,J-1)=1.E-4*(I-2)*(I-I1)/(I2*I2) SCR(I,J,4)=TAUX(I-1,J-1) SCR(I,J,5)=TAUY(I-1,J-1) 450 SCR(I,J,3)=SQRT(TAUX(I-1,J-1)**2+TAUY(I-1,J-1)**2) DO 451 J=3,J2 DO 451 I=3,I2 451 SCR(I,J,6)=.5*((TAUY(I,J-1)-TAUY(I-2,J-1))*ODX(J)) 1 -(TAUX(I-1,J)-TAUX(I-1,J-2))*ODY(J) C Smooth the wind stress curl DO 456 K=1,3 DO 453 J=2,J1 SCR(2,J,6)=SCR(3,J,6) 453 SCR(I1,J,6)=SCR(I2,J,6) DO 454 I=2,I1 SCR(I,2,6)=SCR(I,3,6) 454 SCR(I,J1,6)=SCR(I,J2,6) DO 455 J=3,J2 DO 455 I=3,I2 455 SCR(I,J,6)= 1 .25*(SCR(I+1,J,6)+SCR(I-1,J,6)+SCR(I,J+1,6)+SCR(I,J-1,6)) 456 CONTINUE C Use scale factor to get reasonable curl magnitude TMP=100./ODX(J0/2) DO 457 J=2,J1 DO 457 I=2,I1 457 SCR(I,J,6)=TMP*SCR(I,J,6) CALL XYPLOT('M','wind stress magnitude ',SCR(1,1,3),1,1) CALL XYPLOT('U','long. wind stress(dynes/cm-cm)',SCR(1,1,4),1,1) CALL XYPLOT('V','lat. wind stress (dynes/cm-cm)',SCR(1,1,5),1,1) CALL XYPLOT('N','scaled wind stress curl ',SCR(1,1,6),1,1) 459 WRITE(14,460) (F(J),J=J2,1,-1) 460 FORMAT( 1 /'F(J),J=J2,1,-1) (Latitudinal profile of Coriolis parameter)' 2 /1P,(8(1X,E9.2))) WRITE(14,461) CS 461 FORMAT('METRICS'/'COSLAT at cell centers'/0P,(10(1X,F5.3))) WRITE(14,462) CSV 462 FORMAT('COSLAT at grid lines'/0P,(10(1X,F5.3))) WRITE(14,463) DX 463 FORMAT('DX at cell centers'/1P,(8(1X,E9.2))) WRITE(14,464) DXV 464 FORMAT('DXV at grid lines'/1P,(8(1X,E9.2))) WRITE(14,480) (Z(2*K-1),K=1,K0) 480 FORMAT(/'Interface depths (meters)',-2P/(10(1X,F7.2))) WRITE(14,482) (Z(2*K),K=1,K1) 482 FORMAT(/'Cell center depths (meters)',-2P/(10(1X,F7.2))) WRITE(14,490) (ODZ(K),K=1,K1) 490 FORMAT(/'Inverse grid increments at cell centers (1/cm)'/ 1 (8(1X,1PE9.2))) WRITE(14,492) (ODZW(K),K=1,K0) 492 FORMAT(/'Inverse grid increments at cell interfaces (1/cm)'/ 1 (8(1X,1PE9.2))) C C Barotropic boundary psi PSI=0. C South DO 503 I=I1,2,-1 KM=KB(I,2) DO 502 K=1,KM 502 PSI=PSI-V(I,1,K)*DXV(1)/ODZ(K) 503 SCR(I,1,1)=1.E-12*PSI WRITE(14,504) WRITE(14,8) (SCR(I,1,1),I=2,I1) 504 FORMAT(/20X,'Southern boundary') C West DO 533 J=2,J1 KM=KB(2,J) DO 532 K=1,KM 532 PSI=PSI-U(1,J,K)*DY(J)/ODZ(K) 533 SCR(1,J,1)=1.E-12*PSI WRITE(14,534) 534 FORMAT(/20X,'Western boundary') WRITE(14,8) (SCR(1,J,1),J=2,J1) C C North DO 523 I=2,I1 KM=KB(I,J1) DO 522 K=1,KM 522 PSI=PSI+V(I,J1,K)*DXV(J1)/ODZ(K) 523 SCR(I,1,1)=1.E-12*PSI WRITE(14,524) 524 FORMAT(20X,'Northern boundary') WRITE(14,8) (SCR(I,1,1),I=2,I1) C C East DO 513 J=2,J1 KM=KB(I1,J) DO 512 K=1,KM 512 PSI=PSI-U(I1,J,K)*DY(J)/ODZ(K) 513 SCR(1,J,1)=1.E-12*PSI WRITE(14,514) 514 FORMAT(/20X,'Eastern boundary') WRITE(14,8) (SCR(1,J,1),J=2,J1) C Bathymmetry C WRITE(14,2) ((KB(I,J),I=2,I1),J=J1,2,-1) WRITE(14,540) 540 FORMAT(/'Surface layer masking array') C DO 542 J=J0,1,-1 C542 WRITE(14,543) (IN(I,J,1),I=2,I1) 543 FORMAT(161I1) C IF (LRSTRT.EQ.1) RETURN AV=0. CALL SETR(XCLI,3*I0*J0,0.) CALL SETR(TCLI,I0*J0*K0+2*I0*J0,0.) END subroutine cycle INCLUDE 'resolution.h' PARAMETER(IJ0=I0+J0,I1=I0-1,I2=I0-2,J1=J0-1,J2=J0-2, 1 K1=K0-1,K2=K0-2) IMPLICIT REAL*8 (A-H,O-Z) INTEGER*2 KB,IU0,IV0,IN,IU,IV,IW COMMON/AGRID/ RHO(I0,J0,K1),PX(I0,J0),PY(I0,J0), 1 U1(I0,J0,K1),U2(I0,J0,K1),V1(I0,J0,K1),V2(I0,J0,K1), 2 S1(I0,J0,K1),S2(I0,J0,K1),T1(I0,J0,K1),T2(I0,J0,K1),P(I0,J0,K1), 3 ULF(I0,J0,K1),VLF(I0,J0,K1),SLF(I0,J0,K1),TLF(I0,J0,K1), 4 EV(I2,J2,K2),HV(I2,J2,K2),DAVG(K1),F(J2),TANPHI(J2),SUMIN(K1) COMMON/CGRID/U(I0,J0,K1),V(I0,J1,K1),W(I0,J0,K0) COMMON/KLIMAT/SCR(I0,J0,K1),SC1(IJ0),SC2(IJ0),TW(I0,J0,2), 1 TCLI(I0,J0,K0),UCLI(I0,J0),VCLI(I0,J0),RMSV(I0,J0), 2 XCLI(I0,J0),PCLI(I0,J0,2),HEATOP(I0,J0),SALTOP(I0,J0) COMMON/SCA/ODT,PRN,DAYS,ORXYMX,ORZMX,OFLTW,OM5,OM8,ITF,IT0 COMMON/METRICS/ 1 Y(J0),YV(J0),YVDEG(J0),YDEG(J0), 2 CS(J0),CSV(J1),OCS(J0),OCSV(J1), 3 DX(J0),DXV(J1),ODX(J0),ODXV(J1), 4 DY(J0),DYV(J1),ODY(J0),ODYV(J1) COMMON/TOPOG/KB(I0,J0),IU0(I0,J0),IV0(I0,J0), 1 IN(I0,J0,K1),IU(I0,J0,K1),IV(I0,J1,K1),IW(I0,J0,K0) COMMON/ZFS/Z(K0+K1),ODZ(K1),ODZW(K0) COMMON/months/ > TSURFt(I2,J2,12),SSURFt(I2,J2,12), > TAUXt(I2,J2,12),TAUYt(I2,J2,12) DIMENSION TSURF(I2,J2),SSURF(I2,J2),TAUX(I2,J2),TAUY(I2,J2) character*3 mon(12) character*15 sfc character*21 sfcin data mon /'jan','feb','mar','apr','may','jun', > 'jul','aug','sep','oct','nov','dec'/ data sfc /'PREP/DATA/stop_'/ do l=1,12 sfcin = sfc//mon(l) open(71,file=sfcin,form='unformatted') read(71) ssurf,tsurf,taux,tauy write(6,*) 'read month ',mon(l) taumax = -1e10 do j=1,j2 do i=1,i2 tsqrt = sqrt(taux(i,j)**2 + tauy(i,j)**2) taumax = max(taumax,tsqrt) enddo enddo write(6,*) 'month maxtau ',mon(l),taumax do i=1,i2 do j=1,j2 tsurft(i,j,l) = tsurf(i,j) ssurft(i,j,l) = ssurf(i,j) tauxt(i,j,l) = taux(i,j) tauyt(i,j,l) = tauy(i,j) enddo enddo close(71) enddo return end c new forcing that allows subroutine forcin(daodt,g,iday) INCLUDE 'resolution.h' PARAMETER(IJ0=I0+J0,I1=I0-1,I2=I0-2,J1=J0-1,J2=J0-2, 1 K1=K0-1,K2=K0-2) IMPLICIT REAL*8 (A-H,O-Z) INTEGER*2 KB,IU0,IV0,IN,IU,IV,IW COMMON/AGRID/ RHO(I0,J0,K1),PX(I0,J0),PY(I0,J0), 1 U1(I0,J0,K1),U2(I0,J0,K1),V1(I0,J0,K1),V2(I0,J0,K1), 2 S1(I0,J0,K1),S2(I0,J0,K1),T1(I0,J0,K1),T2(I0,J0,K1),P(I0,J0,K1), 3 ULF(I0,J0,K1),VLF(I0,J0,K1),SLF(I0,J0,K1),TLF(I0,J0,K1), 4 EV(I2,J2,K2),HV(I2,J2,K2),DAVG(K1),F(J2),TANPHI(J2),SUMIN(K1) COMMON/SCA/ODT,PRN,DAYS,ORXYMX,ORZMX,OFLTW,OM5,OM8,ITF,IT0,ITFDAY COMMON/METRICS/ 1 Y(J0),YV(J0),YVDEG(J0),YDEG(J0), 2 CS(J0),CSV(J1),OCS(J0),OCSV(J1), 3 DX(J0),DXV(J1),ODX(J0),ODXV(J1), 4 DY(J0),DYV(J1),ODY(J0),ODYV(J1) COMMON/TOPOG/KB(I0,J0),IU0(I0,J0),IV0(I0,J0), 1 IN(I0,J0,K1),IU(I0,J0,K1),IV(I0,J1,K1),IW(I0,J0,K0) COMMON/ZFS/Z(K0+K1),ODZ(K1),ODZW(K0) COMMON/SURF/SSURF(I2,J2),TSURF(I2,J2) COMMON/WINDS/TAUX(I2,J2),TAUY(I2,J2) COMMON/months/ > TSURFt(I2,J2,12),SSURFt(I2,J2,12), > TAUXt(I2,J2,12),TAUYt(I2,J2,12) dimension vnorm(k1),vsudm(k1),uwstm(k1) real values(12),vx(18),v(12),x(18),y2(18),t(18),s(18) real uw(18),vw(18),sv,uv,tv,rtime,cx(12),cw(12),camp(12) real cws(12),cxs(12),cds(12),camps(12) real ssv,stv,suv,scv,bnn,bss,bww data cx / 124,124,123,122,121,119,112,106,100,93,83,83/ data cw / 20,20,20,20,20,22,25,30,35,42,50,50/ data camp/ 3,0,-5,-10,-13,-13,-11,-8,-6,-4,-1,2.5/ data cxs /144,218,215,210,202,188,178,168,161,154,148,144/ data cws /60,20,22,25,29,35,40,45,49,52,54,56/ data cds /240,70,70,80,100,125,170,200,210,220,230,240/ data camps/ -4,-4,-10,-11,-10,-9,-8,-7,-6,-5,-4,-4/ write(6,*) ' entered forcin iday ',iday iday = int(days+.002) cl = camp(5) do l=1,12 camp(l) = camp(l)/cl enddo dmon = daodt*30. dyear = dmon*12. iyr = int(dyear) + 0.5 itime = mod(itf,iyr) rtime = float(itime) write(6,*) ' daodtiyr,rtime,dmon',daodt,iyr,itime,dmon do i=1,12 v(i) = dmon/2. + (i-1)*dmon x(i+3) = v(i) enddo x(1) = - v(3) x(2) = - v(2) x(3) = - v(1) x(16) = v(12) + (v(12)-v(11)) x(17) = x(16) + (v(12)-v(11)) x(18) = x(17) + (v(12)-v(11)) c surface stresses do i=1,i2 do j=1,j2 do l=1,12 t(l+3) = tauxt(i,j,l) s(l+3) = tauyt(i,j,l) enddo do l = 0,2 t(16+l) = t(4+l) s(16+l) = s(4+l) t(l+1) = t(13+l) s(l+1) = s(13+l) enddo call spline(x,s,18,1.e31,1.e31,y2) call splint(x,s,y2,18,rtime,sv) call spline(x,t,18,1.e31,1.e31,y2) call splint(x,t,y2,18,rtime,tv) taux(i,j) = tv*in(i+1,j+1,1) tauy(i,j) = sv*in(i+1,j+1,1) enddo enddo c sfc temp and salinity do j=1,j2 do i=1,i2 do l=1,12 t(l+3) = tsurft(i,j,l) s(l+3) = ssurft(i,j,l) enddo do l = 0,2 t(16+l) = t(4+l) s(16+l) = s(4+l) t(l+1) = t(13+l) s(l+1) = s(13+l) enddo call spline(x,s,18,1.e31,1.e31,y2) call splint(x,s,y2,18,rtime,sv) call spline(x,t,18,1.e31,1.e31,y2) call splint(x,t,y2,18,rtime,tv) tsurf(i,j) = tv*in(i+1,j+1,1) ssurf(i,j) = sv*in(i+1,j+1,1) enddo enddo return end subroutine spline(x,y,n,yp1,ypn,y2) integer n,nmax real yp1,ypn,x(n),y(n),y2(n) parameter (nmax=500) integer i,k real p,qn,sig,un,u(nmax) if(yp1.gt..99e30) then y2(1)=0. u(1) =0. else y2(1) =-0.5 u(1) = (3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1) endif do i=2,n-1 sig =(x(i)-x(i-1))/(x(i+1)-x(i-1)) p = sig*y2(i-1)+2. y2(i)=(sig-1.)/p u(i) = (6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) / > (x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p enddo if(ypn.gt..99e30) then qn=0. un=0. else qn=0.5 un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1))) endif y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.) do k=n-1,1,-1 y2(k)=y2(k)*y2(k+1)+u(k) enddo return end subroutine splint(xa,ya,y2a,n,x,y) c IMPLICIT REAL*8 (A-H,O-Z) integer n real x,y,xa(n),y2a(n),ya(n) integer k,khi,klo klo=1 khi=n 1 if(khi-klo.gt.1) then k=(khi+klo)/2 if(xa(k).gt.x) then khi= k else klo=k endif goto 1 endif h=xa(khi)-xa(klo) if(h.eq.0.) pause 'bad xa input in splint' a = (xa(khi)-x)/h b = (x-xa(klo))/h y=a*ya(klo) + b*ya(khi) + > ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6. return end