C C*********************************************************************** C pomgcs.f: new generalized vertical coordinate tranformation C with three test grids (ZPS with more layers in BBL): C IF(SIGORZ.EQ.1) CALL MAKSIG(0) ! makes 'pure' sigma-levels C IF(SIGORZ.EQ.2) CALL MAKZLV(0) ! makes 'pure' z-levels C IF(SIGORZ.EQ.3) CALL MAKZPS(0) ! makes mixed z-level, sig system C C Reference: C Mellor, G. L., S. Hakkinen, T. Ezer and R. Patchen, A generalization of C a sigma coordinate ocean model and an intercomparison of model vertical C grids, In: Ocean Forecasting: Conceptual Basis and Applications, N. Pinardi C and J. D. Woods (Eds.), Springer, 55-72, 2001. C C [Preprint copy is available on POM Web page : C http://www.aos.princeton.edu/WWWPUBLIC/htdocs.pom ] C C --------------------------------------------------------- PROGRAM MAIN C************************************************************************ C TEST WITH SIMPLE FLAT BOTTOMED BASIN * C * C THIS IS A VERSION OF THE THREE DIMENSIONAL, TIME DEPENDENT, * C PRIMITIVE EQUATION, OCEAN MODEL DEVELOPED BY ALLAN BLUMBERG AND ME * C WITH SUBSEQUENT CONTRIBUTIONS BY LEO OEY AND OTHERS. TWO REFERENCES * C ARE: * C * C BLUMBERG, A.F. AND G.L. MELLOR, DIAGNOSTIC AND PROGNOSTIC * C NUMERICAL CIRCULATION STUDIES OF THE SOUTH ATLANTIC BIGHT * C J. GEOPHYS. RES. 88, 4579-4592, 1983. * C * C BLUMBERG, A.F. AND G.L. MELLOR, A DESCRIPTION OF A THREE * C COASTAL OCEAN CIRCULATION MODEL, THREE DIMENSIONAL SHELF * C MODELS, COASTAL AND ESTUARINE SCIENCES, 5, N.HEAPS, ED., * C AMERICAN GEOPHYSICAL UNION, 1987. * C * C IN SUBROUTINE PROFQ THE MODEL MAKES USE OF THE TURBULENCE CLOSURE * C SUB-MODEL MOST RECENTLY DESCRIBED IN: * C * C MELLOR, G.L. AND T. YAMADA, DEVELOPMENT OF A TURBULENCE CLOSURE * C MODEL FOR GEOPHYSICAL FLUID PROBLEMS, REV. GEOPHYS. SPACE PHYS., * C 851-879, 1982. * C * C GEORGE MELLOR, 03/87 * C * C IN THIS VERSION THE HORIZONTAL COORDINATES ARE GENERALIZED CURVI- * C LINEAR ORTHOGANAL COORDINATES AND THE MODEL HAS BEEN CONFIGURED * C TO RUN IN HALF PRECISION (32 BITS ) ON THE CYBER 205. * C TO REDUCE ROUNDOFF ERROR, SALINTIES AND TEMPERATURES ARE REDUCED * C BY 35. AND 10. RESPECTIVELY. 64 BIT PRECISION IS USED LOCALLY * C IN SUBROUTINES DENS (WHEREIN THE T AND S OFFSET HAS BEEN TAKEN * C INTO ACCOUNT) AND PROFT. * C GM, 09/87 * C C THIS VERSION HAS BEEN CONVERTED TO STANDARD FORTRAN 77 BY STEVE * C BRENNER WITH SOME ADDITIONAL CHANGES BY ME IN PROFQ WHICH IS NOW * C SOMEWHAT SIMPLER. CODE RESTORED TO SINGLE PRECISION. * C A USER'S GUIDE IS AVAILABLE: * C C MELLOR, G.L., A USER'S GUIDE FOR A THREE-DIMENSIONAL, * C NUMERICAL OCEAN MODEL. PRINCETON UNIVERSITY REPORT, 1990. * C GM, 04/90 * C * C Previous to this date the code was called pmod.f and all changes * C were recorded in pmod.change. Henceforth, the code will be called * C pom.f and changes will be recorded in pom.change. * C GM, 11/92 * C * C THE LAST CODE CHANGE, as recorded in pom.change, WAS ON June 20, 1996 * C************************************************************************ C INCLUDE 'comblkg.h' INTEGER SIGORZ DIMENSION UTB(IM,JM),VTB(IM,JM),UTF(IM,JM),VTF(IM,JM), 1 ADVX(IM,JM,KM),ADVY(IM,JM,KM),ADVUA(IM,JM),ADVVA(IM,JM), 2 TSURF(IM,JM),SSURF(IM,JM), 2 DRHOX(IM,JM,KM),DRHOY(IM,JM,KM),DRX2D(IM,JM),DRY2D(IM,JM), 3 TCLIM(IM,JM,KM),SCLIM(IM,JM,KM),ADX2D(IM,JM),ADY2D(IM,JM), 4 SWRAD(IM,JM),ZDIS(IM,JM),QNET(JM),QSW(JM),QSURF(JM),TC(JM) DIMENSION ITEM(IM) C-------------------------------------------------------------------- REAL ISPI,ISP2I C NAMELIST/PRMTR/MODE,TIME,DTI,ISPLIT,NREAD,IPRTD1,ISWTCH,IPRTD2, C 1 IDAYS,ISPADV,HORCON,TPRNI,SMOTH DATA PI/3.1416/,SMALL/1.E-10/,TIME0/0.0/,BETA/1.98E-11/ IINT=0 TIME=0. NREAD=0 C DTE=240. ISPLIT=30 c ISPLIT=24 DAYS=360. PRTD1=90. ISPADV=5 SMOTH=0.20 HORCON=0.5 TPRNI=.2 RAMP=0. GRAV=9.806 UMOL=2.E-5 MODE=3 SIGORZ=1 ISK=1 JSK=2 C Flow magnitude VEL=0.2 VEL=0.0 C One can create a 'params' file in a runscript to overwrite parameters. INCLUDE 'params' C-------------------------------------------------------------------------- C Herewith is a list of parameters that are user discretionary. C-------------------------------------------------------------------------- C In comblk.h C IM, JM, KM = grid dimension C In MAIN C MODE = 2; 2-D CALCULATION (BOTTOM STRESS CALCULATED IN ADVAVE) C 3; 3-D CALCULATION (BOTTOM STRESS CALCULATED IN PROFU,V) C 4; 3-D CALCULATION WITH T AND S HELD FIXED C DTE = EXTERNAL TIME STEP C ISPLIT = DTI/DTE C NREAD=0, NO RESTART INPUT FILE; NREAD=1, RESTART C PRTD1 = PRINT INTERVAL IN DAYS; C DAYS = LENGTH OF RUN IN DAYS C ISPADV = STEP INTERVAL WHEREIN EXTERNAL MODE ADVECTIVE TERMS ARE C NOT UPDATED C HORCON = CONSTANT IN SMAGORINSKY HORIZONTAL DIFFUSIVITY C TPRNI = HORIZONTAL DIFFUSIVITY INVERSE PRANDTL NUMBER (see S.R. ADVT) C SMOTH = CONSTANT IN TIME SMOOTHER TO PREVENT SOLUTION SPLITTING C TBIAS, SBIAS = Value set in data statement. For 32 bit arithmetic, C setting S=35. and T=10. will reduce round-off error. C RA, HO, DELH = Island/seamount topographical parameters C C In BCOND C HMAX = maximum depth C C In DEPTH C KL1, KL2 = # of logarithmic grid points in surface and bottom C respectively C C In PROFT C NBC = choice of surface boundary conditions C NTP = choice of radiative penetration parameter. C C In SLPMIN C SLMIN = see subroutine C-------------------------------------------------------------------------- DTI=DTE*FLOAT(ISPLIT) DTE2=DTE*2 DTI2=DTI*2 IEND=DAYS*24*3600/DTI IPRINT=PRTD1*24*3600/DTI IFILE=IPRINT write(6,'('' DAYS ='',F10.2,'' SIGORZ ='',I5)') DAYS,SIGORZ WRITE(6,7030) MODE,DTE,DTI,IEND,ISPLIT,ISPADV,IPRINT,SMOTH, 1 HORCON,TPRNI 7030 FORMAT(//,' MODE =',I3, 1 ' DTE =',F7.2,' DTI =',F9.1,' IEND =',I6, 2 ' ISPLIT =',I6,' ISPADV =',I6,' IPRINT =',I6,/, 3 ' SMOTH =',F7.2,' HORCON =',F7.3,' TPRNI ='F7.3,/) WRITE(6,'('' NREAD ='',I5,//)') NREAD C ISPI=1.0/FLOAT(ISPLIT) ISP2I=1.0/(2.0*FLOAT(ISPLIT)) TBIAS=0. SBIAS=0. IP=IM C---------------------------------------------------------------------- C ESTABLISH PROBLEM CHARACTERISTICS C ****** ALL UNITS IN M.K.S. SYSTEM ****** C F,BLANK AND B REFERS TO FORWARD,CENTRAL AND BACKWARD TIME LEVELS. C---------------------------------------------------------------------- C C --- make horizontal grid C C --- domain: 0-20E,60-75N (Winton & Hallberg 97) ALONMN=0. ALONMX=20. ALATMN=60. ALATMX=75. SHW=5. SLX=4. SLY=2. HMAX=3500. HMIN= 500. C ---------------------------------------------- C --- domain: 0-50E, 0-65N ALONMN= 0. ALONMX=50. ALATMN= 0. ALATMX=65. SHW=10. SLX=10. SLY=10. HMAX=4000. HMIN= 500. WRITE(6,'('' ALONMN, ALONMX ='',2F10.3, 1 '' ALATMN, ALATMX ='',2F10.3)') ALONMN,ALONMX,ALATMN,ALATMX WRITE(6,'('' SHW, SLX, SLY ='',3F10.3)') SHW,SLX,SLY C23456 | | | | | | | xxxxxxx| RE=6365.E3 DXLON=(ALONMX-ALONMN)/(IM-1) DYLAT=(ALATMX-ALATMN)/(JM-1) DO J=1,JM DO I=1,IM ALON(I,J)=ALONMN+DXLON*(I-1) ALAT(I,J)=ALATMN+DYLAT*(J-1) ENDDO ENDDO c CALL PRXY(' LAT ',TIME,ALAT,IP,ISK,JM,JSK,1. ) c CALL PRXY(' LON ',TIME,ALON,IP,ISK,JM,JSK,1. ) C C CC1=HMIN CC2=HMAX-HMIN DO J=1,JM DO I=1,IM XPART1=1.-EXP(-1.*((ALON(I,J)-ALONMN)/SLX)**2) XPART2=1.-EXP(-1.*((ALON(I,J)-ALONMX)/SLX)**2) YPART1=1.-EXP(-1.*((ALAT(I,J)-ALATMN)/SLY)**2) YPART1=1. C Leave room for shelf; SHW is shelf width. YPART2=1.-EXP(-1.*((ALAT(I,J)-ALATMX+SHW)/SLY)**2) H(I,J)=CC1 + CC2*XPART1*XPART2*YPART1*YPART2 C shelf with sloping depth IF(ALAT(I,J).GE.(ALATMX-SHW)) 1 H(I,J)=HMIN-(HMIN-20.)*(ALAT(I,J)-ALATMX+SHW)/SHW C ENDDO ENDDO C closed boundaries H(:,1)=1. H(:,JM)=1. H(1,:)=1. H(IM,:)=1. C C SET UP GRID AND INITIAL CONDITIONS FOR STAND-ALONE TEST RUN DO 12 J=2,JM DO 12 I=2,IM RAD=PI/180. DY(I,J)=(ALAT(I,J)-ALAT(I,J-1))*RE*RAD XLAT =(ALAT(I,J)+ALAT(I,J-1))*.5 DX(I,J)=(ALON(I,J)-ALON(I-1,J))*RE*RAD*COS(RAD*XLAT) COR(I,J)=1.454E-4*SIN(RAD*XLAT) 12 CONTINUE DX(1,:)=DX(2,:) DY(1,:)=DY(2,:) DX(:,1)=DX(:,2) DY(:,1)=DY(:,2) DO 13 J=2,JM DO 13 I=2,IM ART(I,J)=DX(I,J)*DY(I,J) ARU(I,J)=.25*(DX(I,J)+DX(I-1,J))*(DY(I,J)+DY(I-1,J)) 13 ARV(I,J)=.25*(DX(I,J)+DX(I,J-1))*(DY(I,J)+DY(I,J-1)) DO 121 J=1,JM COR(1,J)=COR(2,J) ART(1,J)=ART(2,J) ARU(1,J)=ARU(2,J) ARV(1,J)=ARV(2,J) 121 CONTINUE DO 122 I=1,IM COR(I,1)=COR(I,2) ART(I,1)=ART(I,2) ARU(I,1)=ARU(I,2) ARV(I,1)=ARV(I,2) 122 CONTINUE C PERIOD=(2.0*PI)/COR(IM/2,JM/2)/86400. c C IF(SIGORZ.EQ.1) CALL MAKSIG(0) ! makes 'pure' sigma-levels IF(SIGORZ.EQ.2) CALL MAKZLV(0) !makes 'pure' z-levels IF(SIGORZ.EQ.3) CALL MAKZPS(0) !makes mixed z-level, sig system c write(40) ALON,ALAT,Z,ZZ,H,DX,DY c IF(2.GT.1) STOP C CALL PRXY(' TOPOGRAPHY ',TIME, H,IP,ISK,JM,JSK,1. ) C write(6,'('' fsm at j='',I3)')jm/2 do k=1,km write(6,1459) (k),(int(fsm(i,jm/2,k)+0.1),i=1,im) enddo write(6,'('' fsm at i='',I3)') im/2 do k=1,km write(6,'(I4,1x,120I1)') (k),(int(fsm(im/2,j,k)+0.1),j=1,jm) enddo 1458 format(I4,1x,11I6) 1459 format(I4,1x,65I1) DO K=1,KM-1 DO 35 J=2,JM DUM(1,J,K)=0. DO 35 I=1,IM DVM(I,J,K)=FSM(I,J-1,K)*FSM(I,J,K) 35 CONTINUE DO 40 I=2,IM DVM(I,1,K)=0. DO 40 J=1,JM DUM(I,J,K)=FSM(I-1,J,K)*FSM(I,J,K) 40 CONTINUE DUM(1,1,K)=0. DVM(1,1,K)=0. ENDDO write(6,'('' dum at j='',I3)') jm/2 do k=1,km write(6,1459) (k),(int(dum(i,jm/2,k)+0.1),i=1,im) enddo CALL PRXY(' KB ',TIME, FLOAT(KB ),IP,ISK,JM,JSK,1. ) C C A very empirical specification of the bottom roughness C parameter follows DO 45 J=1,JM DO 45 I=1,IM DT(I,J)=H(I,J) Z0B=.01 CBCMIN=.00250 KBB1=KB(I,J)-1 CBC(I,J)=MAX(CBCMIN,.16/LOG(.5*DZB(I,J,KBB1)/Z0B)**2) 45 CONTINUE C C c CALL PRXY(' TOPOGRAPHY 2 ',TIME, H,IP,ISK,JM,JSK,10.) c CALL PRXY(' FSM ',TIME,FSM(1,1,1),IP,ISK,JM,JSK,1.) C DO 47 J=1,JM DO 47 I=1,IM 47 TPS(I,J)=0.5/SQRT(1./DX(I,J)**2+1./DY(I,J)**2) 1 /SQRT(GRAV*H(I,J))*FSM(I,J,1) C CALL PRXY(' EXT. CFL TIME STEP ',TIME,TPS,IP,ISK,JM,JSK,1.) c CALL PRXY(' COR ',TIME,COR,IP,ISK,JM,1,0.) C--------------------------------------------------------------------- C Set initial conditions. C--------------------------------------------------------------------- DO 15 K=1,KM-1 DO 15 J=1,JM DO 15 I=1,IM TB(I,J,K)=5.+15.*EXP(ZZ(I,J,K)/1000.)-TBIAS SB(I,J,K)=35.0-SBIAS TB(I,J,K)=TB(I,J,K)*FSM(I,J,K) SB(I,J,K)=SB(I,J,K)*FSM(I,J,K) TCLIM(I,J,K)=TB(I,J,K) SCLIM(I,J,K)=SB(I,J,K) UB(I,J,K)=0.0 VB(I,J,K)=0.0 AAM(I,J,K)=HVIS 15 CONTINUE c CALL PRXZ(' Z' ,TIME, Z,IM,ISK,JM,10,20,KM,0.,DT) c CALL PRXZ('ZZ' ,TIME,ZZ,IM,ISK,JM,10,20,KM,0.,DT) c CALL GSIGTOZ(ZZ,ZZ,KB ,TLEV,-1500.,TPS,IM,JM,KM) c CALL PRXY(' ZZ AT 1500 M',TIME,TLEV,IP,ISK,JM,JSK,0.E-3) c CALL DENS(SB,TB,RHO) c CALL GSIGTOZ(ZZ,TB,KB ,TLEV,-1500.,TPS,IM,JM,KM) c CALL PRXY(' T AT 1500 M',TIME,TLEV,IP,ISK,JM,JSK,0.E-3) C In this application, the initial condition for density is C a function of z (cartesian). When this is not so, make sure C that RMEAN has been area averaged before transfer to sigma C coordinates. DO 50 K=1,KM-1 DO 50 J=1,JM DO 50 I=1,IM c 50 RMEAN(I,J,K)=0.0 50 RMEAN(I,J,K)=RHO(I,J,K) DO 16 J=1,JM DO 16 I=1,IM TSURF(I,J)=TB(I,J,1) 16 SSURF(I,J)=SB(I,J,1) DO 51 J=1,JM DO 51 I=1,IM UAB(I,J)=0.0 VAB(I,J)=0. ELB(I,J)=0.0 ETB(I,J)=0.0 AAM2D(I,J)=AAM(I,J,1) 51 CONTINUE DO 52 K=1,KM DO 52 J=1,JM DO 52 I=1,IM Q2B(I,J,K)=1.E-8 Q2LB(I,J,K)=1.E-8 KH(I,J,K)=2.E-4 KMT(I,J,K)=KH(I,J,K) L(I,J,K)=Q2LB(I,J,K)/(Q2B(I,J,K)+SMALL) KQ(I,J,K)=0.2*L(I,J,K)*SQRT(Q2B(I,J,K)) 52 CONTINUE C------------------------------------------------ C SET lateral boundary conditions for use in BCOND. C Here we set uniform velocity in the east and C elevation in the east calculated for uniform flow. C------------------------------------------------ ELE(1)=0. DO 60 J=2,JMM1 UABW(J)=VEL*DUM(1,J,1) UABE(J)=VEL*DUM(IM,J,1) ELE(J)=ELE(J-1)-COR(IMM1,J)*UAB(IMM1,J)/GRAV*DY(IMM1,J) 60 CONTINUE DO 61 J=2,JMM1 ELE(J)=(ELE(J)-ELE(JMM1/2))*FSM(IM,J,1) 61 CONTINUE DO K=1,KM-1 DO J=1,JM TBW(J,K)=TB(1,J,K)*FSM(1,J,K) TBE(J,K)=TB(IM,J,K)*FSM(IM,J,K) SBW(J,K)=SB(1,J,K)*FSM(1,J,K) SBE(J,K)=SB(IM,J,K)*FSM(IM,J,K) ENDDO ENDDO C C------------------------------------------------------------------------ C READ IN GRID DATA AND INITIAL CONDITIONS C------------------------------------------------------------------------ C REWIND(40) C REWIND(60) C READ(40) KBB,Z,ZZ,DZ,DZZ,IMM,JMM,ALON,ALAT,DX,DY,H,COR, C 1 ART,ARU,ARV,TCLIM,SCLIM,RMEAN,TBN,TBE,TBS,SBN,SBE,SBS, C 2 ELN,ELE,ELS,VABN,UABE,VABS,FSM,DUM,DVM TIME0=0. C READ(50) TIME0,UB,VB,UAB,VAB,TB,SB,ELB,Q2B,Q2LB,KM,KH C***************************************************************************** DO 62 I=1,IM DO 62 J=1,JM UA(I,J)=UAB(I,J) VA(I,J)=VAB(I,J) EL(I,J)=ELB(I,J) ET(I,J)=ETB(I,J) ETF(I,J)=ET(I,J) D(I,J)=H(I,J)+EL(I,J) DT(I,J)=H(I,J)+ET(I,J) DRX2D(I,J)=0. DRY2D(I,J)=0. 62 CONTINUE DO 63 K=1,KM-1 DO 63 I=1,IM DO 63 J=1,JM DZ(I,J,K)=DZB(I,J,K) DZF(I,J,K)=DZ(I,J,K) Q2(I,J,K)=Q2B(I,J,K) Q2L(I,J,K)=Q2LB(I,J,K) T(I,J,K)=TB(I,J,K) S(I,J,K)=SB(I,J,K) U(I,J,K)=UB(I,J,K) V(I,J,K)=VB(I,J,K) 63 CONTINUE C CALL BAROPG(DRHOX,DRHOY) DO K=1,KM-1 DO J=1,JM DO I=1,IM DRX2D(I,J)=DRX2D(I,J)+DRHOX(I,J,K) DRY2D(I,J)=DRY2D(I,J)+DRHOY(I,J,K) ENDDO ENDDO ENDDO C C THE FOLLOWING DATA ARE NEEDED FOR A SEAMLESS RESTART. DO 10 N=1,NREAD READ(70) TIME0, 1 WUBOT,WVBOT,AAM2D,UA,UAB,VA,VAB,EL,ELB,ET,ETB,EGB,DZ,DZB,Z,ZZ, 2 UTB,VTB,U,UB,W,V,VB,T,TB,S,SB,RHO,ADX2D,ADY2D,ADVUA,ADVVA, 3 KMT,KH,KQ,L,GH,Q2,Q2B,AAM,Q2L,Q2LB TIME=TIME0 WRITE(6,'('' RESTART FILE READ. TIME ='', F8.0, ''DAYS'')') TIME 10 CONTINUE DO 81 J=1,JM DO 81 I=1,IM D(I,J)=H(I,J)+EL(I,J) 81 DT(I,J)=H(I,J)+ET(I,J) C C CALL PRXYZ(' TEMP. ',TIME,TB,IP,ISK,JM,JSK,KM,0.0E0) C CALL PRXYZ('RHO ',TIME,RHO,IP,ISK,JM,1,KM,0.E0) c CALL PRXZ('Z ' ,TIME,Z ,IM,ISK,JM,10,25,KM,0.,DT,ZZND) c CALL PRXZ('DZ' ,TIME,DZ,IM,ISK,JM,10,25,KM,0.,DT,ZZND) c CALL PRXZ('TB' ,TIME,TB,IM,ISK,JM,10,25,KM,0.,DT,ZZND) c CALL PRXZ('RHO ' ,TIME,RHO ,IM,ISK,JM,10,25,KM,0.,DT,ZZND) c CALL PRXZ('RMEAN' ,TIME,RMEAN,IM,ISK,JM,10,25,KM,0.,DT,ZZND) c CALL PRXZ('DRHOX' ,TIME,DRHOX,IM,ISK,JM,10,25,KM,0.,DT,ZZND) c CALL PRXZ('DRHOY' ,TIME,DRHOY,IM,ISK,JM,10,25,KM,0.,DT,ZZND) C CALL PRXY(' DRX2D ',TIME,DRX2D,IP,ISK,JM,JSK,0.) C CALL PRXY(' DRY2D ',TIME,DRY2D,IP,ISK,JM,JSK,0.) c CALL PRXY(' FREE SURFACE ',TIME, ELB,IP,ISK,JM,JSK,1.E-3) C CALL PRXY(' VAB ',TIME, VAB,IP,ISK,JM,JSK,0.) ISA=0 C c IEND=3 c IPRINT=3 c IFILE=10 DO J=2,JMM1 QNET(J)=25. - 50.*ALAT(IM/2,J)/ALATMX QSW (J)=150. - 110.*ALAT(IM/2,J)/ALATMX TC(J)=27.-29.*ALAT(IM/2,J)/ALATMX ENDDO C Check that total heat flux is nil AREA=0. QSUM=0. DO J=2,JMM1 DO I=2,IMM1 AREA=AREA+DX(I,J)*DY(I,J)*FSM(I,J,1) QSUM=QSUM+QNET(J)*DX(I,J)*DY(I,J)*FSM(I,J,1) ENDDO ENDDO QSUM=QSUM/AREA QNET=QNET-QSUM QSURF=QNET-QSW write(6,'('' TOTAL HEAT FLUX ='',E12.4)') QSUM WRITE(6,'('' J,TC,QNET,QSW,QSURF='',I5,4E12.3)') 1 (J,TC(J),QNET(J),QSW(J),QSURF(J),J=2,JMM1) DO 85 J=2,JMM1 DO 85 I=2,IMM1 WUSURF(I,J)=(2.E-4*SIN(2.*PI*ALAT(I,J)/60.)) 1 *.25*(DVM(I,J+1,1)+DVM(I-1,J+1,1)+DVM(I-1,J,1)+DVM(I,J,1)) WVSURF(I,J)=0. WSSURF(I,J)=0.0 SWRAD(I,J)=-QSW(J)/4.093E6*FSM(I,J,1) 85 CONTINUE C*********************************************************************** C * C BEGIN NUMERICAL INTEGRATION * C * C*********************************************************************** DO 9000 IINT=1,IEND c IPRINT=1 c write(6,'('' IINT ='',I7)') IINT C TIME=DTI*FLOAT(IINT)/86400.+TIME0 RAMP=TIME/PERIOD IF(RAMP.GT.1.) RAMP=1. RAMP=1. DO J=2,JMM1 DO I=2,IMM1 WTSURF(I,J)=QSURF(J)-40.*(T(I,J,1)-TC(J)) WTSURF(I,J)=-WTSURF(I,J)/4.093E6*FSM(I,J,1) ENDDO ENDDO C------------------------------------------------------------------------ C------------------------------------------------------------------------ C SET TIME DEPENDENT, SURFACE AND LATERAL BOUNDARY CONDITIONS. C THE LATTER WILL BE USED IN SUBROUTINE BCOND . C------------------------------------------------------------------------ C INTRODUCE SIMPLE WIND STRESS, VALUE IS NEGATIVE FOR WESTERLY WIND C The following wind stress has been tapered along the boundary to suppress C numerically induced oscilations near the boundary (Jamart&Ozer, C J.G.R.,91, 10621-10631) C IF(MODE.NE.2) THEN CALL ADVCT(ADVX,ADVY) CALL BAROPG(DRHOX,DRHOY) c write(6,'('' got here 1'')') C********************************************************************** C HOR VISC = HORCON*DX*DY*SQRT((DU/DX)**2+(DV/DY)**2 C +.5*(DU/DY+DV/DX)**2) C********************************************************************** C If MODE.EQ.2 then initial values of AAM2D are used. If one wishes C Smagorinsky lateral viscosity and diffusion for an external mode C calculation, then appropiate code can be adapted from that below C and installed after s.n 102 and before s.n. 5000 in subroutine ADVAVE. DO 95 K=1,KM-1 DO 95 J=2,JMM1 DO 95 I=2,IMM1 AAM(I,J,K)=HORCON*DX(I,J)*DY(I,J) 1 *SQRT( ((U(I+1,J,K)-U(I,J,K))/DX(I,J))**2 2 +((V(I,J+1,K)-V(I,J,K))/DY(I,J))**2 3 +.50*(.25*(U(I,J+1,K)+U(I+1,J+1,K)-U(I,J-1,K)-U(I+1,J-1,K)) 4 /DY(I,J) 5 +.25*(V(I+1,J,K)+V(I+1,J+1,K)-V(I-1,J,K)-V(I-1,J+1,K)) 6 /DX(I,J)) **2) 95 CONTINUE IF (SIGORZ.EQ.2) AAM=HVIS AAM=HVIS C -- Form vertical averages of 3-D fields for use in external mode -- DO 96 J=1,JM DO 96 I=1,IM ADX2D(I,J)=0. ADY2D(I,J)=0. DRX2D(I,J)=0. DRY2D(I,J)=0. AAM2D(I,J)=0. 96 CONTINUE DO 100 K=1,KM-1 DO 100 J=1,JM DO 100 I=1,IM ADX2D(I,J)=ADX2D(I,J)+ADVX(I,J,K) ADY2D(I,J)=ADY2D(I,J)+ADVY(I,J,K) DRX2D(I,J)=DRX2D(I,J)+DRHOX(I,J,K) DRY2D(I,J)=DRY2D(I,J)+DRHOY(I,J,K) AAM2D(I,J)=AAM2D(I,J)+AAM(I,J,K)*DZ(I,J,K)*FSM(I,J,K) 1 /(H(I,J)+ET(I,J)) 100 CONTINUE C ---------------------------------------------------------------------- c write(6,'('' got here 3'')') CALL ADVAVE(ADVUA,ADVVA,MODE) DO 87 J=1,JM DO 87 I=1,IM ADX2D(I,J)=ADX2D(I,J)-ADVUA(I,J) 87 ADY2D(I,J)=ADY2D(I,J)-ADVVA(I,J) C ENDIF C------------------------------------------------------------------------ DO 399 J=1,JM DO 399 I=1,IM 399 EGF(I,J)=EL(I,J)*ISPI DO 400 J=2,JM DO 400 I=2,IM UTF(I,J)=UA(I,J)*(D(I,J)+D(I-1,J))*ISP2I 400 VTF(I,J)=VA(I,J)*(D(I,J)+D(I,J-1))*ISP2I C********** BEGIN EXTERNAL MODE *************************************** c write(6,'('' got here 4'')') DO 8000 IEXT=1,ISPLIT c write(6,'('' IEXT,TIME ='',I5,F9.2)') IEXT,TIME DO 405 J=2,JM DO 405 I=2,IM FLUXUA(I,J)=.25*(D(I,J)+D(I-1,J))*(DY(I,J)+DY(I-1,J))*UA(I,J) 405 FLUXVA(I,J)=.25*(D(I,J)+D(I,J-1))*(DX(I,J)+DX(I,J-1))*VA(I,J) C DO 410 J=2,JMM1 DO 410 I=2,IMM1 410 ELF(I,J)=ELB(I,J) 1 -DTE2*(FLUXUA(I+1,J)-FLUXUA(I,J)+FLUXVA(I,J+1)-FLUXVA(I,J)) 2 /ART(I,J) C CALL BCOND(1) C c write(6,'('' got here 5'')') IF(MOD(IEXT,ISPADV).EQ.0) CALL ADVAVE(ADVUA,ADVVA,MODE) C Note that ALPHA = 0. is perfectly acceptable. The value, ALPHA = .225 C permits a longer time step. ALPHA=0.225 DO 420 J=2,JMM1 DO 420 I=2,IM C CORFAC=1. C 1 /(DVM(I,J+1,1)+DVM(I,J,1)+DVM(I-1,J+1,1)+DVM(I-1,J,1)+1.E-10) CORFAC=.25 420 UAF(I,J)=ADX2D(I,J)+ADVUA(I,J) 1 -ARU(I,J)*CORFAC*( COR(I,J)*D(I,J)*(VA(I,J+1)+VA(I,J)) 2 +COR(I-1,J)*D(I-1,J)*(VA(I-1,J+1)+VA(I-1,J)) ) 3 +.25*GRAV*(DY(I,J)+DY(I-1,J))*(D(I,J)+D(I-1,J)) 4 *( (1.-2.*ALPHA)*(EL(I,J)-EL(I-1,J)) 4 +ALPHA*(ELB(I,J)-ELB(I-1,J)+ELF(I,J)-ELF(I-1,J)) ) 5 +DRX2D(I,J) 6 +ARU(I,J)*( WUSURF(I,J)-WUBOT(I,J) ) DO 425 J=2,JMM1 DO 425 I=2,IM 425 UAF(I,J)= 1 ((H(I,J)+ELB(I,J)+H(I-1,J)+ELB(I-1,J))*ARU(I,J)*UAB(I,J) 2 -4.*DTE*UAF(I,J)) 3 /((H(I,J)+ELF(I,J)+H(I-1,J)+ELF(I-1,J))*ARU(I,J)) DO 430 J=2,JM DO 430 I=2,IMM1 C CORFAC=1. C 1 /(DUM(I+1,J,1)+DUM(I,J,1)+DUM(I+1,J-1,1)+DUM(I,J-1,1)+1.E-10) CORFAC=.25 430 VAF(I,J)=ADY2D(I,J)+ADVVA(I,J) 1 +ARV(I,J)*CORFAC*( COR(I,J)*D(I,J)*(UA(I+1,J)+UA(I,J)) 2 +COR(I,J-1)*D(I,J-1)*(UA(I+1,J-1)+UA(I,J-1)) ) 3 +.25*GRAV*(DX(I,J)+DX(I,J-1))*(D(I,J)+D(I,J-1)) 4 *( (1.-2.*ALPHA)*(EL(I,J)-EL(I,J-1)) 4 +ALPHA*(ELB(I,J)-ELB(I,J-1)+ELF(I,J)-ELF(I,J-1)) ) 5 +DRY2D(I,J) 6 + ARV(I,J)*( WVSURF(I,J)-WVBOT(I,J) ) DO 435 J=2,JM DO 435 I=2,IMM1 435 VAF(I,J)= 1 ((H(I,J)+ELB(I,J)+H(I,J-1)+ELB(I,J-1))*VAB(I,J)*ARV(I,J) 2 -4.*DTE*VAF(I,J)) 3 /((H(I,J)+ELF(I,J)+H(I,J-1)+ELF(I,J-1))*ARV(I,J)) CALL BCOND(2) c write(6,'('' got here 6'')') C IF(IEXT.LT.(ISPLIT-2)) GO TO 440 IF(IEXT.EQ.(ISPLIT-2))THEN DO 431 J=1,JM DO 431 I=1,IM 431 ETF(I,J)=.25*SMOTH*ELF(I,J) ENDIF IF(IEXT.EQ.(ISPLIT-1)) THEN DO 432 J=1,JM DO 432 I=1,IM 432 ETF(I,J)=ETF(I,J)+.5*(1.-.5*SMOTH)*ELF(I,J) ENDIF IF(IEXT.EQ.(ISPLIT-0)) THEN DO 433 J=1,JM DO 433 I=1,IM 433 ETF(I,J)=(ETF(I,J)+.5*ELF(I,J))*FSM(I,J,1) ENDIF 440 CONTINUE C C TEST FOR CFL VIOLATION. IF SO, PRINT AND STOP C VAMAX=0. DO 442 J=1,JM DO 442 I=1,IM IF(ABS(VAF(I,J)).GE.VAMAX)VAMAX=ABS(VAF(I,J)) 442 CONTINUE IF(VAMAX.GT.100.) GO TO 9001 C C APPLY FILTER TO REMOVE TIME SPLIT. RESET TIME SEQUENCE. DO 445 J=1,JM DO 445 I=1,IM UA(I,J)=UA(I,J)+.5*SMOTH*(UAB(I,J)-2.*UA(I,J)+UAF(I,J)) VA(I,J)=VA(I,J)+.5*SMOTH*(VAB(I,J)-2.*VA(I,J)+VAF(I,J)) EL(I,J)=EL(I,J)+.5*SMOTH*(ELB(I,J)-2.*EL(I,J)+ELF(I,J)) ELB(I,J)=EL(I,J) EL(I,J)=ELF(I,J) D(I,J)=H(I,J)+EL(I,J) UAB(I,J)=UA(I,J) UA(I,J)=UAF(I,J) VAB(I,J)=VA(I,J) VA(I,J)=VAF(I,J) 445 CONTINUE IF(IEXT.EQ.ISPLIT) GO TO 8000 DO 450 J=2,JM DO 450 I=2,IM EGF(I,J)=EGF(I,J)+EL(I,J)*ISPI UTF(I,J)=UTF(I,J)+UA(I,J)*(D(I,J)+D(I-1,J))*ISP2I 450 VTF(I,J)=VTF(I,J)+VA(I,J)*(D(I,J)+D(I,J-1))*ISP2I C 8000 CONTINUE C------create vertical coordinate ---------------------------------- IF(SIGORZ.EQ.1) CALL MAKSIG(1) ! makes 'pure' sigma-levels IF(SIGORZ.EQ.2) CALL MAKZLV(1) !makes 'pure' z-levels IF(SIGORZ.EQ.3) CALL MAKZPS(1) !makes mixed z-level sigma system C--------------------------------------------------------------------- C END EXTERNAL (2-D) MODE CALCULATION C AND CONTINUE WITH INTERNAL (3-D) MODE CALCULATION C--------------------------------------------------------------------- IF(IINT.EQ.1.AND.TIME0.EQ.0.) GO TO 8200 IF(MODE.EQ.2) GO TO 8200 C--------------------------------------------------------------------- C ADJUST U(Z) AND V(Z) SUCH THAT C VERTICAL AVERAGE OF (U,V) = (UA,VA) C--------------------------------------------------------------------- C DO 299 J=1,JM DO 299 I=1,IM ZDIS(I,J)=0. 299 TPS(I,J)=0.0 DO 300 K=1,KM-1 DO 300 J=1,JM DO 300 I=1,IM ZDIST=.5*(DZ(I,J,K)+DZ(I-1,J,K))*DUM(I,J,K) ZDIS(I,J)=ZDIS(I,J)+ZDIST 300 TPS(I,J)=TPS(I,J)+U(I,J,K) *ZDIST DO 302 K=1,KM-1 DO 302 J=1,JM DO 302 I=1,IM if(dum(i,j,1).gt.0.) 1 U(I,J,K)=U(I,J,K) 2 +(.5*(UTB(I,J)+UTF(I,J))-TPS(I,J))/ZDIS(I,J) U(I,J,K)=U(I,J,K)*DUM(I,J,K) 302 CONTINUE DO 303 J=1,JM DO 303 I=1,IM ZDIS(I,J)=0. 303 TPS(I,J)=0. DO 304 K=1,KM-1 DO 304 J=1,JM DO 304 I=1,IM ZDIST =.5*(DZ(I,J,K)+DZ(I,J-1,K))*DVM(I,J,K) ZDIS(I,J)=ZDIS(I,J)+ZDIST 304 TPS(I,J)=TPS(I,J)+V(I,J,K)*ZDIST DO 306 K=1,KM-1 DO 306 J=1,JM DO 306 I=1,IM if(dvm(i,j,1).gt.0.) 1 V(I,J,K)=V(I,J,K) 2 +(.5*(VTB(I,J)+VTF(I,J))-TPS(I,J))/ZDIS(I,J) V(I,J,K)=V(I,J,K)*DVM(I,J,K) 306 CONTINUE C C---------------------------------------------------------------- C VERTVL INPUT = U,V,DT(=H+ET),ETF,ETB; OUTPUT = W C---------------------------------------------------------------- CALL VERTVL(DTI2) C CALL BCOND(5) C DO 307 K=1,KM DO 307 J=1,JM DO 307 I=1,IM UF(I,J,K)=0.0 307 VF(I,J,K)=0.0 C---------------------------------------------------------------- C COMPUTE Q2F AN Q2LF USING UF AND VF AS TEMPORARY VARIABLES C---------------------------------------------------------------- CALL ADVQ(Q2B,Q2,DTI2,UF) CALL ADVQ(Q2LB,Q2L,DTI2,VF) CALL PROFQ(DTI2) CALL BCOND(6) DO 310 K=1,KM DO 310 J=1,JM DO 310 I=1,IM Q2(I,J,K)=Q2(I,J,K) 1 +.5*SMOTH*(UF(I,J,K)+Q2B(I,J,K)-2.*Q2(I,J,K)) Q2L(I,J,K)=Q2L(I,J,K) 1 +.5*SMOTH*(VF(I,J,K)+Q2LB(I,J,K)-2.*Q2L(I,J,K)) Q2B(I,J,K)=Q2(I,J,K) Q2(I,J,K)=UF(I,J,K) Q2LB(I,J,K)=Q2L(I,J,K) Q2L(I,J,K)=VF(I,J,K) 310 CONTINUE C---------------------------------------------------------------- C COMPUTE TF AN SF USING UF AND VF AS TEMPORARY VARIABLES C---------------------------------------------------------------- if (mode.eq.4) goto 360 c C In this calculation, S = constant. However, we calculate C S to demonstrate that the same subroutines can be used to C calculate multiple scalars. C REST=2.*DTI/(10.*86400.) C TCLIM=TCLIM+REST*(T-TCLIM) CALL ADVT(TB,T,TCLIM,DTI2,UF) C CALL ADVT(SB,S,SCLIM,DTI2,VF) CALL PROFT(UF,WTSURF,SWRAD,TSURF,2,DTI2) C CALL PROFT(VF,WSSURF,SWRAD,SSURF,1,DTI2) CALL BCOND(4) DO 355 K=1,KM-1 DO 355 J=2,JM-1 DO 355 I=2,IM-1 IF(DZ(I,J,K).NE.0) THEN T(I,J,K)=T(I,J,K)+.5*SMOTH*(DZF(I,J,K)*UF(I,J,K) 1 +DZB(I,J,K)*TB(I,J,K)-2.*DZ(I,J,K)*T(I,J,K))/DZ(I,J,K) C S(I,J,K)=S(I,J,K)+.5*SMOTH*(DZF(I,J,K)*VF(I,J,K) C 1 +DZB(I,J,K)*SB(I,J,K)-2.*DZ(I,J,K)*S(I,J,K))/DZ(I,J,K) ENDIF TB(I,J,K)=T(I,J,K) T(I,J,K)=UF(I,J,K) C SB(I,J,K)=S(I,J,K) C S(I,J,K)=VF(I,J,K) 355 CONTINUE C 360 CONTINUE C CALL DENS(S,T,RHO) C C---------------------------------------------------------------- C COMPUTE UF AND VF C---------------------------------------------------------------- CALL ADVU(DRHOX,ADVX,DTI2) CALL ADVV(DRHOY,ADVY,DTI2) CALL PROFU(DTI2) CALL PROFV(DTI2) CALL BCOND(3) C DO 369 J=1,JM DO 369 I=1,IM 369 TPS(I,J)=0.0 DO 370 K=1,KM-1 DO 370 J=1,JM DO 370 I=2,IM 370 TPS(I,J)=TPS(I,J)+(UF(I,J,K)+UB(I,J,K)-2.0*U(I,J,K)) 1 *(DZ(I,J,K)+DZ(I-1,J,K))*DUM(I,J,K) DO 371 J=1,JM DO 371 I=2,IM 371 TPS(I,J)=TPS(I,J)/(H(I,J)+ET(I,J)+H(I-1,J)+ET(I-1,J)) DO 372 K=1,KM-1 DO 372 J=1,JM DO 372 I=2,IM 372 U(I,J,K)=U(I,J,K)+.5*SMOTH*(UF(I,J,K)+UB(I,J,K)-2.*U(I,J,K) 1 -TPS(I,J)) DO 3721 J=1,JM DO 3721 I=1,IM 3721 TPS(I,J)=0.0 DO 374 K=1,KM-1 DO 374 J=2,JM DO 374 I=1,IM 374 TPS(I,J)=TPS(I,J)+(VF(I,J,K)+VB(I,J,K)-2.0*V(I,J,K)) 1 *(DZ(I,J,K)+DZ(I,J-1,K))*DVM(I,J,K) DO 375 J=2,JM DO 375 I=1,IM 375 TPS(I,J)=TPS(I,J)/(H(I,J)+ET(I,J)+H(I,J-1)+ET(I,J-1)) DO 376 K=1,KM-1 DO 376 J=2,JM DO 376 I=1,IM 376 V(I,J,K)=V(I,J,K)+.5*SMOTH*(VF(I,J,K)+VB(I,J,K)-2.0*V(I,J,K) 1 -TPS(I,J)) DO 377 K=1,KM DO 377 J=1,JM DO 377 I=1,IM UB(I,J,K)=U(I,J,K) U(I,J,K)=UF(I,J,K) VB(I,J,K)=V(I,J,K) 377 V(I,J,K)=VF(I,J,K) C 8200 CONTINUE C DO 8210 J=1,JM DO 8210 I=1,IM EGB(I,J)=EGF(I,J) ETB(I,J)=ET(I,J) ET(I,J)=ETF(I,J) DT(I,J)=H(I,J)+ET(I,J) UTB(I,J)=UTF(I,J) 8210 VTB(I,J)=VTF(I,J) DO 8220 K=1,KM DO 8220 J=1,JM DO 8220 I=1,IM DZB(I,J,K)=DZ(I,J,K) DZ(I,J,K)=DZF(I,J,K) 8220 CONTINUE C DO I=1,IM C ITEM(I)=1000.*VA(I,4) C ENDDO C IF(MOD(IINT,240).EQ.0) THEN C Write(6,'(30I4)') (ITEM(I),I=1,IM) C ENDIF C-------------------------------------------------------------------- C BEGIN PRINT SECTION C-------------------------------------------------------------------- IF(MOD(IINT,IPRINT).NE.0) GO TO 7000 9001 CONTINUE WRITE(6,'(/,'' TIME ='',F10.2,'' IINT ='',I8, 1 '' IEXT ='',I8,'' IPRINT ='',I5,//)') TIME,IINT,IEXT,IPRINT c CALL PRXZ(' Z' ,TIME, Z,IM,ISK,JM,10,20,KM,0.,DT) c CALL PRXZ('ZZ' ,TIME,ZZ,IM,ISK,JM,10,20,KM,0.,DT) c CALL PRXY(' WTSURF ',TIME,WTSURF,IM,ISK,JM,JSK,0.) c CALL PRXY(' SWRAD ',TIME,SWRAD ,IM,ISK,JM,JSK,0.) c CALL PRXY(' ADVUA ',TIME,ADVUA,IP,ISK,JM,JSK,0.0) c CALL PRXY(' ADVVA ',TIME,ADVVA,IP,ISK,JM,1,0.0) c CALL PRXY(' ADX2D ',TIME,ADX2D,IP,ISK,JM,1,0.0) c CALL PRXY(' ADY2D ',TIME,ADY2D,IP,ISK,JM,1,0.0) c CALL PRXY(' WUSURF ',TIME,WUSURF,IP,ISK,JM,1,0.E0) c CALL PRXY(' WVSURF ',TIME,WVSURF,IP,ISK,JM,1,0.E0) c CALL PRXY(' WUBOT ',TIME,WUBOT,IP,ISK,JM,JSK,0.E0) c CALL PRXY(' WVBOT ',TIME,WVBOT,IP,ISK,JM,JSK,0.E0) c CALL PRXY(' DRX2D ',TIME,DRX2D,IP,ISK,JM,JSK,0.E0) c CALL PRXY(' DRY2D ',TIME,DRY2D,IP,ISK,JM,JSK,0.E0) CALL PRXY(' FREE SURFACE ',TIME, ELB,IP,ISK,JM,JSK,0.E-3) CALL PRXY(' AVERAGE U COMP.',TIME,UAB,IP,ISK,JM,JSK,0.E-3) CALL PRXY(' AVERAGE V COMP.',TIME,VAB,IP,ISK,JM,JSK,0.E-3) CALL FINDPSI C write(6,'(/)') C write(6,'(1x,24I5)') ((INT(H (I,J)*1.E+0),I=1,IM),J=JM,1,-1) C write(6,'(/)') C write(6,'(1x,24I5)') ((INT(ELB(I,J)*1.E+3),I=1,IM),J=JM,1,-1) C write(6,'(/)') C write(6,'(1x,24I5)') ((INT(PSI(I,J)*1.E-5),I=1,IM),J=JM,1,-1) C write(6,'(/)') C I=12 C write(6,'(1x,37I4)') ((INT(TB (I,J,K)*1.E+1),J=1,JM),K=1,KM-1) C write(6,'(/)') IF(MODE.NE.2) THEN c CALL PRXZ('DRHOX' ,TIME,DRHOX,IM,ISK,JM,10,25,KM,0.,DT,ZZND) c CALL PRXZ('DRHOY' ,TIME,DRHOY,IM,ISK,JM,10,25,KM,0.,DT,ZZND) C CALL PRXYZ('AAM ',TIME,AAM,IP,ISK,JM,1,KM,0.E0) c CALL PRXYZ(' RHO ',TIME,RHO,IP,ISK,JM,1,KM,0.E0) c CALL PRXZ('SB' ,TIME,SB,IM,ISK,JM,10,25,KM,1.,DT,ZZND) c CALL PRXZ('RMEAN' ,TIME,RMEAN,IM,ISK,JM,10,25,KM,0.,DT,ZZND) c CALL PRXZ('RHO',TIME,RHO,IM,ISK,JM,10,25,KM,0.0,DT,ZZND) c CALL PRXZ('Q2' ,TIME,Q2,IM,ISK,JM,10,25,KM,0.0,DT) CALL PRXZ('KMT' ,TIME,KMT,IM,ISK,JM,10,25,KM,0.0,DT) C CALL GSIGTOZ(ZZ,TB,KB ,TLEV,-1500.,TPS,IM,JM,KB) C CALL PRXY(' T AT 1500 M',TIME,TLEV,IP,ISK,JM,JSK,0.E-3) CALL PRXZ('U' ,TIME,U,IM,ISK,JM,10,25,KM,0.0,DT) CALL PRXZ('V' ,TIME,V,IM,ISK,JM,10,25,KM,0.0,DT) C CALL PRXZ('TB',TIME,TB,IM,ISK,JM,10,25,KM,1.0,DT ) CALL PRXZ('T' ,TIME,T,IM,ISK,JM,10,25,KM,0.1,DT ) CALL PRYZ('T ',TIME,T,IM, 3,23,JM,JSK,KM,0.1,DT) c CALL PRXZ('TCLIM' ,TIME,TCLIM,IM,ISK,JM,10,25,KM,0.0,KB ) c CALL PRXZ('W' ,TIME,W,IM,ISK,JM,10,25,KM,0.0,KB ) c CALL PRXZ('KB ' ,TIME,FLOAT(KB ),IM,ISK,JM,10,25,1 ,0.0,DT) ENDIF C C WRITE(80) TIME,UB,VB,UAB,VAB,TB,SB,ELB,Q2B,Q2LB,KM,KH C VTOT=0. ATOT=0. TAVER=0. SAVER=0. EAVER=0. DO 8888 K=1,KM-1 DO 8888 J=1,JM DO 8888 I=1,IM DAREA=DX(I,J)*DY(I,J) DVOL =DAREA*DZ(I,J,K)*FSM(I,J,K) VTOT=VTOT+DVOL TAVER=TAVER+TB(I,J,K)*DVOL SAVER=SAVER+SB(I,J,K)*DVOL 8888 CONTINUE DO 8889 J=1,JM DO 8889 I=1,IM DAREA=DX(I,J)*DY(I,J)*FSM(I,J,1) ATOT=ATOT+DAREA EAVER=EAVER+ELB(I,J)*DAREA c EAVER=EAVER+ET(I,J)*DAREA 8889 CONTINUE TAVER=TAVER/VTOT SAVER=SAVER/VTOT EAVER=EAVER/ATOT WRITE(6,'('' VTOT ='',E20.8,'' ATOT ='',E20.8/, 1 '' EAVER ='',E20.8 )') VTOT,ATOT,EAVER WRITE(6,'('' TAVER ='',E20.8, 1 '' SAVER ='',E20.8/)') TAVER,SAVER IF(VAMAX.GT.100.) THEN CALL PRXY(' EL ',TIME, EL ,IP,ISK,JM,1,0.E0) CALL PRXY(' ELF ',TIME, ELF,IP,ISK,JM,1,0.E0) CALL PRXY(' UA ',TIME,UA ,IP,ISK,JM,1,0.E0) CALL PRXY(' VA ',TIME,VA ,IP,ISK,JM,1,0.E0) CALL PRXY(' UAF ',TIME,UAF,IP,ISK,JM,1,0.E0) CALL PRXY(' VAF ',TIME,VAF,IP,ISK,JM,1,0.E0) WRITE(6,'(''**************************************'')') WRITE(6,'(''********** ABNORMAL JOB END **********'')') WRITE(6,'(''*********** USER TERMINATED **********'')') WRITE(6,'('' VAMAX ='',E12.3)') VAMAX WRITE(6,'('' IINT,IEXT ='',2I6)') IINT,IEXT STOP C ENDIF C WRITE(10) TIME,ALON,ALAT,Z,ZZ,H,Q2,KMT,UB,VB,UAB,VAB,TB,SB,ELB,PSI C 7000 CONTINUE C-------------------------------------------------------------------- C END PRINT SECTION C-------------------------------------------------------------------- 9000 CONTINUE C*********************************************************************** C END NUMERICAL INTEGRATION * C*********************************************************************** C WRITE(140) TIME, 1 WUBOT,WVBOT,AAM2D,UA,UAB,VA,VAB,EL,ELB,ET,ETB,EGB,DZ,DZB,Z,ZZ, 2 UTB,VTB,U,UB,W,V,VB,T,TB,S,SB,RHO,ADX2D,ADY2D,ADVUA,ADVVA, 3 KMT,KH,KQ,L,GH,Q2,Q2B,AAM,Q2L,Q2LB C STOP END C SUBROUTINE ADVAVE(ADVUA,ADVVA,MODE) INCLUDE 'comblkg.h' DIMENSION ADVUA(IM,JM),ADVVA(IM,JM),CURV2D(IM,JM) EQUIVALENCE (TPS,CURV2D) C--------------------------------------------------------------------- C CALCULATE U-ADVECTION & DIFFUSION C--------------------------------------------------------------------- C C-------- ADVECTIVE FLUXES ------------------------------------------- C DO 200 J=1,JM DO 200 I=1,IM 200 ADVUA(I,J)=0. DO 300 J=2,JM DO 300 I=2,IMM1 300 FLUXUA(I,J)=.125*((D(I+1,J)+D(I,J))*UA(I+1,J) 1 +(D(I,J)+D(I-1,J))*UA(I,J)) 2 *(UA(I+1,J)+UA(I,J)) DO 400 J=2,JM DO 400 I=2,IM 400 FLUXVA(I,J)=.125*((D(I,J)+D(I,J-1))*VA(I,J) 1 +(D(I-1,J)+D(I-1,J-1))*VA(I-1,J)) 2 *(UA(I,J)+UA(I,J-1)) C----------- ADD VISCOUS FLUXES --------------------------------- DO 460 J=2,JM DO 460 I=2,IMM1 460 FLUXUA(I,J)=FLUXUA(I,J) 1 -D(I,J)*2.*AAM2D(I,J)*(UAB(I+1,J)-UAB(I,J))/DX(I,J) DO 470 J=2,JM DO 470 I=2,IM TPS(I,J)=.25*(D(I,J)+D(I-1,J)+D(I,J-1)+D(I-1,J-1)) 1 *(AAM2D(I,J)+AAM2D(I,J-1)+AAM2D(I-1,J)+AAM2D(I-1,J-1)) 2 *((UAB(I,J)-UAB(I,J-1)) 3 /(DY(I,J)+DY(I-1,J)+DY(I,J-1)+DY(I-1,J-1)) 4 +(VAB(I,J)-VAB(I-1,J)) 5 /(DX(I,J)+DX(I-1,J)+DX(I,J-1)+DX(I-1,J-1)) ) FLUXUA(I,J)=FLUXUA(I,J)*DY(I,J)*FSM(I,J,1) FLUXVA(I,J)=(FLUXVA(I,J)-TPS(I,J)) 1 *.25*(DX(I,J)+DX(I-1,J)+DX(I,J-1)+DX(I-1,J-1)) 2 *DUM(I,J,1)*DUM(I,J-1,1)*DVM(I,J,1)*DVM(I-1,J,1) 470 CONTINUE C---------------------------------------------------------------- DO 480 J=2,JMM1 DO 480 I=2,IMM1 480 ADVUA(I,J)=FLUXUA(I,J)-FLUXUA(I-1,J) 1 +FLUXVA(I,J+1)-FLUXVA(I,J) C---------------------------------------------------------------- C CALCULATE V-ADVECTION & DIFFUSION C---------------------------------------------------------------- C DO 481 J=1,JM DO 481 I=1,IM 481 ADVVA(I,J)=0. C---------ADVECTIVE FLUXES ---------------------------- DO 700 J=2,JM DO 700 I=2,IM 700 FLUXUA(I,J)=.125*((D(I,J)+D(I-1,J))*UA(I,J) 1 +(D(I,J-1)+D(I-1,J-1))*UA(I,J-1))* 2 (VA(I-1,J)+VA(I,J)) DO 800 J=2,JMM1 DO 800 I=2,IM 800 FLUXVA(I,J)=.125*((D(I,J+1)+D(I,J)) 1 *VA(I,J+1)+(D(I,J)+D(I,J-1))*VA(I,J)) 2 *(VA(I,J+1)+VA(I,J)) C------- ADD VISCOUS FLUXES ----------------------------------- DO 860 J=2,JMM1 DO 860 I=2,IM 860 FLUXVA(I,J)=FLUXVA(I,J) 1 -D(I,J)*2.*AAM2D(I,J)*(VAB(I,J+1)-VAB(I,J))/DY(I,J) DO 870 J=2,JM DO 870 I=2,IM FLUXVA(I,J)=FLUXVA(I,J)*DX(I,J)*FSM(I,J,1) 870 FLUXUA(I,J)=(FLUXUA(I,J)-TPS(I,J)) 1 *.25*(DY(I,J)+DY(I-1,J)+DY(I,J-1)+DY(I-1,J-1)) 2 *DUM(I,J,1)*DUM(I,J-1,1)*DVM(I,J,1)*DVM(I-1,J,1) C--------------------------------------------------------------- DO 880 J=2,JMM1 DO 880 I=2,IMM1 880 ADVVA(I,J)=FLUXUA(I+1,J)-FLUXUA(I,J) 1 +FLUXVA(I,J)-FLUXVA(I,J-1) C C--------------------------------------------------------------- IF(MODE.NE.2) GO TO 5000 DO 100 J=2,JMM1 DO 100 I=2,IMM1 WUBOT(I,J)=-0.5*(CBC(I,J)+CBC(I-1,J)) 1 *SQRT(UAB(I,J)**2+(.25*(VAB(I,J) 2 +VAB(I,J+1)+VAB(I-1,J)+VAB(I-1,J+1)))**2)*UAB(I,J) 100 CONTINUE DO 102 J=2,JMM1 DO 102 I=2,IMM1 WVBOT(I,J)=-0.5*(CBC(I,J)+CBC(I,J-1)) 1 *SQRT((.25*(UAB(I,J)+UAB(I+1,J) 2 +UAB(I,J-1)+UAB(I+1,J-1)))**2+VAB(I,J)**2)*VAB(I,J) 102 CONTINUE DO 120 J=2,JMM1 DO 120 I=2,IMM1 CURV2D(I,J)=.25*((VA(I,J+1)+VA(I,J))*(DY(I+1,J)-DY(I-1,J)) 1 -(UA(I+1,J)+UA(I,J))*(DX(I,J+1)-DX(I,J-1)) ) 2 /(DX(I,J)*DY(I,J)) ADVUA(I,J)=ADVUA(I,J) 1 -ARU(I,J)*.25*( CURV2D(I,J)*D(I,J)*(VA(I,J+1)+VA(I,J)) 2 +CURV2D(I-1,J)*D(I-1,J)*(VA(I-1,J+1)+VA(I-1,J)) ) ADVVA(I,J)=ADVVA(I,J) 1 +ARV(I,J)*.25*( CURV2D(I,J)*D(I,J)*(UA(I+1,J)+UA(I,J)) 2 +CURV2D(I,J-1)*D(I,J-1)*(UA(I+1,J-1)+UA(I,J-1)) ) 120 CONTINUE 5000 CONTINUE C RETURN END C SUBROUTINE ADVCT(ADVX,ADVY) C====================================================================== C This routine calculates the horizontal portions of momentum advection C well in advance of their use in ADVU and ADVV so that their vertical C integrals (created in MAIN) may be used in the external mode calculation. C====================================================================== INCLUDE 'comblkg.h' DIMENSION ADVX(IM,JM,KM),ADVY(IM,JM,KM), 1 XFLUX(IM,JM,KM),YFLUX(IM,JM,KM), 1 CURV(IM,JM,KM) EQUIVALENCE (A,XFLUX),(C,YFLUX),(EE,CURV) C DO 60 K=1,KM-1 DO 60 J=2,JMM1 DO 60 I=2,IMM1 60 CURV(I,J,K)= 1 +.25*((V(I,J+1,K)+V(I,J,K))*(DY(I+1,J)-DY(I-1,J)) 2 -(U(I+1,J,K)+U(I,J,K))*(DX(I,J+1)-DX(I,J-1)) ) 3 /(DX(I,J)*DY(I,J)) C----------------------------------------------------------------- C Calculate x-component of velocity advection C----------------------------------------------------------------- DO 59 K=1,KM DO 59 J=1,JM DO 59 I=1,IM ADVX(I,J,K)=0. XFLUX(I,J,K)=0. 59 YFLUX(I,J,K)=0. C C******** HORIZONTAL ADVECTION FLUXES ***************************** DO 100 K=1,KM-1 DO 100 J=1,JM DO 100 I=2,IMM1 100 XFLUX(I,J,K)=.125*((DZ(I+1,J,K)+DZ(I,J,K))* 1 U(I+1,J,K)+(DZ(I,J,K)+DZ(I-1,J,K)) 2 *U(I,J,K))*(U(I+1,J,K)+U(I,J,K)) DO 120 K=1,KM-1 DO 120 J=2,JM DO 120 I=2,IM 120 YFLUX(I,J,K)=.125*((DZ(I,J,K)+DZ(I,J-1,K)) 1 *V(I,J,K)+(DZ(I-1,J,K)+DZ(I-1,J-1,K)) 2 *V(I-1,J,K))*(U(I,J,K)+U(I,J-1,K)) C****** ADD HORIZONTAL DIFFUSION FLUXES **************************** DO 130 K=1,KM-1 DO 130 J=2,JM DO 130 I=2,IMM1 XFLUX(I,J,K)=XFLUX(I,J,K) 1 -DZ(I,J,K)*AAM(I,J,K)*2.*(UB(I+1,J,K)-UB(I,J,K))/DX(I,J) DTAAM=(AAM(I,J,K)+AAM(I-1,J,K)+AAM(I,J-1,K)+AAM(I-1,J-1,K)) 1 *.25*(DZ(I,J,K)+DZ(I-1,J,K)+DZ(I,J-1,K)+DZ(I-1,J-1,K)) YFLUX(I,J,K)=YFLUX(I,J,K) 1 -DTAAM*((UB(I,J,K)-UB(I,J-1,K)) 2 /(DY(I,J)+DY(I-1,J)+DY(I,J-1)+DY(I-1,J-1)) 3 +(VB(I,J,K)-VB(I-1,J,K)) 4 /(DX(I,J)+DX(I-1,J)+DX(I,J-1)+DX(I-1,J-1))) C XFLUX(I,J,K)=DY(I,J)*XFLUX(I,J,K)*FSM(I,J,K) YFLUX(I,J,K) 1 =.25*(DX(I,J)+DX(I-1,J)+DX(I,J-1)+DX(I-1,J-1))*YFLUX(I,J,K) 2 *DUM(I,J,K)*DUM(I,J-1,K)*DVM(I,J,K)*DVM(I-1,J,K) 130 CONTINUE C C IF(IINT.EQ.960) THEN C SCALA=0. C CALL PRXZ('XFLUX,Jslice',TIME,XFLUX,IM,1,JM,20,40,KM,SCALA,DT) C CALL PRYZ('XFLUX,Islice',TIME,XFLUX,IM,10,30,JM,1,KM,SCALA,DT) C ENDIF C C******** DO HORIZ. ADVECTION ******* DO 146 K=1,KM-1 DO 146 J=2,JMM1 DO 146 I=2,IMM1 146 ADVX(I,J,K)= 1 +XFLUX(I,J,K)-XFLUX(I-1,J,K) 2 +YFLUX(I,J+1,K)-YFLUX(I,J,K) 3 -ARU(I,J)*.25*(CURV(I,J,K)*DZ(I,J,K)*(V(I,J+1,K)+V(I,J,K)) 4 +CURV(I-1,J,K)*DZ(I-1,J,K)*(V(I-1,J+1,K)+V(I-1,J,K))) C C----------------------------------------------------------------- C Calculate y-component of velocity advection C----------------------------------------------------------------- DO 299 K=1,KM DO 299 J=1,JM DO 299 I=1,IM ADVY(I,J,K)=0. XFLUX(I,J,K)=0. 299 YFLUX(I,J,K)=0. C C********** HORIZONTAL ADVECTION FLUXES ************************** DO 300 K=1,KM-1 DO 300 J=2,JM DO 300 I=2,IM 300 XFLUX(I,J,K)=.125*((DZ(I,J,K)+DZ(I-1,J,K))*U(I,J,K) 1 +(DZ(I,J-1,K)+DZ(I-1,J-1,K))*U(I,J-1,K)) 2 *(V(I,J,K)+V(I-1,J,K)) DO 320 K=1,KM-1 DO 320 J=2,JMM1 DO 320 I=1,IM 320 YFLUX(I,J,K)=.125*((DZ(I,J+1,K)+DZ(I,J,K))*V(I,J+1,K) 1 +(DZ(I,J,K)+DZ(I,J-1,K))*V(I,J,K)) 2 *(V(I,J+1,K)+V(I,J,K)) C******* ADD HORIZONTAL DIFFUSION FLUXES ************************** DO 700 K=1,KM-1 DO 700 J=2,JMM1 DO 700 I=2,IM DTAAM=(AAM(I,J,K)+AAM(I-1,J,K)+AAM(I,J-1,K)+AAM(I-1,J-1,K)) 1 *.25*(DZ(I,J,K)+DZ(I-1,J,K)+DZ(I,J-1,K)+DZ(I-1,J-1,K)) XFLUX(I,J,K)=XFLUX(I,J,K) 1 -DTAAM*((UB(I,J,K)-UB(I,J-1,K)) 2 /(DY(I,J)+DY(I-1,J)+DY(I,J-1)+DY(I-1,J-1)) 3 +(VB(I,J,K)-VB(I-1,J,K)) 4 /(DX(I,J)+DX(I-1,J)+DX(I,J-1)+DX(I-1,J-1))) YFLUX(I,J,K)=YFLUX(I,J,K) 1 -DZ(I,J,K)*AAM(I,J,K)*2.*(VB(I,J+1,K)-VB(I,J,K))/DY(I,J) C XFLUX(I,J,K) 1 =.25*(DY(I,J)+DY(I-1,J)+DY(I,J-1)+DY(I-1,J-1))*XFLUX(I,J,K) 2 *DUM(I,J,K)*DUM(I,J-1,K)*DVM(I,J,K)*DVM(I-1,J,K) YFLUX(I,J,K)=DX(I,J)*YFLUX(I,J,K)*FSM(I,J,1) 700 CONTINUE C C********** DO HORIZ. ADVECTION ************ DO 400 K=1,KM-1 DO 400 J=2,JMM1 DO 400 I=2,IMM1 400 ADVY(I,J,K)= 1 +XFLUX(I+1,J,K)-XFLUX(I,J,K) 2 +YFLUX(I,J,K)-YFLUX(I,J-1,K) 3 +ARV(I,J)*.25*(CURV(I,J,K)*DZ(I,J,K)*(U(I+1,J,K)+U(I,J,K)) 4 +CURV(I,J-1,K)*DZ(I,J-1,K)*(U(I+1,J-1,K)+U(I,J-1,K))) C RETURN END SUBROUTINE ADVQ(QB,Q,DTI2,QF) INCLUDE 'comblkg.h' DIMENSION QB(IM,JM,KM),Q(IM,JM,KM),QF(IM,JM,KM) DIMENSION XFLUX(IM,JM,KM),YFLUX(IM,JM,KM) EQUIVALENCE (XFLUX,A),(YFLUX,C) C C******* HORIZONTAL ADVECTION ************************************ DO 110 K=2,KM-1 DO 110 J=2,JM DO 110 I=2,IM XFLUX(I,J,K)=.25*(Q(I,J,K)+Q(I-1,J,K)) 1 *(U(I,J,K)+U(I,J,K-1)) 110 YFLUX(I,J,K)=.25*(Q(I,J,K)+Q(I,J-1,K)) 1 *(V(I,J,K)+V(I,J,K-1)) C******* HORIZONTAL DIFFUSION ************************************ DO 315 K=2,KM-1 DO 315 J=2,JM DO 315 I=2,IM XFLUX(I,J,K)=XFLUX(I,J,K) 1 -.5*(AAM(I,J,K)+AAM(I-1,J,K)+AAM(I,J,K-1)+AAM(I-1,J,K-1)) 2 *(QB(I,J,K)-QB(I-1,J,K))*DUM(I,J,K) 3 /(DX(I,J)+DX(I-1,J)) YFLUX(I,J,K)=YFLUX(I,J,K) 1 -.5*(AAM(I,J,K)+AAM(I,J-1,K)+AAM(I,J,K-1)+AAM(I,J-1,K-1)) 2 *(QB(I,J,K)-QB(I,J-1,K))*DVM(I,J,K) 3 /(DY(I,J)+DY(I,J-1)) XFLUX(I,J,K)=.25*(DZZ(I,J,K-1)+DZZ(I-1,J,K-1)) 1 *(DY(I,J)+DY(I-1,J))*XFLUX(I,J,K) YFLUX(I,J,K)=.25*(DZZ(I,J,K-1)+DZZ(I,J-1,K-1)) 1 *(DX(I,J)+DX(I,J-1))*YFLUX(I,J,K) 315 CONTINUE C****** VERTICAL ADVECTION; ADD FLUX TERMS ;THEN STEP FORWARD IN TIME ** DO 230 K=2,KM-1 DO 230 J=2,JMM1 DO 230 I=2,IMM1 QF(I,J,K)=.5*(W(I,J,K-1)*Q(I,J,K-1)-W(I,J,K+1)*Q(I,J,K+1)) 1 *ART(I,J) 2 +XFLUX(I+1,J,K)-XFLUX(I,J,K) 3 +YFLUX(I,J+1,K)-YFLUX(I,J,K) 230 QF(I,J,K)=((DZB(I,J,K)+DZB(I,J,K-1))*ART(I,J)*QB(I,J,K) 1 -2.*DTI2*QF(I,J,K)) 2 /((DZF(I,J,K)+DZF(I,J,K-1))*ART(I,J)) C RETURN END C SUBROUTINE ADVT(FB,F,FCLIM,DTI2,FF) C C THIS SUBROUTINE INTEGRATES CONSERVATIVE SCALAR EQUATIONS C INCLUDE 'comblkg.h' DIMENSION FB(IM,JM,KM),F(IM,JM,KM),FF(IM,JM,KM),FCLIM(IM,JM,KM) DIMENSION XFLUX(IM,JM,KM),YFLUX(IM,JM,KM) EQUIVALENCE (XFLUX,A),(YFLUX,C) C DO 529 J=1,JM DO 529 I=1,IM kbb=kb(i,J) kbb1=kb(i,J)-1 F(I,J,KBB)=F(I,J,KBB1) 529 CONTINUE C C******* DO ADVECTION FLUXES ************************************** DO 530 K=1,KM-1 DO 530 J=2,JM DO 530 I=2,IM FF(I,J,K)=0. XFLUX(I,J,K)=.5*(F(I,J,K)+F(I-1,J,K))*U(I,J,K) 530 YFLUX(I,J,K)=.5*(F(I,J,K)+F(I,J-1,K))*V(I,J,K) C****** ADD DIFFUSIVE FLUXES ************************************* C FB=FB-FCLIM c DO 100 K=1,KM-1 DO 100 J=2,JM DO 100 I=2,IM XFLUX(I,J,K)=XFLUX(I,J,K) 1 -TPRNI*(AAM(I,J,K)+AAM(I-1,J,K)) 2 *(FB(I,J,K)-FB(I-1,J,K))*DUM(I,J,K)/(DX(I,J)+DX(I-1,J)) YFLUX(I,J,K)=YFLUX(I,J,K) 1 -TPRNI*(AAM(I,J,K)+AAM(I,J-1,K)) 2 *(FB(I,J,K)-FB(I,J-1,K))*DVM(I,J,K)/(DY(I,J)+DY(I,J-1)) XFLUX(I,J,K) 1 =.25*(DZ(I,J,K)+DZ(I-1,J,K))*(DY(I,J)+DY(I-1,J))*XFLUX(I,J,K) YFLUX(I,J,K) 1 =.25*(DZ(I,J,K)+DZ(I,J-1,K))*(DX(I,J)+DX(I,J-1))*YFLUX(I,J,K) 100 CONTINUE c FB=FB+FCLIM c C****** DO VERTICAL ADVECTION ************************************* DO 505 J=2,JMM1 DO 505 I=2,IMM1 505 FF(I,J,1)=-.5*(F(I,J,1)+F(I,J,2))*W(I,J,2)*ART(I,J) DO 520 K=2,KM-1 DO 520 J=2,JMM1 DO 520 I=2,IMM1 520 FF(I,J,K)=.5*((F(I,J,K-1)+F(I,J,K))*W(I,J,K) 1 -(F(I,J,K)+F(I,J,K+1))*W(I,J,K+1))*ART(I,J) C****** ADD NET HORIZONTAL FLUXES; THEN STEP FORWARD IN TIME ********** DO 120 K=1,KM-1 DO 120 J=2,JMM1 DO 120 I=2,IMM1 FF(I,J,K)=FF(I,J,K) 1 +XFLUX(I+1,J,K)-XFLUX(I,J,K) 2 +YFLUX(I,J+1,K)-YFLUX(I,J,K) FF(I,J,K)=(FB(I,J,K)*DZB(I,J,K)*ART(I,J)-DTI2*FF(I,J,K)) 1 /(DZF(I,J,K)*ART(I,J)) 120 CONTINUE RETURN END SUBROUTINE ADVU(DRHOX,ADVX,DTI2) INCLUDE 'comblkg.h' DIMENSION DRHOX(IM,JM,KM),ADVX(IM,JM,KM) C C Do vertical advection DO 100 K=1,KM DO 100 J=1,JM DO 100 I=1,IM 100 UF(I,J,K)=0. DO 140 K=2,KM-1 DO 140 J=1,JM DO 140 I=2,IM 140 UF(I,J,K)=.25*(W(I,J,K)+W(I-1,J,K))*(U(I,J,K)+U(I,J,K-1)) C****COMBINE HOR. and VERT. ADVECTION with C -FVD + GDEG/DX + BAROCLINIC TERM ********************** DO 150 K=1,KM-1 DO 150 J=2,JMM1 DO 150 I=2,IMM1 CORFAC=1. 1 /(DVM(I,J+1,K)+DVM(I,J,K)+DVM(I-1,J+1,K)+DVM(I-1,J,K)+1.E-10) CORFAC=.25 150 UF(I,J,K)=ADVX(I,J,K)+(UF(I,J,K)-UF(I,J,K+1))*ARU(I,J) 1 -ARU(I,J)*CORFAC*(COR(I,J)*DZ(I,J,K)*(V(I,J+1,K)+V(I,J,K)) 2 +COR(I-1,J)*DZ(I-1,J,K)*(V(I-1,J+1,K)+V(I-1,J,K))) 3 +GRAV*.25*(DZ(I,J,K)+DZ(I-1,J,K)) 4 *.5*(EGF(I,J)-EGF(I-1,J)+EGB(I,J)-EGB(I-1,J)) 5 *(DY(I,J)+DY(I-1,J)) 6 +DRHOX(I,J,K) C******* STEP FORWARD IN TIME *********************************** DO 190 K=1,KM-1 DO 190 J=2,JMM1 DO 190 I=2,IMM1 190 UF(I,J,K)= 1 ((DZB(I,J,K)+DZB(I-1,J,K))*ARU(I,J)*UB(I,J,K) 2 -2.*DTI2*UF(I,J,K)) 3 /((DZF(I,J,K)+DZF(I-1,J,K))*ARU(I,J)) RETURN END SUBROUTINE ADVV(DRHOY,ADVY,DTI2) INCLUDE 'comblkg.h' DIMENSION DRHOY(IM,JM,KM),ADVY(IM,JM,KM) C C Do vertical advection DO 100 K=1,KM DO 100 J=1,JM DO 100 I=1,IM 100 VF(I,J,K)=0. DO 140 K=2,KM-1 DO 140 J=2,JM DO 140 I=1,IM 140 VF(I,J,K)=.25*(W(I,J,K)+W(I,J-1,K))*(V(I,J,K)+V(I,J,K-1)) C****COMBINE HOR. and VERT. ADVECTION with C +FUD + GDEG/DY + BAROCLINIC TERM ********************** DO 340 K=1,KM-1 DO 340 J=2,JMM1 DO 340 I=2,IMM1 CORFAC=1. 1 /(DUM(I+1,J,K)+DUM(I,J,K)+DUM(I+1,J-1,K)+DUM(I,J-1,K)+1.E-10) CORFAC=.25 340 VF(I,J,K)=ADVY(I,J,K)+(VF(I,J,K)-VF(I,J,K+1))*ARV(I,J) 1 +ARV(I,J)*CORFAC*(COR(I,J)*DZ(I,J,K)*(U(I+1,J,K)+U(I,J,K)) 2 +COR(I,J-1)*DZ(I,J-1,K)*(U(I+1,J-1,K)+U(I,J-1,K))) 3 +GRAV*.25*(DZ(I,J,K)+DZ(I,J-1,K)) 4 *.5*(EGF(I,J)-EGF(I,J-1)+EGB(I,J)-EGB(I,J-1)) 5 *(DX(I,J)+DX(I,J-1)) 6 +DRHOY(I,J,K) C******* STEP FORWARD IN TIME *************************************** DO 390 K=1,KM-1 DO 390 J=2,JMM1 DO 390 I=2,IMM1 390 VF(I,J,K)= 1 ((DZB(I,J,K)+DZB(I,J-1,K))*ARV(I,J)*VB(I,J,K) 2 -2.*DTI2*VF(I,J,K)) 3 /((DZF(I,J,K)+DZF(I,J-1,K))*ARV(I,J)) RETURN END C SUBROUTINE BAROPG(DRHOX,DRHOY) INCLUDE 'comblkg.h' DIMENSION DRHOX(IM,JM,KM),DRHOY(IM,JM,KM) C DO 299 K=1,KM DO 299 J=1,JM DO 299 I=1,IM 299 RHO(I,J,K)=RHO(I,J,K)-RMEAN(I,J,K) C C X COMPONENT OF BAROCLINIC PRESSURE GRADIENT C DO 300 J=2,JM DO 300 I=2,IM 300 DRHOX(I,J,1)=.5*GRAV*(-ZZ(I,J,1)-ZZ(I-1,J,1)) 1 *(RHO(I,J,1)-RHO(I-1,J,1)) DO 310 K=2,KM-1 DO 310 J=2,JM DO 310 I=2,IM 310 DRHOX(I,J,K)=DRHOX(I,J,K-1) 1 +GRAV*.25*(ZZ(I,J,K-1)-ZZ(I,J,K)+ZZ(I-1,J,K-1)-ZZ(I-1,J,K)) 2 *(RHO(I,J,K)-RHO(I-1,J,K)+RHO(I,J,K-1)-RHO(I-1,J,K-1)) 3 -GRAV*.25*(ZZ(I,J,K)-ZZ(I-1,J,K)+ZZ(I,J,K-1)-ZZ(I-1,J,K-1)) 4 *(RHO(I,J,K-1)+RHO(I-1,J,K-1)-RHO(I,J,K)-RHO(I-1,J,K)) C DO 360 K=1,KM-1 DO 360 J=2,JM DO 360 I=2,IM 360 DRHOX(I,J,K)=.25*(DZ(I,J,K)+DZ(I-1,J,K))*DRHOX(I,J,K) 1 *(DY(I,J)+DY(I-1,J))*DUM(I,J,K)*RAMP C C Y COMPONENT OF BAROCLINIC PRESSURE GRADIENT C DO 500 J=2,JM DO 500 I=2,IM 500 DRHOY(I,J,1)=.5*GRAV*(-ZZ(I,J,1)-ZZ(I,J-1,1)) 1 *(RHO(I,J,1)-RHO(I,J-1,1)) DO 510 K=2,KM-1 DO 510 J=2,JM DO 510 I=2,IM DRHOY(I,J,K)=DRHOY(I,J,K-1) 1 +GRAV*.25*(ZZ(I,J,K-1)-ZZ(I,J,K)+ZZ(I,J-1,K-1)-ZZ(I,J-1,K)) 2 *(RHO(I,J,K)-RHO(I,J-1,K)+RHO(I,J,K-1)-RHO(I,J-1,K-1)) 3 -GRAV*.25*(ZZ(I,J,K)-ZZ(I,J-1,K)+ZZ(I,J,K-1)-ZZ(I,J-1,K-1)) 4 *(RHO(I,J,K-1)+RHO(I,J-1,K-1)-RHO(I,J,K)-RHO(I,J-1,K)) 510 CONTINUE C DO 560 K=1,KM-1 DO 560 J=2,JM DO 560 I=2,IM 560 DRHOY(I,J,K)=.25*(DZ(I,J,K)+DZ(I,J-1,K))*DRHOY(I,J,K) 1 *(DX(I,J)+DX(I,J-1))*DVM(I,J,K)*RAMP C DO 571 K=1,KM DO 571 J=1,JM DO 571 I=1,IM 571 RHO(I,J,K)=RHO(I,J,K)+RMEAN(I,J,K) C RETURN END SUBROUTINE BCOND(IDX) INCLUDE 'comblkg.h' C C Closed boundary conditions are automatically enabled through C specification of the masks, DUM, DVM and FSM. For problems with C open boundaries, dependent variables must be specified on the C boundary grid points. C C C SPECIFICATION OF OPEN BOUNDARY CONDITIONS C BRACKETS OF * INDICATE INTERIOR (NON-BOUNDARY) GRID POINTS. C HORIZONTAL LOCATIONS OF EL, T AND S ARE COINCIDENT. C C I-1 I I+1 C C C U(JM) EL(JM) =. U(JM) BC C . C V(JM) . BC C . C *U(JMM1)* *EL(JMM1)* *U(JMM1)* C C *V(JMM1)* C . *V(IMM1)* V(IM) C . C . *U(IMM1)* *EL(IMM1)* U(IM) EL(IM) C . ". . . . . . . ." C . C . *V(IMM1)* V(IM) C . C . BC BC C C *V(3)* C C *U(2)* *EL(2)* *U(2)* C . C V(2) . BC C . C U(1) EL(1)=. U(1) BC C -------------------------------------------------------------------- C V(1) IS NOT USED C C DATA PI/3.14167/,GEE/9.807/ C GO TO (10,20,30,40,50,60), IDX C C----------------------------------------------------------------------- C EXTERNAL B.C.'S C In this example the governing boundary conditions are a radition C condition on UAF(IM,J) in the east and an inflow UAF(2,J) in the C west. The tangential velocities are set to zero on both boundaries. C These are only one set of possibilities and may not represent a choice C which yields the most physically realistic result. C----------------------------------------------------------------------- C C 10 CONTINUE C---------ELEVATION -------------------- C In this application, elevation is not a primary B.C. DO 110 J=1,JM ELF(1,J)=ELF(2,J) ELF(IM,J)=ELF(IMM1,J) 110 CONTINUE C DO 111 J=1,JM DO 111 I=1,IM 111 ELF(I,J)=ELF(I,J)*FSM(I,J,1) RETURN C C 20 CONTINUE C---------- VELOCITY -------------- DO 210 J=2,JMM1 C--- Governing B.C.s ------------- c UAF(2,J)=RAMP*UABW(J) c 1 *(H(2,J)+H(1,J))/(H(2,J)+ELF(2,J)+H(1,J)+ELF(1,J)) c UAF(IM,J)=RAMP*UABE(J)+SQRT(GRAV/H(IMM1,J))*(EL(IMM1,J)-ELE(J)) VAF(1,J)=0.0 VAF(IM,J)=0.0 C--------------------------- c UAF(1,J)=UAF(2,J) 210 CONTINUE C DO 131 J=1,JM DO 131 I=1,IM UAF(I,J)=UAF(I,J)*DUM(I,J,1) 131 VAF(I,J)=VAF(I,J)*DVM(I,J,1) RETURN C C 30 CONTINUE C----------------------------------------------------------------------- C INTERNAL VEL B.C.'S C----------------------------------------------------------------------- C Eastern and western radiation boundary conditions. Some smoothing C is used in the J-direction; it may not be necessary. C c DO 140 K=1,KM-1 c DO 140 J=2,JMM1 c GA=SQRT(H(IM,J)/HMAX) c UF(IM,J,K) c 1 =GA*(.25*U(IMM1,J-1,K)+.5*U(IMM1,J,K)+.25*U(IMM1,J+1,K)) c 2 +(1.-GA)*(.25*U(IM,J-1,K)+.5*U(IM,J,K)+.25*U(IM,J+1,K)) c GA=SQRT(H(1,J)/HMAX) c UF(2,J,K) c 1 =GA*(.25*U(3,J-1,K)+.5*U(3,J,K)+.25*U(3,J+1,K)) c 2 +(1.-GA)*(.25*U(2,J-1,K)+.5*U(2,J,K)+.25*U(2,J+1,K)) c 140 CONTINUE C DO 141 K=1,KM DO 141 J=1,JM VF(1,J,K)=0. VF(IM,J,K)=0. 141 CONTINUE C DO 160 K=1,KM-1 DO 160 J=1,JM DO 160 I=1,IM UF(I,J,K)=UF(I,J,K)*DUM(I,J,K) VF(I,J,K)=VF(I,J,K)*DVM(I,J,K) 160 CONTINUE C RETURN C C 40 CONTINUE C----------------------------------------------------------------------- C TEMP(UF) & SAL(VF) B.C.'S C----------------------------------------------------------------------- DO 220 K=1,KM-1 DO 220 J=1,JM U1=U(2,J,K)+ABS(U(2,J,K)) U2=U(2,J,K)-ABS(U(2,J,K)) UF(1,J,K)=T(1,J,K)-DTI/(DX(1,J)+DX(2,J)) 1 *( U2*(T(2,J,K)-T(1,J,K))+U1*(T(1,J,K)-TBW(J,K)) ) VF(1,J,K)=S(1,J,K)-DTI/(DX(1,J)+DX(2,J)) 1 *( U2*(S(2,J,K)-S(1,J,K))+U1*(S(1,J,K)-SBW(J,K)) ) C U1=U(IM,J,K)+ABS(U(IM,J,K)) U2=U(IM,J,K)-ABS(U(IM,J,K)) UF(IM,J,K)=T(IM,J,K)-DTI/(DX(IM,J)+DX(IMM1,J)) 1 *( U1*(T(IM,J,K)-T(IMM1,J,K))+U2*(TBE(J,K) -T(IM,J,K)) ) VF(IM,J,K)=S(IM,J,K)-DTI/(DX(IM,J)+DX(IMM1,J)) 1 *( U1*(S(IM,J,K)-S(IMM1,J,K))+U2*(SBE(J,K) -S(IM,J,K)) ) 220 CONTINUE C DO 240 K=1,KM-1 DO 240 J=1,JM DO 240 I=1,IM UF(I,J,K)=UF(I,J,K)*FSM(I,J,K) VF(I,J,K)=VF(I,J,K)*FSM(I,J,K) 240 CONTINUE RETURN C C 50 CONTINUE C---------------VERTICAL VEL. B. C.'S -------------------------------- DO 250 J=1,JM DO 250 I=1,IM KBB=KB(I,j)-1 DO 250 K=1,KBB W(I,J,K)=W(I,J,K)*FSM(I,J,K) 250 CONTINUE C RETURN C C 60 CONTINUE C---------------- Q2 AND Q2L B.C.'S ----------------------------------- DO 300 K=1,KM C DO 295 J=1,JM UF(IM,J,K)=1.E-10 VF(IM,J,K)=1.E-10 UF(1,J,K)=1.E-10 VF(1,J,K)=1.E-10 295 CONTINUE DO 297 J=1,JM DO 297 I=1,IM UF(I,J,K)=UF(I,J,K)*FSM(I,J,K) 297 VF(I,J,K)=VF(I,J,K)*FSM(I,J,K) 300 CONTINUE RETURN END SUBROUTINE DENS(SI,TI,RHOO) INCLUDE 'comblkg.h' DIMENSION SI(IM,JM,KM),TI(IM,JM,KM),RHOO(IM,JM,KM) C If using 32 bit precision, it is recommended that C TR, SR, P, RHOR , CR be made double precision. C and the E's in the constants should be changed C to D's. C C THIS SUBROUTINE COMPUTES DENSITY- 1.000 C T = POTENTIAL TEMPERATURE C ( See: Mellor, 1991, J. Atmos. Oceanic Tech., 609-611) C DO 1 K=1,KM-1 DO 1 J=1,JM DO 1 I=1,IM TR=TI(I,J,K)+TBIAS SR=SI(I,J,K)+SBIAS TR2=TR*TR TR3=TR2*TR TR4=TR3*TR C Approximate pressure in units of bars P=-GRAV*1.025*ZZ(I,J,K)*0.01 C RHOR = 0.157406 + 6.793952E-2*TR $ - 9.095290E-3*TR2 + 1.001685E-4*TR3 $ - 1.120083E-6*TR4 + 6.536332E-9*TR4*TR C RHOR = RHOR + (0.824493 - 4.0899E-3*TR $ + 7.6438E-5*TR2 - 8.2467E-7*TR3 $ + 5.3875E-9*TR4) * SR $ + (-5.72466E-3 + 1.0227E-4*TR $ - 1.6546E-6*TR2) * ABS(SR)**1.5 $ + 4.8314E-4 * SR*SR C CR=1449.1+.0821*P+4.55*TR-.045*TR2 $ +1.34*(SR-35.) RHOR=RHOR + 1.E5*P/(CR*CR)*(1.-2.*P/(CR*CR)) C RHOO(I,J,K)=RHOR*1.E-3*FSM(I,J,K) 1 CONTINUE C RETURN END C SUBROUTINE FINDPSI INCLUDE 'comblkg.h' DO 9004 J=1,JM DO 9004 I=1,IM 9004 PSI(I,J)=0.0 DO 9005 J=2,JM DO 9005 I=1,IMM1 9005 PSI(I+1,J)=PSI(I,J)-VAB(I,J)*.25*(D(I,J)+D(I,J-1)) 1 *(DX(I,J)+DX(I,J-1)) CALL PRXY(' PSI FROM V ',TIME,PSI,IM,ISK,JM,JSK,1.E6) DO 9006 J=2,JMM1 DO 9006 I=2,IM 9006 PSI(I,J+1)=PSI(I,J)+.25*UAB(I,J)*(D(I,J)+D(I-1,J)) 1 *(DY(I,J)+DY(I-1,J)) CALL PRXY(' PSI FROM U ',TIME,PSI,IM,ISK,JM,JSK,1.E6) RETURN END C SUBROUTINE MAKZLV(INITLZ) include 'comblkg.h' DIMENSION KDZ(12) DATA KDZ/1,1,2,4,8,16,32,64,128,256,512,1024/ C *** C THIS SUBROUTINE ESTABLISHES THE VERTICAL Z-LEVEL GRID WITH LOG C DISTRIBUTIONS AT THE TOP AND A LINEAR DISTRIBUTION C BELOW. CHOOSE KL1. THE VALUE, KL1 = 2, C YIELDS AN EQUALLY SPACED DISTRIBUTION WITH NO LOG PORTION. C *** if(INITLZ.eq.0) then KL1=8 C LZ(1)=0. DO K=2,KL1 LZ(K)=LZ(K-1)+KDZ(K-1) ENDDO DELZ=LZ(KL1)-LZ(KL1-1) DO K=KL1+1,KM LZ(K)=LZ(K-1)+DELZ ENDDO C Normalize LZ. DO K=1,KM LZ(K)=-LZ(K)/LZ(KM) ENDDO DO K=1,KM-1 LZZ(K)=0.5*(LZ(K)+LZ(K+1)) ENDDO LZZ(KM)=2.*LZZ(KM-1)-LZZ(KM-2) WRITE(6,'(/,2X,''K'',7X,''LZ'',8X,''LZZ'',/)') WRITE(6,'('' '',I5,2F10.3)') (K,LZ(K),LZZ(K),K=1,KM) C c DO I=1,IM DO J=1,JM KB(I,J)=KM fsm(i,j,1)=0. DO K=2,KM fsm(i,j,k)=0. if(h(i,j).ge.-lz(k)*hmax) then FSM(I,J,K-1)=1. KB(I,J)=k endif ENDDO C PROFS need at least 2 levels to work. if(KB(I,J).EQ.2) THEN KB(I,J)=3 FSM(I,J,2)=1. endif ENDDO ENDDO C Reset H to constant levels and restore basin volume VOL1=0 DO I=1,IM DO J=1,JM AREA=DX(I,J)*DY(I,J)*FSM(I,J,1) VOL1=VOL1+H(I,J)*AREA ENDDO ENDDO VOL2=0. DO I=1,IM DO J=1,JM kbot=kb(i,j) h(i,j)=-fsm(i,j,1)*lz(kbot)*HMAX h(i,j)=amax1(h(i,j),1.) AREA=DX(I,J)*DY(I,J)*FSM(I,J,1) VOL2=VOL2+H(I,J)*AREA ENDDO ENDDO DO I=1,IM DO J=1,JM H(I,J)=H(I,J)*VOL1/VOL2 KBB=KB(I,J) IF(H(I,J).LT.100.) THEN H(I,J)=1. FSM(I,J,:)=0. KB(I,J)=KM ENDIF DO K=1,KM Z(I,J,K)=-LZ(K)*(H(I,J)+ET(I,J))/LZ(KBB) ENDDO DO K=1,KM-1 DZB(I,J,K)=(Z(I,J,K)-Z(I,J,K+1)) ENDDO ENDDO ENDDO C RETURN endif c DO J=1,JM DO I=1,IM KBB=KB(I,J) DO K=1,KM Z(I,J,K)=-LZ(K)*(H(I,J)+ET(I,J))/LZ(KBB) ZF(I,J,K)=-LZ(K)*(H(I,J)+ETF(I,J))/LZ(KBB) ZZ(I,J,K)=-LZZ(K)*(H(I,J)+ET(I,J))/LZ(KBB) ZZF(I,J,K)=-LZZ(K)*(H(I,J)+ETF(I,J))/LZ(KBB) ENDDO ENDDO ENDDO DO K=1,KM-1 DO I=1,IM DO J=1,JM DZZ(I,J,K)=(ZZ(I,J,K)-ZZ(I,J,K+1)) DZF(I,J,K)=(ZF(I,J,K)-ZF(I,J,K+1)) DZZF(I,J,K)=(ZZF(I,J,K)-ZZF(I,J,K+1)) ENDDO ENDDO ENDDO RETURN END C SUBROUTINE MAKSIG(INITLZ) include 'comblkg.h' DIMENSION KDZ(12) DATA KDZ/1,1,2,4,8,16,32,64,128,256,512,1024/ C *** C THIS SUBROUTINE ESTABLISHES THE VERTICAL SIGMA GRID WITH LOG C DISTRIBUTIONS AT THE TOP AND BOTTOM AND A LINEAR DISTRIBUTION C BETWEEN. CHOOSE KL1 AND KL2. THE VALUES, KL1 = 2 AND KL2 = KM-2, C YIELD AN EQUALLY SPACED DISTRIBUTION WITH NO LOG PORTIONS. C *** if(INITLZ.eq.0) then KL1=8 KL2=KM-2 C LZ(1)=0. DO K=2,KL1 LZ(K)=LZ(K-1)+KDZ(K-1) ENDDO DELZ=LZ(KL1)-LZ(KL1-1) DO K=KL1+1,KL2 LZ(K)=LZ(K-1)+DELZ ENDDO DO K=KL2+1,KM LZ(K)=LZ(K-1)+FLOAT(KDZ(KM-K+1))*DELZ/FLOAT(KDZ(KM-KL2)) ENDDO C Normalize LZ. DO K=1,KM LZ(K)=-LZ(K)/LZ(KM) ENDDO DO K=1,KM-1 LZZ(K)=0.5*(LZ(K)+LZ(K+1)) ENDDO LZZ(KM)=2.*LZZ(KM-1)-LZZ(KM-2) WRITE(6,'(/,2X,''K'',7X,''LZ'',8X,''LZZ'',/)') WRITE(6,'('' '',I5,2F10.3)') (K,LZ(K),LZZ(K),K=1,KM) DO K=1,KM DO I=1,IM DO J=1,JM fsm(i,j,k)=0. Z(I,J,K)=LZ(K)*H(I,J) ZZ(I,J,K)=LZZ(K)*H(I,J) ENDDO ENDDO ENDDO DO I=1,IM DO J=1,JM if(H(I,J).LT.1.0) H(I,J)=1. ENDDO ENDDO DO I=1,IM DO J=1,JM KB(I,J)=KM DO K=1,KM-1 DZB(I,J,K)=Z(I,J,K)-Z(I,J,K+1) if(h(i,j).gt.1.) FSM(I,J,K)=1. ENDDO ENDDO ENDDO RETURN endif DO I=1,IM DO J=1,JM DO K=1,KM Z(I,J,K)=LZ(K)*(H(I,J)+ET(I,J)) ZF(I,J,K)=LZ(K)*(H(I,J)+ETF(I,J)) ZZ(I,J,K)=LZZ(K)*(H(I,J)+ET(I,J)) ZZF(I,J,K)=LZZ(K)*(H(I,J)+ETF(I,J)) ENDDO ENDDO ENDDO DO K=1,KM-1 DO I=1,IM DO J=1,JM DZZ(I,J,K)=(ZZ(I,J,K)-ZZ(I,J,K+1)) DZF(I,J,K)=(ZF(I,J,K)-ZF(I,J,K+1)) DZZF(I,J,K)=(ZZF(I,J,K)-ZZF(I,J,K+1)) ENDDO ENDDO ENDDO RETURN END C SUBROUTINE MAKZPS(INITLZ) include 'comblkg.h' C23456 | | | | | | | xxxxxxx| DIMENSION KDZ(31) DATA KDZ /1,1,2,4,8,16,32,64,128,256,512,768,1014,1270,1526, 1 1782,2038,2294,2550,2806,3062,3318,3574,3830,4086,4342, 2 4596,4854,5110,5366,5622/ IF(INITLZ.EQ.0) then KL1=8 BOTSW=0.9 BOTDZ=0.05 write(6,'(/,''HMAX ='',F10.1)') HMAX write(6,'(''BOTSW ='',F10.3,'' BOTDZ ='',F10.3)')BOTSW,BOTDZ DO I=1,IM DO J=1,JM if(H(I,J).LT.1.0) H(I,J)=1. ENDDO ENDDO C CALL PRXZ('DZ',TIME,DZ,IM,1,JM,10,73,KM,0.0,DT) DO 50 I=1,IM DO 50 J=1,JM Z(I,J,1)=0. DO K=2,KL1 Z(I,J,K)=Z(I,J,K-1)+KDZ(K-1) ENDDO DELZ=Z(I,J,KL1)-Z(I,J,KL1-1) DO K=KL1+1,KM Z(I,J,K)=Z(I,J,K-1)+DELZ ENDDO DO K=1,KM Z(I,J,K)=-Z(I,J,K)*HMAX/Z(I,J,KM) ENDDO 50 CONTINUE DO 100 I=1,IM DO 100 J=1,JM DO K=1,KM IF(Z(I,J,K).LT.(-BOTSW*H(I,J))) GOTO 10 KE=K ENDDO 10 CONTINUE C Search to find starting bottom KDZ so that Z(KM-1)-Z(KM) C is not smaller than BOTDZ DO KGS=2,19 FAC=(H(I,J)+Z(I,J,KE))/FLOAT(KDZ(KM-KE+KGS)-KDZ(KGS)) IF(FAC*FLOAT(KDZ(KGS)).GT.BOTDZ) GOTO 20 ENDDO 20 CONTINUE C Z(I,J,KM)=-H(I,J) DO K=1,KM-KE IF((K+KGS).LE.21) THEN Z(I,J,KM-K)=-H(I,J)+FAC*FLOAT(KDZ(K+KGS)) ENDIF ENDDO DO K=1,KM IF(H(I,J).LE.1.0) Z(I,J,K)=-H(I,J)*FLOAT(K-1)/FLOAT(KM-1) ENDDO Z(I,J,KE)=0.5*(Z(I,J,KE+1)+Z(I,J,KE-1)) C write(6,'('' I,KE,KST ='',3I5)') I,KE,KST IF(J.EQ.10) write(6,'(''I,KE,KGS,FAC='',3I5,F10.3)') 1 I,KE,KGS,FAC 100 CONTINUE DO 200 I=1,IM DO 200 J=1,JM DO 200 K=1,KM fsm(i,j,k)=0. LLZ(I,J,K)=-Z(I,J,K)/Z(I,J,KM) Z(I,J,K)=LLZ(I,J,K)*H(I,J) 200 CONTINUE C DO K=1,KM-1 DO I=1,IM DO J=1,JM ZZ(I,J,K)=0.5*(LLZ(I,J,K)+LLZ(I,J,K+1))*H(I,J) ENDDO ENDDO ENDDO CALL PRXZ(' Z' ,TIME, Z,IM, 1,JM,10,20,KM,0.,DT) DO I=1,IM DO J=1,JM KB(I,J)=KM DO K=1,KM-1 DZB(I,J,K)=Z(I,J,K)-Z(I,J,K+1) if(h(i,j).gt.1.) fsm(I,J,K)=1. ENDDO TPS(I,J)=DZB(I,J,KB(I,J)-1) ENDDO ENDDO C RETURN endif DO I=1,IM DO J=1,JM DO K=1,KM Z(I,J,K)=LLZ(I,J,K)*(H(I,J)+ET(I,J)) ZF(I,J,K)=LLZ(I,J,K)*(H(I,J)+ETF(I,J)) ENDDO ENDDO ENDDO DO I=1,IM DO J=1,JM DO K=1,KM-1 ZZ(I,J,K)=0.5*(Z(I,J,K)+Z(I,J,K+1)) ZZF(I,J,K)=0.5*(ZF(I,J,K)+ZF(I,J,K+1)) ENDDO ENDDO ENDDO DO K=1,KM-1 DO I=1,IM DO J=1,JM DZZ(I,J,K)=(ZZ(I,J,K)-ZZ(I,J,K+1)) DZF(I,J,K)=(ZF(I,J,K)-ZF(I,J,K+1)) DZZF(I,J,K)=(ZZF(I,J,K)-ZZF(I,J,K+1)) ENDDO ENDDO ENDDO RETURN END C SUBROUTINE PROFQ(DT2) INCLUDE 'comblkg.h' REAL KN,KAPPA DIMENSION SM(IM,JM,KM),SH(IM,JM,KM) DIMENSION PROD(IM,JM,KM),KN(IM,JM,KM),BOYGR(IM,JM,KM) DIMENSION DH(IM,JM),CC(IM,JM,KM) EQUIVALENCE (A,SM),(C,SH),(PROD,KN),(TPS,DH) EQUIVALENCE (DTEF,CC) DATA A1,B1,A2,B2,C1/0.92,16.6,0.74,10.1,0.08/ DATA E1/1.8/,E2/1.33/,E3/1.0/ DATA KAPPA/0.40/,SQ/0.20/,SEF/1./ DATA GEE/9.806/ DATA SMALL/1.E-8/ C DO K=1,KM-1 DO J=1,JM DO I=1,IM A(I,J,K)=0. C(I,J,K)=0. UF(I,J,K)=UF(I,J,K)*FSM(I,J,K) VF(I,J,K)=VF(I,J,K)*FSM(I,J,K) ENDDO ENDDO ENDDO DO 100 K=2,KM-1 DO 100 J=1,JM DO 100 I=1,IM A(I,J,K)=-DT2*(KQ(I,J,K+1)+KQ(I,J,K)+2.*UMOL)*.5 1 /(DZZF(I,J,K-1)*DZF(I,J,K)) C(I,J,K)=-DT2*(KQ(I,J,K-1)+KQ(I,J,K)+2.*UMOL)*.5 1 /(DZZF(I,J,K-1)*DZF(I,J,K-1)) 100 CONTINUE C*********************************************************************** C * C THE FOLLOWING SECTION SOLVES THE EQUATION * C DT2*(KQ*Q2')' - Q2*(2.*DT2*DTEF+1.) = -Q2B * C * C*********************************************************************** C------ SURFACE AND BOTTOM B.C.S ------------ CONST1=16.6**.6666667*SEF DO 90 J=1,JMM1 DO 90 I=1,IMM1 KBB=KB(I,J) EE(I,J,1)=0. GG(I,J,1)=SQRT( (.5*(WUSURF(I,J)+WUSURF(I+1,J)))**2 1 +(.5*(WVSURF(I,J)+WVSURF(I,J+1)))**2 )*CONST1 UF(I,J,KBB)=SQRT( (.5*(WUBOT(I,J)+WUBOT(I+1,J)))**2 1 +(.5*(WVBOT(I,J)+WVBOT(I,J+1)))**2 )*CONST1 90 CONTINUE C--------------------------------------------------------------------- C RESTORE 10 DEG, 35 PPT AND 1.025 BIAS IN T, S AND RHO C--------------------------------------------------------------------- DO 101 K=1,KM-1 DO 101 J=1,JM DO 101 I=1,IM TP=T(I,J,K)+TBIAS SP=S(I,J,K)+SBIAS C Calculate pressure in units of decibars P=-GEE*1.025*ZZF(I,J,K)*.1 CC(I,J,K)=1449.1+.00821*P+4.55*TP -.045*TP**2 $ +1.34*(SP-35.) CC(I,J,K)=CC(I,J,K)/SQRT((1.-.01642*P/CC(I,J,K)) $ *(1.-0.40*P/CC(I,J,K)**2)) 101 CONTINUE DO 102 K=2,KM-1 DO 102 J=1,JM DO 102 I=1,IM Q2B(I,J,K)=ABS(Q2B(I,J,K)) Q2LB(I,J,K)=ABS(Q2LB(I,J,K)) BOYGR(I,J,K)=GEE*((RHO(I,J,K-1)-RHO(I,J,K))/DZZF(I,J,K-1)) 1 +GEE**2*2.*1.025/(CC(I,J,K-1)**2+CC(I,J,K)**2) 102 CONTINUE C------ CALC. T.K.E. PRODUCTION ----------------------------------- DO 120 K=2,KM-1 DO 120 J=1,JMM1 DO 120 I=1,IMM1 PROD(I,J,K)=KMT(I,J,K)*.25*SEF 1 *( (U(I,J,K)-U(I,J,K-1)+U(I+1,J,K)-U(I+1,J,K-1))**2 2 +(V(I,J,K)-V(I,J,K-1)+V(I,J+1,K)-V(I,J+1,K-1))**2 ) 3 /DZZF(I,J,K-1)**2 PROD(I,J,K)=PROD(I,J,K)+KH(I,J,K)*BOYGR(I,J,K) 120 CONTINUE DO 110 K=2,KM-1 DO 110 J=1,JM DO 110 I=1,IM STF=1. GHC=-2.5 IF(GH(I,J,K).LT.0.) STF=1.0-0.9*(GH(I,J,K)/GHC)**1.5 IF(GH(I,J,K).LT.GHC) STF=0.1 110 DTEF(I,J,K)=Q2B(I,J,K)*SQRT(Q2B(I,J,K))/(B1*Q2LB(I,J,K)+SMALL)*STF DO 140 K=2,KM-1 DO 140 J=1,JM DO 140 I=1,IM GG(I,J,K)=1./(A(I,J,K)+C(I,J,K)*(1.-EE(I,J,K-1)) 1 -(2.*DT2*DTEF(I,J,K)+1.) ) EE(I,J,K)=A(I,J,K)*GG(I,J,K)*FSM(I,J,K-1) GG(I,J,K)=(-2.*DT2*PROD(I,J,K) 1 +C(I,J,K)*GG(I,J,K-1)-UF(I,J,K))*GG(I,J,K) GG(I,J,K)=GG(I,J,K)*FSM(I,J,K-1) 140 CONTINUE DO 150 K=1,KM-1 KI=KM-K DO 150 J=1,JM DO 150 I=1,IM if(ki.le.kb(i,j)-1) +UF(I,J,KI)=EE(I,J,KI)*UF(I,J,KI+1)+GG(I,J,KI) 150 CONTINUE C*********************************************************************** C * C THE FOLLOWING SECTION SOLVES THE EQUATION * C DT2(KQ*Q2L')' - Q2L*(DT2*DTEF+1.) = -Q2LB * C * C*********************************************************************** DO 155 J=1,JM DO 155 I=1,IM EE(I,J,1)=0.0 GG(I,J,1)=0.0 KBB=KB(I,J) VF(I,J,KBB)=0. 155 CONTINUE DO 160 K=2,KM-1 DO 160 J=1,JM DO 160 I=1,IM kbb=kb(i,j) if(k.le.(kbb-1) ) +DTEF(I,J,K) =DTEF(I,J,K)*(1.+E2*((1./ABS(ZF(I,J,K)-ZF(I,J,1))+ 1 1./ABS(ZF(I,J,K)-ZF(I,J,KBB))) *L(I,J,K)/KAPPA)**2) GG(I,J,K)=1./(A(I,J,K)+C(I,J,K)*(1.-EE(I,J,K-1)) 1 -(DT2*DTEF(I,J,K)+1.)) EE(I,J,K)=A(I,J,K)*GG(I,J,K)*FSM(I,J,K-1) GG(I,J,K)=(DT2*(-PROD(I,J,K) 1 *L(I,J,K)*E1)+C(I,J,K)*GG(I,J,K-1)-VF(I,J,K))*GG(I,J,K) GG(I,J,K)=GG(I,J,K)*FSM(I,J,K-1) 160 CONTINUE DO 170 K=1,KM-1 KI=KM-K DO 170 J=1,JM DO 170 I=1,IM if(ki.le.kb(i,j)-1) +VF(I,J,KI)=EE(I,J,KI)*VF(I,J,KI+1)+GG(I,J,KI) 170 CONTINUE DO 180 K=2,KM-1 DO 180 J=1,JM DO 180 I=1,IM IF(UF(I,J,K).LT.SMALL.OR.VF(I,J,K).LT.SMALL) THEN UF(I,J,K)=SMALL VF(I,J,K)=SMALL ENDIF 180 CONTINUE C*********************************************************************** C * C THE FOLLOWING SECTION SOLVES FOR KM AND KH * C * C*********************************************************************** COEF1=A2*(1.-6.*A1/B1) COEF2=3.*A2*B2+18.*A1*A2 COEF3=A1*(1.-3.*C1-6.*A1/B1) COEF4=18.*A1*A1+9.*A1*A2 COEF5=9.*A1*A2 C NOTE THAT SM,SH LIMIT TO INFINITY WHEN GH APPROACHES 0.0288 DO 220 K=2,KM-1 DO 220 J=1,JM DO 220 I=1,IM L(I,J,K)=VF(I,J,K)/UF(I,J,K)*FSM(I,J,K-1) 220 CONTINUE DO 225 J=1,JM DO 225 I=1,IM L(I,J,1)=0. L(I,J,KM)=0. GH(I,J,1)=0. GH(I,J,KM)=0. 225 CONTINUE DO K=2,KM-1 DO J=1,JM DO I=1,IM GH(I,J,K)=L(I,J,K)**2/UF(I,J,K)*BOYGR(I,J,K) GH(I,J,K)=MIN(GH(I,J,K),.028) ENDDO ENDDO ENDDO DO 230 K=1,KM DO 230 J=1,JM DO 230 I=1,IM SH(I,J,K)=COEF1/(1.-COEF2*GH(I,J,K)) SM(I,J,K)=COEF3+SH(I,J,K)*COEF4*GH(I,J,K) SM(I,J,K)=SM(I,J,K)/(1.-COEF5*GH(I,J,K)) 230 CONTINUE C DO 280 K=1,KM DO 280 J=1,JM DO 280 I=1,IM KN(I,J,K)=L(I,J,K)*SQRT(ABS(Q2(I,J,K))) KQ(I,J,K)=(KN(I,J,K)*.41*SH(I,J,K)+KQ(I,J,K))*.5 KMT(I,J,K)=(KN(I,J,K)*SM(I,J,K)+KMT(I,J,K))*.5 KH(I,J,K)=(KN(I,J,K)*SH(I,J,K)+KH(I,J,K))*.5 280 CONTINUE RETURN END SUBROUTINE PROFT(F,WFSURF,SWRAD,FSURF,NBC,DT2) C INCLUDE 'comblkg.h' DIMENSION F(IM,JM,KM),WFSURF(IM,JM),FSURF(IM,JM) DIMENSION SWRAD(IM,JM),RAD(IM,JM,KM),TR(5),EXTC(5) C NTP=2 C NTP = 1 2 3 4 5 C JERLOV TYPE = I IA IB II III DATA TR/ .32 , .31 , .29 , .26 , .24 / DATA EXTC/ .037 , .042 , .056 , .073 , .127 / C UMOLPR=UMOL C*********************************************************************** C * C THE FOLLOWING SECTION SOLVES THE EQUATION * C DTI2*(KH*F')'-F=-FB * C * C*********************************************************************** DO 20 K=2,KM-1 DO 20 J=1,JM DO 20 I=1,IM A(I,J,K-1)=-DT2*(KH(I,J,K)+UMOLPR)/(DZF(I,J,K-1)*DZZF(I,J,K-1)) C(I,J,K)=-DT2*(KH(I,J,K)+UMOLPR)/(DZF(I,J,K)*DZZF(I,J,K-1)) 20 CONTINUE C----------------------------------------------------------------------- C NBC=1: SURF. B.C. IS WFSURF. NO SW RADIATIVE PENETRATION. C NBC=2; SURF. B.C. IS WFSURF+SWRAD*(1.-TR). TR*SWRAD PENETRATES WATER COLUMN C NBC=3; SURF. B.C. IS TSURF. NO SW RADIATIVE PENETRATION. C NBC=4; SURF. B.C. IS TSURF. TR*SWRAD PENETRATES WATER COLUMN C C NOTE THAT WTSURF (=WFSURF) AND SWRAD ARE NEG. VALUES WHEN WATER COLUMN IS C WARMING. C----------------------------------------------------------------------- C------------------------------------------------------------------ C Penetrative Radiation Calculation. At the bottom any C unattenuated radiation is deposited in the bottom layer. C------------------------------------------------------------------ DO 512 K=1,KM DO 512 J=1,JM DO 512 I=1,IM 512 RAD(I,J,K)=0. IF(NBC.EQ.2.OR.NBC.EQ.4.) THEN DO 511 K=1,KM-1 DO 511 J=1,JM DO 511 I=1,IM 511 RAD(I,J,K)=SWRAD(I,J)*TR(NTP)*EXP(EXTC(NTP)*ZF(I,J,K)) ENDIF GO TO (50,51,52,52), NBC 50 CONTINUE DO 500 J=1,JM DO 500 I=1,IM IF(FSM(I,J,1).gt.0.) then EE(I,J,1)=A(I,J,1)/(A(I,J,1)-1.0) GG(I,J,1)=-DT2*WFSURF(I,J)/(-DZF(I,J,1))-F(I,J,1) GG(I,J,1)=GG(I,J,1)/(A(I,J,1)-1.0) endif 500 CONTINUE GO TO 53 C 51 CONTINUE DO 510 J=1,JM DO 510 I=1,IM IF(FSM(I,J,1).gt.0.) then EE(I,J,1)=A(I,J,1)/(A(I,J,1)-1.0) GG(I,J,1)=DT2*(WFSURF(I,J)+SWRAD(I,J)*(1.0-TR(NTP)) 2 +RAD(I,J,1)-RAD(I,J,2)) 1 /(DZF(I,J,1))-F(I,J,1) GG(I,J,1)=GG(I,J,1)/(A(I,J,1)-1.0) endif 510 CONTINUE GO TO 53 C 52 CONTINUE DO 520 J=1,JM DO 520 I=1,IM EE(I,J,1)=0. 520 GG(I,J,1)=FSURF(I,J) C---------------------------------------------------------------------- 53 CONTINUE C C---------------------------------------------------------------------- DO 101 K=2,KM-2 DO 101 J=1,JM DO 101 I=1,IM GG(I,J,K)=1./(A(I,J,K)+C(I,J,K)*(1.-EE(I,J,K-1))-1.) EE(I,J,K)=A(I,J,K)*GG(I,J,K)*FSM(I,J,K) GG(I,J,K)=(C(I,J,K)*GG(I,J,K-1)-F(I,J,K) 1 +DT2*(RAD(I,J,K)-RAD(I,J,K+1))/DZF(I,J,K))*GG(I,J,K) GG(I,J,K)=GG(I,J,K)*FSM(I,J,K) 101 CONTINUE C----- BOTTOM ADIABATIC B.C. ------------------------------------------ DO 102 J=1,JM DO 102 I=1,IM KBB1=KB(I,J)-1 KBB2=KB(I,J)-2 KBB=KB(I,J) F(I,J,KBB1)=((C(I,J,KBB1)*GG(I,J,KBB2)-F(I,J,KBB1) 1 +DT2*(RAD(I,J,KBB1)-RAD(I,J,KBB))/DZF(I,J,KBB1)) 2 /(C(I,J,KBB1)*(1.-EE(I,J,KBB2))-1.)) 102 CONTINUE C---------------------------------------------------------------------- DO 105 K=2,KM-1 KI=KM-K DO 105 J=1,JM DO 105 I=1,IM IF(ki.lt.kb(i,j)-1) +F(I,J,KI)=(EE(I,J,KI)*F(I,J,KI+1)+GG(I,J,KI)) 105 CONTINUE RETURN END SUBROUTINE PROFU(DT2) INCLUDE 'comblkg.h' DATA UMOLB/1.36E-6 / C*********************************************************************** C * C THE FOLLOWING SECTION SOLVES THE EQUATION * C DT2*(KM*U')'-U=-UB * C * C*********************************************************************** C DO 90 K=1,KM DO 90 J=2,JM DO 90 I=2,IM UF(I,J,K)=UF(I,J,K)*DUM(I,J,K) 90 C(I,J,K)=(KMT(I,J,K)+KMT(I-1,J,K))*.5 DO 100 K=2,KM-1 DO 100 J=1,JM DO 100 I=2,IM A(I,J,K-1)=-DT2*(C(I,J,K)+UMOL ) 1 /(0.25*(DZF(I,J,K-1)+DZF(I-1,J,K-1)) 2 *(DZZF(I,J,K-1)+DZZF(I-1,J,K-1))) C(I,J,K)=-DT2*(C(I,J,K)+UMOL ) 1 /(0.25*(DZF(I,J,K)+DZF(I-1,J,K)) 2 *(DZZF(I,J,K-1)+DZZF(I-1,J,K-1))) 100 CONTINUE C DO 1011 J=1,JM DO 1011 I=2,IM EE(I,J,1)=A(I,J,1)/(A(I,J,1)-1.) GG(I,J,1)=(-DT2*WUSURF(I,J)/(-0.5*(DZF(I,J,1)+ 1 DZF(I-1,J,1)))-UF(I,J,1)) 2 /(A(I,J,1)-1.) 1011 CONTINUE DO 101 K=2,KM-2 DO 101 J=1,JM DO 101 I=2,IM GG(I,J,K)=1./(A(I,J,K)+C(I,J,K)*(1.-EE(I,J,K-1))-1.) EE(I,J,K)=A(I,J,K)*GG(I,J,K)*DUM(I,J,K) GG(I,J,K)=(C(I,J,K)*GG(I,J,K-1)-UF(I,J,K))*GG(I,J,K) GG(I,J,K)=GG(I,J,K)*DUM(I,J,K) 101 CONTINUE DO 102 J=2,JMM1 DO 102 I=2,IMM1 KBB1=KB(I,J)-1 KBB2=KB(I,J)-2 IF(FSM(I,J,KBB1).gt.0.) then TPS(I,J)=1.E-7+0.5*(CBC(I,J)+CBC(I-1,J)) 1 *SQRT(UB(I,J,KBB1)**2+(.25*(VB(I,J,KBB1) 2 +VB(I,J+1,KBB1)+VB(I-1,J,KBB1)+VB(I-1,J+1,KBB1)))**2) UF(I,J,KBB1)=(C(I,J,KBB1)*GG(I,J,KBB2)-UF(I,J,KBB1))/(TPS(I,J) 1 *DT2/(.5*(-DZF(I,J,KBB1)-DZF(I-1,J,KBB1))) 2 -1.-(EE(I,J,KBB2)-1.)*C(I,J,KBB1)) endif 102 CONTINUE DO 103 K=2,KM-1 KI=KM-K DO 103 J=2,JMM1 DO 103 I=2,IMM1 IF(ki.lt.kb(i,j)-1) +UF(I,J,KI)=(EE(I,J,KI)*UF(I,J,KI+1)+GG(I,J,KI))*DUM(I,J,KI) 103 CONTINUE C DO 104 J=2,JMM1 DO 104 I=2,IMM1 KBB1=KB(I,J)-1 104 WUBOT(I,J)=-TPS(I,J)*UF(I,J,KBB1) RETURN END SUBROUTINE PROFV(DT2) INCLUDE 'comblkg.h' DATA UMOLB/1.36E-6/ C*********************************************************************** C * C THE FOLLOWING SECTION SOLVES THE EQUATION * C DT2*(KM*U')'-U=-UB * C * C*********************************************************************** C DO 90 K=1,KM DO 90 J=2,JM DO 90 I=2,IM VF(I,J,K)=VF(I,J,K)*DVM(I,J,K) 90 C(I,J,K)=(KMT(I,J,K)+KMT(I,J-1,K))*.5 DO 100 K=2,KM-1 DO 100 J=2,JM DO 100 I=1,IM A(I,J,K-1)=-DT2*(C(I,J,K)+UMOL ) 1 /(.25*(DZF(I,J,K-1)+DZF(I,J-1,K-1)) 2 *(DZZF(I,J,K-1)+DZZF(I,J-1,K-1))) C(I,J,K)=-DT2*(C(I,J,K)+UMOL ) 1 /(.25*(DZF(I,J,K)+DZF(I,J-1,K)) 2 *(DZZF(I,J,K-1)+DZZF(I,J-1,K-1))) 100 CONTINUE C DO 1001 J=2,JM DO 1001 I=1,IM EE(I,J,1)=A(I,J,1)/(A(I,J,1)-1.) GG(I,J,1)=(-DT2*WVSURF(I,J)/(-0.5*(DZF(I,J,1)+ 1 DZF(I,J-1,1)))-VF(I,J,1)) 2 /(A(I,J,1)-1.) 1001 CONTINUE DO 101 K=2,KM-2 DO 101 J=2,JM DO 101 I=1,IM GG(I,J,K)=1./(A(I,J,K)+C(I,J,K)*(1.-EE(I,J,K-1))-1.) EE(I,J,K)=A(I,J,K)*GG(I,J,K)*DVM(I,J,K) GG(I,J,K)=(C(I,J,K)*GG(I,J,K-1)-VF(I,J,K))*GG(I,J,K) GG(I,J,K)=GG(I,J,K)*DVM(I,J,K) 101 CONTINUE DO 102 J=2,JMM1 DO 102 I=2,IMM1 KBB1=KB(I,J)-1 KBB2=KB(I,J)-2 if(FSM(I,J,KBB1).GT.0.) THEN TPS(I,J)=0.5*(CBC(I,J)+CBC(I,J-1)) 1 *SQRT((.25*(UB(I,J,KBB1)+UB(I+1,J,KBB1) 2 +UB(I,J-1,KBB1)+UB(I+1,J-1,KBB1)))**2+VB(I,J,KBB1)**2) VF(I,J,KBB1)=(C(I,J,KBB1)*GG(I,J,KBB2)-VF(I,J,KBB1))/(TPS(I,J) 1 *DT2/(.5*(-DZF(I,J,KBB1)-DZF(I,J-1,KBB1))) 2 -1.-(EE(I,J,KBB2)-1.)*C(I,J,KBB1)) endif 102 CONTINUE DO 103 K=2,KM-1 KI=KM-K DO 103 J=2,JMM1 DO 103 I=2,IMM1 IF(ki.lt.kb(i,j)-1) +VF(I,J,KI)=(EE(I,J,KI)*VF(I,J,KI+1)+GG(I,J,KI))*DVM(I,J,KI) 103 CONTINUE C DO 104 J=2,JMM1 DO 104 I=2,IMM1 KBB1=KB(I,J)-1 104 WVBOT(I,J)=-TPS(I,J)*VF(I,J,KBB1) RETURN END SUBROUTINE PRXY (LABEL,TIME,A,IM,ISKP,JM,JSKP,SCALA) C >>> C THIS WRITES A 2-D FIELD C TIME=TIME IN DAYS C A = ARRAY(IM,JM) TO BE PRINTED C ISKP=PRINT SKIP FOR I C JSKP=PRINT SKIP FOR J C SCALE=DIVISOR FOR VALUES OF A C C IMPLICIT HALF PRECISION (A-H,O-Z) DIMENSION A(IM,JM),NUM(150),LINE(150) CHARACTER LABEL*(*) DATA ZERO /1.E-12/ C SCALE=SCALA IF (SCALE.GT.ZERO) GO TO 160 AMX=ZERO DO 150 J=1,JM,JSKP DO 150 I=1,IM,ISKP AMX=MAX(ABS(A(I,J)),AMX) 150 CONTINUE IF(AMX.EQ.0.) THEN SCALEI=0. GOTO 165 ENDIF SCALE=10.0**(INT(LOG10(AMX)+1.E2)-103) 160 CONTINUE SCALEI=1.0/SCALE 165 CONTINUE WRITE(6,170) LABEL 170 FORMAT(1X,A40) WRITE(6,180) TIME,SCALE 180 FORMAT(' TIME =',F9.2,' DAYS MULTIPLY ALL VALUES BY',1PE10.3) DO 190 I=1,IM 190 NUM(I)=I IB=1 C 200 CONTINUE IE=IB+23*ISKP IF(IE.GT.IM) IE=IM WRITE(6,210) (NUM(I),I=IB,IE,ISKP) 210 FORMAT(/,2X,24I5,/) DO 260 J=1,JM,JSKP JWR=JM+1-J DO 220 I=IB,IE,ISKP 220 LINE(I)=INT(SCALEI*A(I,JWR)) WRITE(6,240) JWR,(LINE(I),I=IB,IE,ISKP) 240 FORMAT(1X,I3,24I5) 260 CONTINUE WRITE(6,280) 280 FORMAT(//) IF(IE.GE.IM) RETURN IB=IB+24*ISKP GO TO 200 END SUBROUTINE PRXYZ (LABEL,TIME,A,IM,ISKP,JM,JSKP,KM,SCALA) C >>> C C THIS WRITES HORIZONTAL LAYERS OF A 3-D FIELD C TIME=TIME IN DAYS C A = ARRAY(IM,JM,KM) TO BE PRINTED C ISKP=PRINT SKIP FOR I C JSKP=PRINT SKIP FOR J C SCALE=DIVISOR FOR VALUES OF A C DIMENSION A(IM,JM,KM),NUM(150),LINE(150),KP(4) CHARACTER LABEL*(*) DATA KE,KP /4, 1,2,3,6 / DATA ZERO /1.E-10/ C SCALE=SCALA IF (SCALE.GT.ZERO) GO TO 150 AMX=ZERO DO 140 KK=1,KE K=KP(KK) DO 140 J=1,JM,JSKP DO 140 I=1,IM,ISKP AMX=MAX(ABS(A(I,J,K)),AMX) 140 CONTINUE IF(AMX.EQ.0.) THEN SCALEI=0. GOTO 165 ENDIF SCALE=10.0**(INT(LOG10(AMX)+1.E2)-103) 150 CONTINUE SCALEI=1.0/SCALE 165 CONTINUE WRITE(6,160) LABEL 160 FORMAT(1X,A40) WRITE(6,170) TIME,SCALE 170 FORMAT(' TIME = ',F9.2,' DAYS MULTIPLY ALL VALUES BY',1PE10.3) DO 180 I=1,IM 180 NUM(I)=I C DO 500 KK=1,KE K=KP(KK) WRITE(6,190) K 190 FORMAT(3X,/7H LAYER ,I2) IB=1 C 200 CONTINUE IE=IB+23*ISKP IF(IE.GT.IM) IE=IM WRITE(6,220) (NUM(I),I=IB,IE,ISKP) 220 FORMAT(/,2X,24I5,/) DO 260 J=1,JM,JSKP JWR=JM+1-J DO 230 I=IB,IE,ISKP 230 LINE(I)=INT(SCALEI*A(I,JWR,K)) WRITE(6,240) JWR,(LINE(I),I=IB,IE,ISKP) 240 FORMAT(1X,I3,24I5) 260 CONTINUE WRITE(6,280) 280 FORMAT(//) IF(IE.GE.IM) GO TO 500 IB=IB+24*ISKP GO TO 200 500 CONTINUE RETURN END SUBROUTINE PRXZ(LABEL,TIME,A,IM,ISKP,JM,J1,J2,KM,SCALA,DT) DIMENSION A(IM,JM,KM),NUM(100),LINE(100),IDT(100), 1 ZZ(KM),DT(IM,JM) CHARACTER LABEL*(*) DATA ZERO/1.E-10/ C C THIS SUBROUTINE WRITES VERTICAL SECTIONS OF A 3-D FIELD C C TIME=TIME IN DAYS C A= ARRAY(IM,JM,KM) TO BE PRINTED C ISPL=PRINT SKIP FOR I C JSPL=PRINT SKIP FOR J C SCALE=DIVISOR FOR VALUES OF A C SCALE=SCALA IF (SCALE.GT.ZERO) GO TO 150 AMX=ZERO DO 140 K=1,KM DO 140 J=1,JM DO 140 I=1,IM,ISKP AMX=MAX(ABS(A(I,J,K)),AMX) 140 CONTINUE IF(AMX.EQ.0.) THEN SCALEI=0. GOTO 165 ENDIF SCALE=10.0**(INT(LOG10(AMX)+1.E2)-103) 150 CONTINUE SCALEI=1./SCALE 165 CONTINUE WRITE(6,9980) LABEL 9980 FORMAT(' ',A30) WRITE(6,9981) TIME,SCALE 9981 FORMAT(' TIME =',F9.3,' MULTIPLY ALL VALUES BY ',E8.1) DO 10 JPP=1,2 J=J1 IF(JPP.EQ.2) J=J2 IF(J.EQ.0) RETURN WRITE(6,30) J 30 FORMAT(3X,/' SECTION J =',I3) IB=1 IE=IB+23*ISKP 50 CONTINUE IF(IE.GT.IM) IE=IM DO 100 I=IB,IE,ISKP IDT(I)=DT(I,J) 100 NUM(I)=I WRITE(6,9990) (NUM(I),I=IB,IE,ISKP) 9990 FORMAT(/,' I = ',24I5,/) WRITE(6,9991) (IDT(I),I=IB,IE,ISKP) 9991 FORMAT(8X,'D =',24I5.0,/,' ') DO 200 K=1,KM DO 190 I=IB,IE,ISKP 190 LINE(I)=SCALEI*A(I,J,K) WRITE(6,9966) K,(LINE(I),I=IB,IE,ISKP) 9966 FORMAT(1X,I2,2X,6X,24I5) 200 CONTINUE WRITE(6,9984) IF(IE.EQ.IM) GO TO 10 IB=IB+24*ISKP IE=IB+23*ISKP GO TO 50 9984 FORMAT(//) 10 CONTINUE RETURN END SUBROUTINE PRYZ(LABEL,TIME,A,IM,I1,I2,JM,JSKP,KM,SCALA,DT) DIMENSION A(IM,JM,KM),NUM(200),LINE(200),DT(IM,JM),ZZ(KM) DIMENSION IP(2) CHARACTER LABEL*(*) DATA ZERO/1.E-10/ C C THIS SUBROUTINE WRITES VERTICAL SECTIONS OF A 3-D FIELD C C TIME=TIME IN DAYS C A= ARRAY(IM,JM,KM) TO BE PRINTED C ISKP=PRINT SKIP FOR I C JSKP=PRINT SKIP FOR J C SCALE=DIVISOR FOR VALUES OF A C SCALE=SCALA IF (SCALE.GT.ZERO) GO TO 150 AMX=ZERO DO 140 K=1,KM DO 140 J=1,JM,JSKP DO 140 IPP=1,2 I=I1 IF(IPP.EQ.2) I=I2 AMX=MAX(ABS(A(I,J,K)),AMX) 140 CONTINUE IF(AMX.EQ.0.) THEN SCALEI=0. GOTO 165 ENDIF SCALE=10.0**(INT( LOG10(AMX)+1.E2)-103) 150 CONTINUE SCALEI=1./SCALE 165 CONTINUE WRITE(6,9980) LABEL 9980 FORMAT(' ',A30) WRITE(6,9981) TIME,SCALE 9981 FORMAT(' TIME =',F9.2,' MULTIPLY ALL VALUES BY ',E8.1) DO 10 IPP=1,2 I=I1 IF(IPP.EQ.2) I=I2 WRITE(6,30) I 30 FORMAT(3X,/' SECTION I =',I3) JB=1 JE=JB+23*JSKP 50 CONTINUE IF(JE.GT.JM) JE=JM DO 100 J=JB,JE,JSKP 100 NUM(J)=J WRITE(6,9990) (NUM(J),J=JB,JE,JSKP) 9990 FORMAT(/,' J = ',24I5,/) WRITE(6,9991) (DT(I,J),J=JB,JE,JSKP) 9991 FORMAT(7X,'D =',24F5.0,/,' ') DO 200 K=1,KM DO 190 J=JB,JE,JSKP 190 LINE(J)=SCALEI*A(I,J,K) WRITE(6,9966) K,(LINE(J),J=JB,JE,JSKP) 9966 FORMAT(1X,I2,7X,24I5) 200 CONTINUE WRITE(6,9984) IF(JE.EQ.JM) GO TO 10 JB=JB+24*JSKP JE=JB+23*JSKP GO TO 50 9984 FORMAT(//) 10 CONTINUE RETURN END SUBROUTINE SLPMIN(H,IM,JM,FSM,SL) DIMENSION H(IM,JM),FSM(IM,JM) DIMENSION SL(IM,JM) C This subroutine limits the maximum difference in H divided C by twice the average of H of adjacent cells. C The maximum possible value is unity. C SLMIN=0.2 C DO 100 LOOP=1,10 C sweep right DO 3 J=2,JM-1 DO 1 I=2,IM-1 IF(FSM(I,J).EQ.0..OR.FSM(I+1,J).EQ.0.) GOTO 1 SL(I,J)=ABS(H(I+1,J)-H(I,J))/(H(I,J)+H(I+1,J)) IF(SL(I,J).LT.SLMIN) GOTO 1 DH=0.5*(SL(I,J)-SLMIN)*(H(I,J)+H(I+1,J)) SN=-1. IF(H(I+1,J).GT.H(I,J)) SN=1. H(I+1,J)=H(I+1,J)-SN*DH H(I,J)=H(I,J)+SN*DH 1 CONTINUE C sweep left DO 2 I=IM-1,2,-1 IF(FSM(I,J).EQ.0..OR.FSM(I+1,J).EQ.0.) GOTO 2 SL(I,J)=ABS(H(I+1,J)-H(I,J))/(H(I,J)+H(I+1,J)) IF(SL(I,J).LT.SLMIN) GOTO 2 DH=0.5*(SL(I,J)-SLMIN)*(H(I,J)+H(I+1,J)) SN=-1. IF(H(I+1,J).GT.H(I,J)) SN=1. H(I+1,J)=H(I+1,J)-SN*DH H(I,J)=H(I,J)+SN*DH 2 CONTINUE 3 CONTINUE C sweep up DO 13 I=2,IM-1 DO 11 J=2,JM-1 IF(FSM(I,J).EQ.0..OR.FSM(I,J+1).EQ.0.) GOTO 11 SL(I,J)=ABS(H(I,J+1)-H(I,J))/(H(I,J)+H(I,J+1)) IF(SL(I,J).LT.SLMIN) GOTO 11 DH=0.5*(SL(I,J)-SLMIN)*(H(I,J)+H(I,J+1)) SN=-1. IF(H(I,J+1).GT.H(I,J)) SN=1. H(I,J+1)=H(I,J+1)-SN*DH H(I,J)=H(I,J)+SN*DH 11 CONTINUE C sweep down DO 12 J=JM-1,2,-1 IF(FSM(I,J).EQ.0..OR.FSM(I,J+1).EQ.0.) GOTO 12 SL(I,J)=ABS(H(I,J+1)-H(I,J))/(H(I,J)+H(I,J+1)) IF(SL(I,J).LT.SLMIN) GOTO 12 DH=0.5*(SL(I,J)-SLMIN)*(H(I,J)+H(I,J+1)) SN=-1. IF(H(I,J+1).GT.H(I,J)) SN=1. H(I,J+1)=H(I,J+1)-SN*DH H(I,J)=H(I,J)+SN*DH 12 CONTINUE 13 CONTINUE C 100 CONTINUE RETURN END SUBROUTINE VERTVL(DTI2) INCLUDE 'comblkg.h' DIMENSION XFLUX(IM,JM,KM),YFLUX(IM,JM,KM) EQUIVALENCE (XFLUX,A),(YFLUX,C) C C CALCULATE NEW VERTICAL VELOCITY C C REESTABLISH BOUNDARY CONDITIONS do k=1,km do j=1,jm do i=1,im xflux(i,j,k) = 0.0 yflux(i,j,k) = 0.0 enddo enddo enddo DO 100 K=1,KM-1 DO 100 J=2,JM DO 100 I=2,IM XFLUX(I,J,K) 1 =.25*(DY(I,J)+DY(I-1,J))*(DZ(I,J,K)+DZ(I-1,J,K))*U(I,J,K) 100 CONTINUE DO 120 K=1,KM-1 DO 120 J=2,JM DO 120 I=1,IM YFLUX(I,J,K) 1 =.25*(DX(I,J)+DX(I,J-1))*(DZ(I,J,K)+DZ(I,J-1,K))*V(I,J,K) 120 CONTINUE DO 125 J=1,JM DO 125 I=1,IM 125 W(I,J,1)=0.0 DO 710 K=1,KM-1 DO 710 J=2,JMM1 DO 710 I=2,IMM1 W(I,J,K+1)=W(I,J,K) 1 +(XFLUX(I+1,J,K)-XFLUX(I,J,K) 2 +YFLUX(I,J+1,K)-YFLUX(I,J,K))/(DX(I,J)*DY(I,J)) 3 +FSM(I,J,K)*(DZF(I,J,K)-DZB(I,J,K))/DTI2 710 CONTINUE RETURN END C SUBROUTINE PRXY1(LABEL,TIME,A,IM,ISKP,JM,JSKP,SCALA) C >>> C THIS WRITES A 2-D FIELD C TIME=TIME IN DAYS C A = ARRAY(IM,JM) TO BE PRINTED C ISKP=PRINT SKIP FOR I C JSKP=PRINT SKIP FOR J C SCALE=DIVISOR FOR VALUES OF A C C IMPLICIT HALF PRECISION (A-H,O-Z) DIMENSION A(IM,JM),NUM(150),LINE(150) CHARACTER LABEL*(*) DATA ZERO /1.E-12/ C SCALE=SCALA IF (SCALE.GT.ZERO) GO TO 160 AMX=ZERO DO 150 J=1,JM,JSKP DO 150 I=1,IM,ISKP AMX=MAX(ABS(A(I,J)),AMX) 150 CONTINUE IF(AMX.EQ.0.) THEN SCALEI=0. GOTO 165 ENDIF SCALE=10.0**(INT(LOG10(AMX)+1.E2)-103) 160 CONTINUE SCALEI=1.0/SCALE 165 CONTINUE WRITE(6,170) LABEL 170 FORMAT(1X,A40) WRITE(6,180) TIME,SCALE 180 FORMAT(' TIME =',F9.2,' DAYS MULTIPLY ALL VALUES BY',1PE10.3) DO 190 I=1,IM 190 NUM(I)=mod(I,10) IB=1 C 200 CONTINUE IE=IB+70*ISKP IF(IE.GT.IM) IE=IM WRITE(6,210) (NUM(I),I=IB,IE,ISKP) 210 FORMAT(/,4X,70I1,/) DO 260 J=1,JM,JSKP JWR=JM+1-J DO 220 I=IB,IE,ISKP 220 LINE(I)=INT(SCALEI*A(I,JWR)) WRITE(6,240) JWR,(LINE(I),I=IB,IE,ISKP) 240 FORMAT(I3,1x,70I1) 260 CONTINUE WRITE(6,280) 280 FORMAT(//) IF(IE.GE.IM) RETURN IB=IB+70*ISKP GO TO 200 END SUBROUTINE PRXY2(LABEL,TIME,A,IM,ISKP,JM,JSKP,SCALA) C >>> C THIS WRITES A 2-D FIELD C TIME=TIME IN DAYS C A = ARRAY(IM,JM) TO BE PRINTED C ISKP=PRINT SKIP FOR I C JSKP=PRINT SKIP FOR J C SCALE=DIVISOR FOR VALUES OF A C C IMPLICIT HALF PRECISION (A-H,O-Z) DIMENSION A(IM,JM),NUM(150),LINE(150) CHARACTER LABEL*(*) DATA ZERO /1.E-12/ C SCALE=SCALA IF (SCALE.GT.ZERO) GO TO 160 AMX=ZERO DO 150 J=1,JM,JSKP DO 150 I=1,IM,ISKP AMX=MAX(ABS(A(I,J)),AMX) 150 CONTINUE IF(AMX.EQ.0.) THEN SCALEI=0. GOTO 165 ENDIF SCALE=10.0**(INT(LOG10(AMX)+1.E2)-103) 160 CONTINUE SCALEI=1.0/SCALE 165 CONTINUE WRITE(6,170) LABEL 170 FORMAT(1X,A40) WRITE(6,180) TIME,SCALE 180 FORMAT(' TIME =',F9.2,' DAYS MULTIPLY ALL VALUES BY',1PE10.3) DO 190 I=1,IM 190 NUM(I)=mod(I,10) IB=1 C 200 CONTINUE IE=IB+35*ISKP IF(IE.GT.IM) IE=IM WRITE(6,210) (NUM(I),I=IB,IE,ISKP) 210 FORMAT(/,4X,35I2,/) DO 260 J=1,JM,JSKP JWR=JM+1-J DO 220 I=IB,IE,ISKP 220 LINE(I)=INT(SCALEI*A(I,JWR)) WRITE(6,240) JWR,(LINE(I),I=IB,IE,ISKP) 240 FORMAT(I3,1x,35I2) 260 CONTINUE WRITE(6,280) 280 FORMAT(//) IF(IE.GE.IM) RETURN IB=IB+35*ISKP GO TO 200 END SUBROUTINE GSIGTOZ(ZZ,T,KB,TLEV,ZLEV,FSM,IM,JM,KM) C------------------------------------------------------------------- C THIS ROUTINE LINERLY INTERPOLATES TLEV AT THE LEVEL, C ZLEV, FROM T LOCATED ON SIGMA LEVELS, ZZ. C NOTE THAT A NEW MASK ,FSM, IS CREATED. C------------------------------------------------------------------- INTEGER KB(IM,JM) REAL ZZ(IM,JM,KM),T(IM,JM,KM),TLEV(IM,JM),FSM(IM,JM) DO j=1,JM do i=1,im TLEV(I,J)=0.0 FSM(I,J)=0.0 enddo enddo K=1 DO 200 J=1,JM DO 200 I=1,IM 100 CONTINUE IF(K.LT.1) GO TO 200 IF(K.GT.KB(I,J)-1) THEN K=KB(I,J)-2 GO TO 200 ENDIF C C FIND K AND K+1 INTERVAL THAT BRACKETS ZLEV, THEN INTERPOLATE C IF(ZLEV.LT.ZZ(I,J,K)) THEN IF(ZLEV.GE.ZZ(I,J,K+1)) THEN TLEV(I,J)=T(I,J,K)+(ZLEV-ZZ(I,J,K)) 1 *(T(I,J,K+1)-T(I,J,K))/(ZZ(I,J,K+1)-ZZ(I,J,K)) FSM(I,J)=1.0 IF(K.EQ.KB(I,J)-1) TLEV(I,J)=0 GO TO 200 ELSE K=K+1 GO TO 100 ENDIF ELSE K=K-1 GO TO 100 ENDIF 200 CONTINUE RETURN END C SUBROUTINE SMOOTH(A,B,FSM,IM,JM,NITS) C------------------------------------------------------------------- C THIS ROUTINE SMOOTHS DATA WITH A FIVE POINT LAPLACIAN FILTER. C------------------------------------------------------------------- DIMENSION A(IM,JM),FSM(IM,JM),B(IM,JM) DO 100 N=1,NITS do j=1,jm do i=1,im b(i,j)=a(i,j) enddo enddo DO 200 J=2,JM-1 DO 200 I=2,IM-1 SMFAC=FSM(I+1,J)+FSM(I,J-1)+FSM(I-1,J)+FSM(I,J+1)+1.E-5 A(I,J)=B(I,J)+(.5/SMFAC) 1 *(B(I+1,J)*FSM(I+1,J)+B(I,J-1)*FSM(I,J-1) 2 +B(I-1,J)*FSM(I-1,J)+B(I,J+1)*FSM(I,J+1) 3 -SMFAC*B(I,J)) A(I,J)=A(I,J)*FSM(I,J) 200 CONTINUE 100 CONTINUE RETURN END C SUBROUTINE ARMEAN C GDEM VERTICAL LEVELS INCLUDE 'comblkg.h' PARAMETER(KS=35) DIMENSION ZL(KS),TLEV(IM,JM) DATA ZL /0.,10.,20.,30.,50.,75.,100.,125.,150.,200., 2 250.,300.,400.,500.,600.,700.,800.,900.,1000.,1100., 3 1200.,1300.,1400.,1500.,1750.,2000.,2500.,3000.,4000.,5000., 4 6000.,7000.,8000.,9000.,10000./ DO K=1,KS-1 ZSM=-.5*(ZL(K)+ZL(K+1)) CALL GSIGTOZ(ZZ,TB,KB ,TLEV,ZSM,TPS,IM,JM,KM) ENDDO RETURN END