C C C 2d-channel model, alongshore pressure C gradient or wind stress driven. C T = RHO so no eqn of state. C No diffusion for T in this model C Hua 6/23/95 C insert changes in DO 96 & 199 loops to conform to 3-D model C correct error in PSI calculation C glm 1/96 C C This was last used to do the research described in: C Mellor, G. L. and X. H. Wang, Pressure compensation and the bottom C boundary layer, J. Phys. Oceanogr., 26(10), 2214-2222, 1996. C However, at this date, I am no longer familiar with the details C of the problem setup. C C PROGRAM MAIN C REAL KM,KH,KQ,L,IRX,IRHO,MKE,MPE PARAMETER (IM=41,KB=30,KS=20,IMM1=IM-1,KBM1=KB-1) PARAMETER (IMM2=IM-2,KBM2=KB-2) PARAMETER (IM2=IM*2,IMK=IM*KB,IMKM2=IM*KBM2) PARAMETER (IMKM1=IM*KBM1) COMMON/BLKCON/ 1 INT,IPRINT,DTE,DTI,TPRNU,UMOL,PR, 2 GRAV,TIME,RAMP C---------------- 1-D ARRAYS -------------------------------------- COMMON/BLK1D/ 1 DZR(KB),Z(KB),ZZ(KB),DZ(KB),DZZ(KB) C---------------- 2-D ARRAYS -------------------------------------- COMMON/BLK2D/H(IM),DX(IM),D(IM),DT(IM), 1 ART(IM),ARU(IM),ARV(IM),CBC(IM), 2 ALON(IM),ALAT(IM),ANG(IM), 3 DUM(IM),DVM(IM),FSM(IM),COR4(IM),CURV42D(IM), 4 WUSURF(IM),WVSURF(IM),WUBOT(IM),WVBOT(IM), 5 WTSURF(IM),WSSURF(IM),TPS(IM),AAM2D(IM), 6 UAF(IM),UA(IM),UAB(IM),VAF(IM),VA(IM), 7 VAB(IM),ELF(IM),EL(IM),ELB(IM), 8 ETF(IM),ET(IM),ETB(IM),FLUXUA(IM),FLUXVA(IM), 9 EGF(IM),EGB(IM), ELBB(IM) C---------------- 3-D ARRAYS -------------------------------------- COMMON/BLK3D/ 1 A(IM,KB),C(IM,KB),VH(IM,KB),VHP(IM,KB), 1 UF(IM,KB),VF(IM,KB), 2 KM(IM,KB),KH(IM,KB),KQ(IM,KB),L(IM,KB), 3 Q2(IM,KB),Q2B(IM,KB),AAM(IM,KB), 4 Q2L(IM,KB),Q2LB(IM,KB), 5 U(IM,KB),UB(IM,KB),W(IM,KB), 6 V(IM,KB),VB(IM,KB), 7 T(IM,KB),TB(IM,KB), 8 S(IM,KB),SB(IM,KB), 9 RHO(IM,KB),DTEF(IM,KB),RMEAN(IM,KB) C----------- 1 AND 2-D BOUNDARY VALUE ARRAYS ------------------------ COMMON/BDRY/ 1 TBE(KB),TBN(IM,KB),TBS(IM,KB),SBN(IM,KB),SBE(KB), 2 SBS(IM,KB),VABN(IM),VABS(IM),UGS(KB),VGS(IM,KB), 3 ELN(IM),ELS(IM),COVRHN(IM) DIMENSION UTB(IM),VTB(IM),UTF(IM),VTF(IM), 1 ADVUA(IM),ADVVA(IM),TSURF(IM),SSURF(IM), 2 DRHOX(IM,KB),TRNU(IM),ADVX(IM,KB),ADVY(IM,KB), 3 TCLIM(IM,KB),SCLIM(IM,KB),ADVUU(IM),ADVVV(IM), 4 IRHO(IM,KS),IRX(KS),X(IM),RX(KB),RHOX(IM,KB) DIMENSION ZLEV(KS),RHOLEV(KS),RHOZ(KB), 1 ZD(KB),ZZD(KB) DIMENSION SAVRHO(IM,KB),PSI(IM,KB), 1 ADVTR(IM),CORTR(IM),STRTR(IM),PGRTR(IM),TENTR(IM) DIMENSION ADVECY(IM) REAL ISPI,ISP2I C MODE=3 TIME=0. DTI=1800. NREAD=0 DAYS=1 PRTD1=1 PRTD2=1 SWTCH=3600 ISPLIT=30 ISPADV=1 HORCON=0.2 TPRNU=1. SMOTH=.1 INT=0 AAA=0. INCLUDE 'params' C-------------------------------------------------------------------------- C READ NAMEIMST FOR PROBLEM PARAMETERS@D C DTI = INTERNAL TIME STEP C ISPLIT = DTI/DTE C DTI = INTERNAL TIME STEP 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 READ(5,PRMTR) C-------------------------------------------------------------------------- DTE=DTI/FLOAT(ISPLIT) DTE2=DTE*2 DTI2=DTI*2 IPRINT=PRTD1*24*3600/IFIX(DTI)+.5 ISWTCH=SWTCH*24*3600/IFIX(DTI) IEND=DAYS*24*3600/IFIX(DTI)+.5 WRITE(6,7030) MODE,DTE,DTI,IEND,ISPLIT,ISPADV,IPRINT,SMOTH,HORCON, 1 AAA,TPRNU 7030 FORMAT(//,' MODE =',I3, 1 ' DTE =',F7.2,' DTI =',F7.2,' IEND =',I6, 2 ' ISPLIT =',I6,' ISPADV =',I6,' IPRINT =',I6,/, 3 ' SMOTH =',F7.2,' HORCON =',F7.3,' AAA =',F8.2,' TPRNU = 4',F6.2,/) WRITE(6,'('' UMOL ='',E12.3,'' PR ='',F10.1)') UMOL,PR WRITE(6,'('' NREAD ='',I5,//)') NREAD DATA UAF/IM*0.E0/,UA/IM*0.E0/,UAB/IM*0.E0/ DATA VAF/IM*0.E0/,VA/IM*0.E0/,VAB/IM*0.E0/ DATA UTB/IM*0.E0/,VTB/IM*0.E0/UTF/IM*0.E0/,VTF/IM*0.E0/ DATA ELF/IM*0.E0/,EL/IM*0.E0/,ELB/IM*0.E0/,ETB/IM*0.E0/ DATA ET/IM*0.E0/,ETF/IM*0.E0/,ADVUA/IM*0.E0/,AAM2D/IM*50.E0/ DATA ADVVA/IM*0.E0/,TRNU/IM*0.E0/ DATA ADVUU/IM*0.E0/,ADVVV/IM*0.E0/,FLUXUA/IM*0./ DATA FLUXVA/IM*0.E0/,WUSURF/IM*0.E0/ c DATA WVSURF/IM*1.0E-4/,WTSURF/IM*0.E0/,WSSURF/IM*0.E0/ DATA WVSURF/IM*0.0E0/,WTSURF/IM*0.E0/,WSSURF/IM*0.E0/ DATA WUBOT/IM*0.E0/,WVBOT/IM*0.E0/ DATA A/IMK*0.E0/,C/IMK*0.E0/,VH/IMK*0.E0/,VHP/IMK*0.E0/ DATA UF/IMK*0.E0/,U/IMK*0.E0/,UB/IMK*0.E0/,W/IMK*0.E0/ DATA VF/IMK*0.E0/,V/IMK*0.E0/,VB/IMK*0.E0/ DATA ADVX/IMK*0.0/,ADVY/IMK*0.0/ DATA T/IMK*0.E0/,TB/IMK*0.E0/,S/IMK*0.E0/,SB/IMK*0.E0/ DATA RHO/IMK*0.E0/,Q2B/IMK*0.E0/,Q2LB/IMK*0.E0/ DATA TCLIM/IMK*0.E0/,SCLIM/IMK*0.E0/,RMEAN/IMK*0.E0/ DATA DRHOX/IMK*0.E0/ DATA PI/3.1416E0/,RAMP/1.E0/,SMALL/1.E-10/,TIME0/0.E0/ C C---------------------------------------------------------------------- C ESTABIMSH PROBLEM CHARACTERISTICS C ****** ALL UNITS IN M.K.S. SYSTEM ****** C F,BLANK AND B REFERS TO FORWARD,CENTRAL AND BACKWARD TIME LEVELS. C---------------------------------------------------------------------- DATA BETA/1.98E-11/,GRAV/9.806E0/,UMOL/1.E-4/,PR/100./ C DAYI=1.E0/86400.E0 ISPI=1.E0/FLOAT(ISPLIT) ISP2I=1.E0/(2.E0*FLOAT(ISPLIT)) C------------------------------------------------------------------------ C READ IN GRID DATA AND INITIAL CONDITIONS C------------------------------------------------------------------------ C C C REWIND(40) C REWIND(50) C REWIND(60) C READ(40) KBB,Z,ZZ,DZ,DZZ,IMM,MM,ALON,ALAT,DX,DY,H,COR4,ANG, C 1 ART,ARU,ARV,TCLIM,SCLIM,RMEAN,TBN,TBE,TBS,SBN,SBE,SBS, C 2 ELN,ELE,ELS,VABN,UABE,VABS C CALL VABFIX C***************************************************************************** C READ(60) WUSURF,WVSURF,WTSURF C READ(50) TIME0,UB,VB,UAB,VAB,TB,SB,ELB,Q2B,Q2LB,KM,KH C TSURF(*,*)=TB(*,*,1) C SSURF(*,*)=SB(*,*,1) C*SB* SETUP GRID AND TEST INITIAL CONDITIONS WO=100.E3 DO=4000. CALL DEPTH(Z,ZZ,DZ,DZZ,DZR,KB,KB-1) DO 100 I=1,IM DX(I)=20.E3 100 X(I)=DX(I)*FLOAT(I-(IM+1)/2) DO 150 I=1,IM 150 H(I)=DO*(1.-EXP(-((X(I)-4.E5)/WO)**2) * -EXP(-((X(I)+4.E5)/WO)**2)) H(2)=H(3)/2. H(40)=H(39)/2. H(1)=1. H(IM)=1. DO 30 I=1,IM DT(I)=H(I) FSM(I)=1.E0 DUM(I)=1.E0 DVM(I)=1.E0 IF(H(I).GT.1.E0) GO TO 30 FSM(I)=0.E0 DUM(I)=0.E0 DVM(I)=0.E0 30 CONTINUE DO 35 I=1,IM IF(FSM(I).EQ.0.E0.AND.FSM(I).NE.0.E0) DVM(I)=0.E0 35 CONTINUE DO 40 I=1,IMM1 IF(FSM(I).EQ.0.E0.AND.FSM(I+1).NE.0.E0) DUM(I+1)=0.E0 40 CONTINUE C CALL SLPMIN(H,IM,FSM) C WRITE(6,'('' I X H'')') WRITE(6,'(I5,F12.2,F10.3)') (I,X(I),H(I),I=1,IM) C23456 | | | | | | | xxxxxxx| LAB=1 CALL PRXZ(' FSM ',TIME,FSM,IM,1,1,0.,DT,0.,LAB) CALL PRXZ(' DUM ',TIME,DUM,IM,1,1,0.,DT,0.,LAB) CALL PRXZ(' DVM ',TIME,DVM,IM,1,1,0.,DT,0.,LAB) DO 11 K=1,KB DO 11 I=1,IM Q2B(I,K)=1.E-4 KH(I,K)=0.1E-2 KH(I,K)=0. KM(I,K)=KH(I,K) 11 AAM(I,K)=AAA DO 12 I=1,IM ART(I)=DX(I) ARU(I)=DX(I) ARV(I)=DX(I) 12 CONTINUE DO 14 I=1,IM COR4(I)=1.E-4/4. 14 CONTINUE DO 1 K=1,KS DO 1 I=1,IM ZLEV(K)=FLOAT(K-1)*DO/FLOAT(KS-1) IRHO(I,K)=.028-.004*EXP(-ZLEV(K)/2000.) c IRHO(I,K)=.026 C::North Atlantic Ocean Density Distribution C C IRHO(I,K)=.028-.003*(1.-ZLEV(K)/1000.) 1 CONTINUE DO I=1,IM VAB(I)=0.5 DO K=1,KB VB(I,K)=0.5 ENDDO ENDDO FLOW=0. ELBB(1)=1 DO I=2,IMM1 FLOW=FLOW+VAB(I)*DX(I)*H(I) ELBB(I+1)=ELBB(I)+4.*COR4(I)*VB(I,1)*DX(I)/GRAV ENDDO DO I=2,IMM1 ELBB(I)=ELBB(I)-ELBB(IM/2) ENDDO C TRY THIS FLOWB=FLOW DO I=1,IM ELB(I)=0. VAB(I)=0. DO K=1,KB VB(I,K)=0. DRHOX(I,K)=0. ENDDO ENDDO CALL ZTOSIG(ZLEV,IRHO,ZZ,H,RHO,IM,KS,KB) CALL PRXZ(' RHO ',0.,RHO,IM,1,KB,0.00001,DT,ZZ,LAB) DO 15 K=1,KBM1 DO 15 I=1,IM RHO(I,K)=RHO(I,K)*FSM(I) SAVRHO(I,K) =RHO(I,K) TB(I,K)=RHO(I,K) SB(I,K)=0.E0 15 CONTINUE DO 16 I=1,IM DT(I)=H(I) TSURF(I)=TB(I,1) 16 SSURF(I)=SB(I,1) DO 17 K=1,KBM1 DO 17 I=1,IM RMEAN(I,K)=RHO(I,K) TCLIM(I,K)=TB(I,K) 17 SCLIM(I,K)=SB(I,K) RHOBOT=RHO(21,KBM1) WRITE(6,'('' RHOBOT ='',E14.6 )') RHOBOT C INTRODUCE SIMPLE WIND STRESS, VALUE IS NEGATIVE FOR WESTERLY WIND DO 18 I=1,IM WUSURF(I)=0. 18 CONTINUE C END OF SETUP C DO 20 K=1,KBM1 20 DZR(K)=1./DZ(K) PERIOD=DAYI*(2.*PI)/(COR4((IM+1)/2)*4.) DO 21 I=1,IM 21 CURV42D(I)=COR4(I) C CALL PRXZ(' FSM ',TIME,FSM,IM,1,1,0.,DT,0.,LAB) CALL PRXZ(' DUM ',TIME,DUM,IM,1,1,0.,DT,0.,LAB) CALL PRXZ(' DVM ',TIME,DVM,IM,1,1,0.,DT,0.,LAB) DO 45 I=1,IM Z0B=.01E0 CBCMIN=.0025E0 45 CBC(I)=AMAX1(CBCMIN,.16E0/ALOG((ZZ(KBM1)-Z(KB))*H(I) 1 /Z0B)**2) C DO 46 I=1,IM 46 CBC(I)=CBC(I)*FSM(I) IP=IM ISK=1 DO 47 I=1,IM 47 TPS(I)=0.5E0/SQRT(1.E0/DX(I)**2) 1 /SQRT(GRAV*H(I))*FSM(I) CALL PRXZ(' CFL DEL TIME',TIME,TPS,IM,ISK,1 ,0. ,DT,ZZ,LAB) CALL PRXZ(' ELB ',TIME,ELB,IM,1,1,0.,DT,0.,LAB) CALL PRXZ(' TB ',TIME,TB,IM,ISK,KB,.0001,DT,ZZ,LAB) C CALL PRXZ(' SB ',TIME,SB,IM,ISK,KB,.0001,DT,ZZ,LAB) C------------------------------------------------------------------------ C DEFINE LATERAL B.C.'S C------------------------------------------------------------------------ CSB** c DO 48 I=1,IM c 48 COVRHN(I)=.1E0*SQRT(GRAV/H(I)) DO 50 I=1,IM L(I,1)=0.E0 50 L(I,KB)=0.E0 DO 51 I=1,IM UA(I)=UAB(I) VA(I)=VAB(I) EL(I)=ELB(I) ET(I)=ETB(I) ETF(I)=ET(I) D(I)=H(I)+EL(I) 51 DT(I)=H(I)+ET(I) DO 52 K=1,KB DO 52 I=1,IM L(I,K)=Q2LB(I,K)/(Q2B(I,K)+SMALL) KQ(I,K)=0.2E0*L(I,K)*SQRT(Q2B(I,K)) Q2(I,K)=Q2B(I,K) Q2L(I,K)=Q2LB(I,K) T(I,K)=TB(I,K) S(I,K)=SB(I,K) U(I,K)=UB(I,K) 52 V(I,K)=VB(I,K) C c DO 10 N=1,NREAD c READ(60) TIME0, c 1 WUBOT,WVBOT,AAM2D,UA,UAB,VA,VAB,EL,ELB,ET,ETB,EGB, c 2 UTB,VTB,U,UB,W,V,VB,T,TB,S,SB,RHO,ADVUU,ADVVV,ADVUA,ADVVA, c 3 KM,KH,KQ,L,Q2,Q2B,AAM,Q2L,Q2LB C READ(70) TIME0,UB,VB,UAB,VAB,TB,SB,ELB,Q2B,Q2LB,KM,KH 10 CONTINUE DO 53 I=1,IM D(I)=H(I)+EL(I) 53 DT(I)=H(I)+ET(I) C CALL VERTVL(DTI2) CALL BAROPG(DRHOX,TRNU) C CALL PRXZ(' DRHOX ',0.,DRHOX,IM,1,KB,0.,DT,ZZ) C CALL PRXZ(' TRNU ',0.,TRNU ,IM,1,1 ,0.,DT,ZZ) CALL ADVAVE(ADVUA,ADVVA,MODE) C C*********************************************************************** C * C BEGIN NUMERICAL INTEGRATION * C * C*********************************************************************** C DO 9000 INT=1,IEND C TIME=DAYI*DTI*FLOAT(INT)+TIME0 C IF(TIME.GT.3600.) CBC=0. RAMP=TIME/PERIOD IF(RAMP.GT.1.E0) RAMP=1.E0 FLOW=RAMP*FLOWB IF(RAMP.LT.0.99) THEN DO I=1,IM ELB(I)=ELBB(I)*RAMP ENDDO ENDIF C------------------------------------------------------------------------ C SET TIME DEPENDENT, SURFACE AND LATERAL BOUNDARY CONDITIONS. C THE LATTER WILL BE USED IN SUBROUTINE BCOND . C------------------------------------------------------------------------ DO 85 I=1,IM 85 WTSURF(I)=0.E0 C C------------------------------------------------------------------------ C IF(MODE.NE.2) THEN CALL BAROPG(DRHOX) ENDIF C********************************************************************** C HOR VISC = CONST*DX(I)*DY*SQRT((DU/DX(I))**2+(DV/DY)**2 C +.5*(DU/DY+DV/DX(I))**2) C********************************************************************** DO 95 K=1,KBM1 DO 95 I=2,IM-1 AAM(I,K)=HORCON*DX(I)**2 1 *SQRT( ((U(I+1,K)-U(I,K))/DX(I))**2 5 +.5E0*(.5E0*(V(I+1,K)-V(I-1,K)) 6 /DX(I)) **2) 95 CONTINUE DO 96 I=1,IM ADVUU(I)=0. ADVVV(I)=0. TRNU(I)=0. AAM2D(I)=0.E0 96 CONTINUE DO 199 K=1,KBM1 DO 199 I=1,IM ADVUU(I)=ADVUU(I)+ADVX(I,K)*DZ(K) ADVVV(I)=ADVVV(I)+ADVY(I,K)*DZ(K) TRNU(I)=TRNU(I)+DRHOX(I,K)*DZ(K) AAM2D(I)=AAM2D(I)+AAM(I,K)*DZ(K) 199 CONTINUE C CALL ADVAVE(ADVUA,ADVVA,MODE) DO 200 I=1,IM ADVUU(I)=ADVUU(I)-ADVUA(I) ADVVV(I)=ADVVV(I)-ADVVA(I) 200 CONTINUE DO 399 I=1,IM 399 EGF(I)=EL(I)*ISPI DO 999 I=2,IM UTF(I)=UA(I)*(D(I)+D(I-1))*ISP2I 999 VTF(I)=2.E0*VA(I)*D(I)*ISP2I C********** BEGIN EXTERNAL MODE *************************************** DO 8000 IEXT=1,ISPLIT DO 405 I=2,IM 405 FLUXUA(I)=.5*(D(I)+D(I-1))*UA(I) C DO 410 I=1,IMM1 410 ELF(I)=ELB(I) 1 -DTE2*(FLUXUA(I+1)-FLUXUA(I)) 2 /ART(I) C CALL BCOND(1) C IF(MOD(IEXT,ISPADV).EQ.0) CALL ADVAVE(ADVUA,ADVVA,MODE) C ALPHA=0.225 DO 420 I=3,IM 420 UAF(I)=ADVUU(I)+ADVUA(I) 1 -ARU(I)*(2.*CURV42D(I)*D(I)*VA(I) 2 +CURV42D(I-1)*D(I-1)*2.*VA(I-1) ) 3 +.5*GRAV*(D(I)+D(I-1)) 4 *( (1.-2.*ALPHA)*(EL(I)-EL(I-1)) 4 +ALPHA*(ELB(I)-ELB(I-1)+ELF(I)-ELF(I-1)) ) 5 +RAMP*TRNU(I) 6 -ARU(I)*(-.5*(WUSURF(I)+WUSURF(I-1))+WUBOT(I) ) DO 425 I=3,IM 425 UAF(I)= 1 ((H(I)+ELB(I)+H(I-1)+ELB(I-1))*ARU(I)*UAB(I) 2 -4.*DTE*UAF(I)) 3 /((H(I)+ELF(I)+H(I-1)+ELF(I-1))*ARU(I)) DO 430 I=1,IMM1 430 VAF(I)=ADVVV(I)+ADVVA(I) 1 +ARV(I)*( CURV42D(I)*D(I)*(UA(I+1)+UA(I)) 2 +CURV42D(I)*D(I)*(UA(I+1)+UA(I)) ) C 6 + ARV(I)*(WVSURF(I)-WVBOT(I) ) 6 + ARV(I)*(0.0 -WVBOT(I) ) DO 435 I=1,IM 435 VAF(I)= 1 ((H(I)+ELB(I))*VAB(I)*ARV(I) 2 -2.*DTE*VAF(I)) 3 /((H(I)+ELF(I))*ARV(I)) CALL BCOND(2) C IF(IEXT.LT.(ISPLIT-2)) GO TO 440 IF(IEXT.EQ.(ISPLIT-2))THEN DO 431 I=1,IM 431 ETF(I)=.25*SMOTH*ELF(I) ENDIF IF(IEXT.EQ.(ISPLIT-1)) THEN DO 432 I=1,IM 432 ETF(I)=ETF(I)+.5*(1.-.5*SMOTH)*ELF(I) ENDIF IF(IEXT.EQ.(ISPLIT-0)) THEN DO 433 I=1,IM 433 ETF(I)=(ETF(I)+.5*ELF(I))*FSM(I) ENDIF 440 CONTINUE C TEST FOR CFL VIOLATION. IF SO, PRINT AND STOP C VAMAX=-1.E10 DO 442 I=1,IM IF(ABS(VAF(I)).GE.VAMAX)VAMAX=VAF(I) 442 CONTINUE IF(ABS(VAMAX).GT.100.) GO TO 9001 C C SOLVE FOR ALONG-SHORE PRESSURE GRADIENT C TO MAINTAIN CONSTANT FLOW USING SPLIT TIME STEP TECHNIQUE C FLI=0. AREA=0. DO 446 I=1,IM FLI=FLI+VAF(I)*DX(I)*D(I) AREA=AREA+DX(I)*D(I) 446 CONTINUE DFL=(FLOW-FLI)/AREA DO I=1,IM VAF(I)=VAF(I)+DFL ENDDO PGRAD=-DFL/(GRAV*DTE2) C PGRAD=0. IF(INT.EQ.5) 1 write(6,'(1x,3E12.3)') FLOW,FLI,DFL c DO 447 I=1,IM c WVSURF(I)=-DFL*D(I)/DTE2 c IF(INT.EQ.5.and.IEXT.EQ.ISPLIT) c 1 write(6,'(1x,I5,2E12.3)') I,D(I),WVSURF(I) c 447 CONTINUE C C-------------------DIAGNOSTIC------------------------------------ IF(MOD(INT,IPRINT).EQ.0) THEN IF(IEXT.EQ.ISPLIT) THEN WRITE(6,'('' I ADVTR CORTR STRTR PGRTR 1 TENTR SUM ADVECY'')') C23456 | | | | | | | xxxxxxx| DO 500 I=2,IMM1 ADVTR(I)=(ADVVV(I)+ADVVA(I))/ARV(I) CORTR(I)= (CURV42D(I)*D(I)*(UA(I+1)+UA(I)) 2 +CURV42D(I)*D(I)*(UA(I+1)+UA(I))) STRTR(I)= (WVSURF(I)-WVBOT(I)) PGRTR(I)= D(I)*GRAV*PGRAD TENTR(I)= 1 ((H(I)+ELF(I))*VAF(I) 2 - (H(I)+ELB(I))*VAB(I))/DTE2 SUM=ADVTR(I)+CORTR(I)+STRTR(I)+PGRTR(I)+TENTR(I) WRITE(6,'(1x,I5,7E12.3)') 1 I,ADVTR(I),CORTR(I),STRTR(I),PGRTR(I),TENTR(I),SUM,ADVECY(I) 500 CONTINUE ENDIF ENDIF C------------------------------------------------------------------- C APPLY FILTER TO REMOVE TIME SPIMT AND RESET TIME SEQUENCE DO 445 I=1,IM UA(I)=UA(I)+.5E0*SMOTH*(UAB(I)-2.E0*UA(I)+UAF(I)) VA(I)=VA(I)+.5E0*SMOTH*(VAB(I)-2.E0*VA(I)+VAF(I)) EL(I)=EL(I)+.5E0*SMOTH*(ELB(I)-2.E0*EL(I)+ELF(I)) ELB(I)=EL(I) EL(I)=ELF(I) D(I)=H(I)+EL(I) UAB(I)=UA(I) UA(I)=UAF(I) VAB(I)=VA(I) VA(I)=VAF(I) 445 CONTINUE C IF(IEXT.EQ.ISPLIT) GO TO 8000 DO 450 I=1,IM EGF(I)=EGF(I)+EL(I)*ISPI UTF(I)=UTF(I)+UA(I)*(D(I)+D(I-1))*ISP2I 450 VTF(I)=VTF(I)+VA(I)*(D(I)+D(I))*ISP2I C C 8000 CONTINUE C--------------------------------------------------------------------- C END EXTERNAL (2-D) MODE CALCULATION C AND CONTINUE WITH INTERNAL (3-D) MODE CALCULATION C--------------------------------------------------------------------- IF(INT.EQ.1) GO TO 8200 IF(MODE.EQ.2) GO TO 8200 DO 777 I=1,IM 777 TPS(I)=0.E0 DO 299 K=1,KBM1 DO 299 I=1,IM 299 TPS(I)=TPS(I)+U(I,K)*DZ(K) DO 302 K=1,KBM1 DO 302 I=2,IMM1 302 U(I,K)=(U(I,K)-TPS(I))+ 1 (UTB(I)+UTF(I))/(DT(I)+DT(I-1)) DO 3021 I=1,IM 3021 TPS(I)=0.E0 DO 304 K=1,KBM1 DO 304 I=2,IMM1 304 TPS(I)=TPS(I)+V(I,K)*DZ(K) DO 306 K=1,KBM1 DO 306 I=2,IMM1 306 V(I,K)=(V(I,K)-TPS(I))+ 2 (VTB(I)+VTF(I))/(2.E0*DT(I)) C C---------------------------------------------------------------- C VERTVL INPUT = U,V,DT(=H+ET),ETF,ETB; OUTPUT = W C---------------------------------------------------------------- CALL VERTVL(DTI2) CALL BCOND(5) C CSB** DO 307 K=1,KB DO 307 I=1,IM UF(I,K)=0.E0 307 VF(I,K)=0.E0 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 367 I=1,IM DO 367 K=1,KB Q2(I,K)=Q2(I,K)+.5*SMOTH*(UF(I,K)+Q2B(I,K)-2.E0*Q2(I,K)) Q2L(I,K)=Q2L(I,K)+.5*SMOTH*(VF(I,K)+Q2LB(I,K)-2.E0*Q2L(I,K)) Q2B(I,K)=Q2(I,K) Q2(I,K)=UF(I,K) Q2LB(I,K)=Q2L(I,K) Q2L(I,K)=VF(I,K) 367 CONTINUE C---------------------------------------------------------------- C COMPUTE TF AN SF USING UF AND VF AS TEMPORARY VARIABLES C---------------------------------------------------------------- IF(MODE.EQ.4) GO TO 360 CALL ADVT(TB,T,TCLIM,DTI2,UF) C CALL ADVT(SB,S,SCLIM,DTI2,VF) CALL PROFT(UF,WTSURF,TSURF,1,DTI2) C CALL PROFT(VF,WSSURF,SSURF,1,DTI2) CALL BCOND(4) C DO 355 K=1,KBM1 DO 355 I=1,IM T(I,K)=T(I,K)+.5*SMOTH*(UF(I,K) 1 +TB(I,K)-2.E0*T(I,K)) S(I,K)=S(I,K)+.5*SMOTH*(VF(I,K) 1 +SB(I,K)-2.E0*S(I,K)) TB(I,K)=T(I,K) T(I,K)=UF(I,K) SB(I,K)=S(I,K) S(I,K)=VF(I,K) 355 RHO(I,K)=T(I,K) C C CALL DENS C 360 CONTINUE DO 678 I=1,IM 678 EGB(I)=EGF(I) C---------------------------------------------------------------- C COMPUTE UF AND VF C---------------------------------------------------------------- CALL ADVU(DRHOX,ADVX,DTI2) CALL ADVV(PGRAD,ADVY,DTI2,ADVECY) CALL PROFU(DTI2) CALL PROFV(DTI2) CALL BCOND(3) C CSB** DO 369 I=1,IM 369 TPS(I)=0.E0 DO 370 K=1,KBM1 DO 370 I=1,IM 370 TPS(I)=TPS(I)+(UF(I,K)+UB(I,K)-2.E0*U(I,K))*DZ(K) DO 372 K=1,KBM1 DO 372 I=1,IM 372 U(I,K)=U(I,K)+.5*SMOTH*(UF(I,K)+UB(I,K)-2.E0*U(I,K) 1 -TPS(I)) DO 3721 I=1,IM 3721 TPS(I)=0.E0 DO 374 K=1,KBM1 DO 374 I=1,IM 374 TPS(I)=TPS(I)+(VF(I,K)+VB(I,K)-2.E0*V(I,K))*DZ(K) DO 376 K=1,KBM1 DO 376 I=1,IM 376 V(I,K)=V(I,K)+.5*SMOTH*(VF(I,K)+VB(I,K)-2.E0*V(I,K) 1 -TPS(I)) C DO 377 K=1,KB DO 377 I=1,IM UB(I,K)=U(I,K) U(I,K)=UF(I,K) VB(I,K)=V(I,K) 377 V(I,K)=VF(I,K) C 8200 CONTINUE C DO 8210 I=1,IM EGB(I)=EGF(I) ETB(I)=ET(I) ET(I)=ETF(I) DT(I)=H(I)+ET(I) UTB(I)=UTF(I) 8210 VTB(I)=VTF(I) C-------------------------------------------------------------------- C BEGIN PRINT SECTION C-------------------------------------------------------------------- IF(MOD(INT,IPRINT/5).EQ.0) THEN c WRITE(6,'('' TIME, PGRAD,WVBOT(21) ='',F10.2,2E12.3)') c 1 TIME,PGRAD,WVBOT(21) WRITE(6,'('' TIME, WVSURF(21),WVBOT(21) ='',F10.2,2E12.3)') 1 TIME,WVSURF(21),WVBOT(21) ENDIF C IF(INT.GE.ISWTCH) IPRINT=PRTD2*24*3600/IFIX(DTI) IF(MOD(INT,IPRINT).NE.0) GO TO 7000 C 9001 CONTINUE C c WRITE(6,'(/,'' TIME ='',F10.2,'' INT ='',I8, c 1 '' IEXT ='',I8,'' IPRINT ='',I5,'' MODE ='',I5,//)') c 2 TIME,INT,IEXT,IPRINT,MODE C WRITE(6,'('' TIME ='',F10.2)') TIME C CALL PRXZSH(' DRHOX ',TIME,DRHOX ,IM,1,KB,0.,DT,ZZ) C CALL PRXZ(' TRNU ',TIME,TRNU,IM,1,1,0.,DT,0.) c CALL PRXZSH(' AAM ',TIME,AAM,IM,1,KB,0.,DT,ZZ,1) LAB=1 CALL PRXZ(' ELB ',TIME,ELB,IM,1,1,0.,DT,0.1,LAB) CALL PRXZ(' UAB ',TIME,UAB,IM,1,1,0.,DT,0.,LAB) CALL PRXZ(' VAB ',TIME,VAB,IM,1,1,0.,DT,0.,LAB) WRITE(6,'('' DFL ='',E12.3,'' PGRAD ='',E12.3)') DFL,PGRAD CALL PRXZ('WVSURF',TIME,WVSURF,IM,1,1 ,0. ,DT,ZZ,LAB) CALL PRXZ('WUBOT',TIME,WUBOT,IM,1,1 ,0. ,DT,ZZ,LAB) CALL PRXZ('WVBOT',TIME,WVBOT,IM,1,1 ,0. ,DT,ZZ,LAB) C CALL PRXZ('ADVVA ',TIME, ADVVA ,IM,1,1,0.,DT,Z) C CALL PRXZ('ADVVV ',TIME, ADVVV ,IM,1,1,0.,DT,Z) CALL PRXZ(' VB ',TIME,VB ,IM,1,KB,1.E-3,DT,ZZ,LAB) CALL PRXZ(' UB ',TIME,UB ,IM,1,KB,1.E-3,DT,ZZ,LAB) c DO I=1,IM c DO K=1,KBM1 c RHO(I,K)=(RHO(I,K)-RHOBOT)*FSM(I) c ENDDO c ENDDO CALL PRXZ(' RHO ',TIME,RHO,IM,1,KB,1.E-4,DT,ZZ,LAB) c DO I=1,IM c DO K=1,KBM1 c RHO(I,K)=(RHO(I,K)+RHOBOT)*FSM(I) c ENDDO c ENDDO CALL PRXZ(' KM ',TIME, KM,IM,1,KB,1.E-3,DT,Z,LAB) C DO I=1,IM PSI(I,1)=0 DO K=2,KB PSI(I,K)=PSI(I,K-1)-U(I,K)*DZ(K)*0.5*(DT(I)+DT(I-1)) ENDDO ENDDO CALL PRXZ(' PSI ',TIME,PSI,IM,1,KB,0. ,DT,ZZ,LAB) WRITE(6,'('' SAVE TO TAPE. TIME ='',F10.2)') TIME WRITE(80) TIME,DFL,Z,ZZ,H,ELB,WVBOT,UB,VB,UAB,VAB,TB,SB,ELB, 1 PSI,KH,KM 7000 CONTINUE C-------------------------------------------------------------------- C END PRINT SECTION C-------------------------------------------------------------------- IF(ABS(VAMAX).GT.100.) THEN WRITE(6,'(///,''ABNORMAL END'')') STOP ENDIF 9000 CONTINUE C CALL PRXZ(' W ',TIME,W ,IM,1,KB,0.,DT,Z) C C REWIND(70) C WRITE(70) TIME, C 1 WUBOT,WVBOT,AAM2D,UA,UAB,VA,VAB,EL,ELB,ET,ETB,EGB, C 2 UTB,VTB,U,UB,W,V,VB,T,TB,S,SB,RHO,ADVUU,ADVVV,ADVUA,ADVVA, C 3 KM,KH,KQ,L,Q2,Q2B,AAM,Q2L,Q2LB C WRITE(6,'('' RUN COMPLETED, TIME = '', F10.2)') TIME STOP END C ########################################################### SUBROUTINE ADVAVE(ADVUA,ADVVA,MODE) C REAL KM,KH,KQ,L PARAMETER (IM=41,KB=30,KS=20,IMM1=IM-1,KBM1=KB-1) PARAMETER (IMM2=IM-2,KBM2=KB-2) PARAMETER (IM2=IM*2,IMK=IM*KB,IMKM2=IM*KBM2) PARAMETER (IMKM1=IM*KBM1) COMMON/BLKCON/ 1 INT,IPRINT,DTE,DTI,TPRNU,UMOL,PR, 2 GRAV,TIME,RAMP C---------------- 1-D ARRAYS -------------------------------------- COMMON/BLK1D/ 1 DZR(KB),Z(KB),ZZ(KB),DZ(KB),DZZ(KB) C---------------- 2-D ARRAYS -------------------------------------- COMMON/BLK2D/H(IM),DX(IM),D(IM),DT(IM), 1 ART(IM),ARU(IM),ARV(IM),CBC(IM), 2 ALON(IM),ALAT(IM),ANG(IM), 3 DUM(IM),DVM(IM),FSM(IM),COR4(IM),CURV42D(IM), 4 WUSURF(IM),WVSURF(IM),WUBOT(IM),WVBOT(IM), 5 WTSURF(IM),WSSURF(IM),TPS(IM),AAM2D(IM), 6 UAF(IM),UA(IM),UAB(IM),VAF(IM),VA(IM), 7 VAB(IM),ELF(IM),EL(IM),ELB(IM), 8 ETF(IM),ET(IM),ETB(IM),FLUXUA(IM),FLUXVA(IM), 9 EGF(IM),EGB(IM), ELBB(IM) C---------------- 3-D ARRAYS -------------------------------------- COMMON/BLK3D/ 1 A(IM,KB),C(IM,KB),VH(IM,KB),VHP(IM,KB), 1 UF(IM,KB),VF(IM,KB), 2 KM(IM,KB),KH(IM,KB),KQ(IM,KB),L(IM,KB), 3 Q2(IM,KB),Q2B(IM,KB),AAM(IM,KB), 4 Q2L(IM,KB),Q2LB(IM,KB), 5 U(IM,KB),UB(IM,KB),W(IM,KB), 6 V(IM,KB),VB(IM,KB), 7 T(IM,KB),TB(IM,KB), 8 S(IM,KB),SB(IM,KB), 9 RHO(IM,KB),DTEF(IM,KB),RMEAN(IM,KB) C----------- 1 AND 2-D BOUNDARY VALUE ARRAYS ------------------------ COMMON/BDRY/ 1 TBE(KB),TBN(IM,KB),TBS(IM,KB),SBN(IM,KB),SBE(KB), 2 SBS(IM,KB),VABN(IM),VABS(IM),UGS(KB),VGS(IM,KB), 3 ELN(IM),ELS(IM),COVRHN(IM) DIMENSION ADVUA(IM),ADVVA(IM) C--------------------------------------------------------------------- C CALCULATE U-ADVECTION & DIFFUSION C--------------------------------------------------------------------- C C-------- ADVECTIVE FLUXES ------------------------------------------- DO 299 I=1,IM 299 TPS(I)=0.E0 DO 300 I=2,IMM1 300 FLUXUA(I)=.125E0*((D(I+1)+D(I))*UA(I+1) 1 +(D(I)+D(I-1))*UA(I)) 2 *(UA(I+1)+UA(I)) C----------- ADD VISCOUS FLUXES --------------------------------- DO 460 I=1,IMM1 460 FLUXUA(I)=FLUXUA(I) 1 -D(I)*2.E0*AAM2D(I)*(UAB(I+1)-UAB(I))/DX(I) DO 470 I=2,IM TPS(I)=.5*(D(I)+D(I-1))*(AAM2D(I)+AAM2D(I-1)) 4 *(VAB(I)-VAB(I-1)) 5 /(DX(I)+DX(I-1)) 470 CONTINUE C---------------------------------------------------------------- DO 480 I=2,IM 480 ADVUA(I)=FLUXUA(I)-FLUXUA(I-1) C---------------------------------------------------------------- C CALCULATE V-ADVECTION & DIFFUSION C---------------------------------------------------------------- DO 481 I=1,IM 481 ADVVA(I)=0.0E0 C---------ADVECTIVE FLUXES ---------------------------- DO 700 I=2,IM 700 FLUXUA(I)=.25E0*((D(I)+D(I-1))*UA(I)) 2 *(VA(I-1)+VA(I)) C------- ADD VISCOUS FLUXES ----------------------------------- DO 870 I=1,IM 870 FLUXUA(I)=FLUXUA(I)-TPS(I) C--------------------------------------------------------------- DO 880 I=1,IMM1 880 ADVVA(I)=FLUXUA(I+1)-FLUXUA(I) C--------------------------------------------------------------- C COMPUTE BOTTOM FRICTION IF RUN IS EXTERNAL MODE ONLY C--------------------------------------------------------------- IF(MODE.NE.2) GO TO 5000 DO 100 I=2,IM WUBOT(I)=0.E0 WUBOT(I)=-0.5E0*(CBC(I)+CBC(I-1)) 1 *SQRT(UAB(I)**2+(.5E0*(VAB(I) 2 +VAB(I-1)))**2)*UAB(I) 100 CONTINUE DO 102 I=1,IMM1 WVBOT(I)=0.E0 WVBOT(I)=-1.E0*(CBC(I)) 1 *SQRT((.5E0*(UAB(I)+UAB(I+1)))**2 2 +VAB(I)**2)*VAB(I) 102 CONTINUE 5000 CONTINUE RETURN END C############################################################# SUBROUTINE ADVQ(QB,Q,DTI2,QF) C REAL KM,KH,KQ,L PARAMETER (IM=41,KB=30,KS=20,IMM1=IM-1,KBM1=KB-1) PARAMETER (IMM2=IM-2,KBM2=KB-2) PARAMETER (IM2=IM*2,IMK=IM*KB,IMKM2=IM*KBM2) PARAMETER (IMKM1=IM*KBM1) COMMON/BLKCON/ 1 INT,IPRINT,DTE,DTI,TPRNU,UMOL,PR, 2 GRAV,TIME,RAMP C---------------- 1-D ARRAYS -------------------------------------- COMMON/BLK1D/ 1 DZR(KB),Z(KB),ZZ(KB),DZ(KB),DZZ(KB) C---------------- 2-D ARRAYS -------------------------------------- COMMON/BLK2D/H(IM),DX(IM),D(IM),DT(IM), 1 ART(IM),ARU(IM),ARV(IM),CBC(IM), 2 ALON(IM),ALAT(IM),ANG(IM), 3 DUM(IM),DVM(IM),FSM(IM),COR4(IM),CURV42D(IM), 4 WUSURF(IM),WVSURF(IM),WUBOT(IM),WVBOT(IM), 5 WTSURF(IM),WSSURF(IM),TPS(IM),AAM2D(IM), 6 UAF(IM),UA(IM),UAB(IM),VAF(IM),VA(IM), 7 VAB(IM),ELF(IM),EL(IM),ELB(IM), 8 ETF(IM),ET(IM),ETB(IM),FLUXUA(IM),FLUXVA(IM), 9 EGF(IM),EGB(IM), ELBB(IM) C---------------- 3-D ARRAYS -------------------------------------- COMMON/BLK3D/ 1 A(IM,KB),C(IM,KB),VH(IM,KB),VHP(IM,KB), 1 UF(IM,KB),VF(IM,KB), 2 KM(IM,KB),KH(IM,KB),KQ(IM,KB),L(IM,KB), 3 Q2(IM,KB),Q2B(IM,KB),AAM(IM,KB), 4 Q2L(IM,KB),Q2LB(IM,KB), 5 U(IM,KB),UB(IM,KB),W(IM,KB), 6 V(IM,KB),VB(IM,KB), 7 T(IM,KB),TB(IM,KB), 8 S(IM,KB),SB(IM,KB), 9 RHO(IM,KB),DTEF(IM,KB),RMEAN(IM,KB) C----------- 1 AND 2-D BOUNDARY VALUE ARRAYS ------------------------ COMMON/BDRY/ 1 TBE(KB),TBN(IM,KB),TBS(IM,KB),SBN(IM,KB),SBE(KB), 2 SBS(IM,KB),VABN(IM),VABS(IM),UGS(KB),VGS(IM,KB), 3 ELN(IM),ELS(IM),COVRHN(IM) DIMENSION QB(IM,KB),Q(IM,KB),QF(IM,KB) 3 DIMENSION XFLUX(IM,KB),YFLUX(IM,KB) EQUIVALENCE (XFLUX,A),(YFLUX,C) C C******* HORIZONTAL ADVECTION ************************************ DO 110 K=2,KBM1 DO 110 I=2,IM 110 XFLUX(I,K)=.125E0*(Q(I,K)+Q(I-1,K)) 1 *(DT(I)+DT(I-1))*(U(I,K)+U(I,K-1)) C******* HORIZONTAL DIFFUSION ************************************ DO 315 K=2,KBM1 DO 315 I=2,IM XFLUX(I,K)=XFLUX(I,K) 1 -.5E0*(AAM(I,K)+AAM(I-1,K))*(H(I)+H(I-1)) 2 *(QB(I,K)-QB(I-1,K))*DUM(I) 3 /(DX(I)+DX(I-1)) 315 CONTINUE C****** VERTICAL ADVECTION; ADD FLUX TERMS ;THEN STEP FORWARD IN TIME ** DO 230 K=2,KBM1 DO 230 I=1,IMM1 QF(I,K)=(W(I,K-1)*Q(I,K-1)-W(I,K+1)*Q(I,K+1)) 1 /(DZ(K)+DZ(K-1))*ART(I) 2 +XFLUX(I+1,K)-XFLUX(I,K) 230 QF(I,K)=((H(I)+ETB(I))*ART(I)*QB(I,K)-DTI2*QF(I,K)) 1 /((H(I)+ETF(I))*ART(I)) RETURN END C############################################################# SUBROUTINE ADVT(FB,F,FCLIM,DTI2,FF) C C THIS SUBROUTINE INTEGRATES CONSERVATIVE CONSTITUENT EQUATIONS C PASSIVE AND ACTIVE...... C REAL KM,KH,KQ,L PARAMETER (IM=41,KB=30,KS=20,IMM1=IM-1,KBM1=KB-1) PARAMETER (IMM2=IM-2,KBM2=KB-2) PARAMETER (IM2=IM*2,IMK=IM*KB,IMKM2=IM*KBM2) PARAMETER (IMKM1=IM*KBM1) COMMON/BLKCON/ 1 INT,IPRINT,DTE,DTI,TPRNU,UMOL,PR, 2 GRAV,TIME,RAMP C---------------- 1-D ARRAYS -------------------------------------- COMMON/BLK1D/ 1 DZR(KB),Z(KB),ZZ(KB),DZ(KB),DZZ(KB) C---------------- 2-D ARRAYS -------------------------------------- COMMON/BLK2D/H(IM),DX(IM),D(IM),DT(IM), 1 ART(IM),ARU(IM),ARV(IM),CBC(IM), 2 ALON(IM),ALAT(IM),ANG(IM), 3 DUM(IM),DVM(IM),FSM(IM),COR4(IM),CURV42D(IM), 4 WUSURF(IM),WVSURF(IM),WUBOT(IM),WVBOT(IM), 5 WTSURF(IM),WSSURF(IM),TPS(IM),AAM2D(IM), 6 UAF(IM),UA(IM),UAB(IM),VAF(IM),VA(IM), 7 VAB(IM),ELF(IM),EL(IM),ELB(IM), 8 ETF(IM),ET(IM),ETB(IM),FLUXUA(IM),FLUXVA(IM), 9 EGF(IM),EGB(IM), ELBB(IM) C---------------- 3-D ARRAYS -------------------------------------- COMMON/BLK3D/ 1 A(IM,KB),C(IM,KB),VH(IM,KB),VHP(IM,KB), 1 UF(IM,KB),VF(IM,KB), 2 KM(IM,KB),KH(IM,KB),KQ(IM,KB),L(IM,KB), 3 Q2(IM,KB),Q2B(IM,KB),AAM(IM,KB), 4 Q2L(IM,KB),Q2LB(IM,KB), 5 U(IM,KB),UB(IM,KB),W(IM,KB), 6 V(IM,KB),VB(IM,KB), 7 T(IM,KB),TB(IM,KB), 8 S(IM,KB),SB(IM,KB), 9 RHO(IM,KB),DTEF(IM,KB),RMEAN(IM,KB) C----------- 1 AND 2-D BOUNDARY VALUE ARRAYS ------------------------ COMMON/BDRY/ 1 TBE(KB),TBN(IM,KB),TBS(IM,KB),SBN(IM,KB),SBE(KB), 2 SBS(IM,KB),VABN(IM),VABS(IM),UGS(KB),VGS(IM,KB), 3 ELN(IM),ELS(IM),COVRHN(IM) DIMENSION FB(IM,KB),F(IM,KB),FF(IM,KB),FCLIM(IM,KB) DIMENSION XFLUX(IM,KB),YFLUX(IM,KB) EQUIVALENCE (XFLUX,A),(YFLUX,C) DO K=1,KB DO I=1,IM XFLUX(I,K)=0. ENDDO ENDDO DO 529 I=1,IM F(I,KB)=F(I,KBM1) 529 FB(I,KB)=FB(I,KBM1) C******* DO ADVECTION FLUXES ************************************** DO 530 K=1,KBM1 DO 530 I=2,IM 530 XFLUX(I,K)=.25*((DT(I)+DT(I-1)) 2 *(F(I,K)+F(I-1,K))*U(I,K)) C****** ADD DIFFUSIVE FLUXES ************************************* C********* DELETE DIFFUSIVE FLUXES ******************************* IF(2..GT.1.) GOTO 300 DO 99 K=1,KB DO 99 I=1,IM 99 FB(I,K)=FB(I,K)-FCLIM(I,K) DO 100 K=1,KB-1 DO 100 I=2,IM XFLUX(I,K)=XFLUX(I,K) 1 -.5*(AAM(I,K)+AAM(I-1,K)) 2 *(H(I)+H(I-1))*( (FB(I,K)-FB(I-1,K)) 5 )*DUM(I)/(TPRNU*(DX(I)+DX(I-1))) 100 CONTINUE DO 102 K=1,KB DO 102 I=1,IM 102 FB(I,K)=FB(I,K)+FCLIM(I,K) 300 CONTINUE C********* DELETE DIFFUSIVE FLUXES ******************************* C****** DO VERTICAL ADVECTION ************************************* DO 505 I=1,IM 505 FF(I,1)=-.5E0*DZR(1)*(F(I,1)+F(I,2))*W(I,2)*ART(I) DO 520 K=2,KBM1 DO 520 I=1,IM 520 FF(I,K)=.5E0*DZR(K)*((F(I,K-1)+F(I,K))*W(I,K) 1 -(F(I,K)+F(I,K+1))*W(I,K+1))*ART(I) C****** ADD NET HORIZONTAL FLUXES; THEN STEP FORWARD IN TIME ********** DO 120 K=1,KBM1 DO 120 I=1,IM FF(I,K)=FF(I,K) 1 +XFLUX(I+1,K)-XFLUX(I,K) FF(I,K)=(FB(I,K)*(H(I)+ETB(I))*ART(I)-DTI2*FF(I,K)) 1 /((H(I)+ETF(I))*DX(I)) 120 CONTINUE RETURN END C################################################################# SUBROUTINE ADVU(DRHOX,ADVX,DTI2) C REAL KM,KH,KQ,L PARAMETER (IM=41,KB=30,KS=20,IMM1=IM-1,KBM1=KB-1) PARAMETER (IMM2=IM-2,KBM2=KB-2) PARAMETER (IM2=IM*2,IMK=IM*KB,IMKM2=IM*KBM2) PARAMETER (IMKM1=IM*KBM1) COMMON/BLKCON/ 1 INT,IPRINT,DTE,DTI,TPRNU,UMOL,PR, 2 GRAV,TIME,RAMP C---------------- 1-D ARRAYS -------------------------------------- COMMON/BLK1D/ 1 DZR(KB),Z(KB),ZZ(KB),DZ(KB),DZZ(KB) C---------------- 2-D ARRAYS -------------------------------------- COMMON/BLK2D/H(IM),DX(IM),D(IM),DT(IM), 1 ART(IM),ARU(IM),ARV(IM),CBC(IM), 2 ALON(IM),ALAT(IM),ANG(IM), 3 DUM(IM),DVM(IM),FSM(IM),COR4(IM),CURV42D(IM), 4 WUSURF(IM),WVSURF(IM),WUBOT(IM),WVBOT(IM), 5 WTSURF(IM),WSSURF(IM),TPS(IM),AAM2D(IM), 6 UAF(IM),UA(IM),UAB(IM),VAF(IM),VA(IM), 7 VAB(IM),ELF(IM),EL(IM),ELB(IM), 8 ETF(IM),ET(IM),ETB(IM),FLUXUA(IM),FLUXVA(IM), 9 EGF(IM),EGB(IM), ELBB(IM) C---------------- 3-D ARRAYS -------------------------------------- COMMON/BLK3D/ 1 A(IM,KB),C(IM,KB),VH(IM,KB),VHP(IM,KB), 1 UF(IM,KB),VF(IM,KB), 2 KM(IM,KB),KH(IM,KB),KQ(IM,KB),L(IM,KB), 3 Q2(IM,KB),Q2B(IM,KB),AAM(IM,KB), 4 Q2L(IM,KB),Q2LB(IM,KB), 5 U(IM,KB),UB(IM,KB),W(IM,KB), 6 V(IM,KB),VB(IM,KB), 7 T(IM,KB),TB(IM,KB), 8 S(IM,KB),SB(IM,KB), 9 RHO(IM,KB),DTEF(IM,KB),RMEAN(IM,KB) C----------- 1 AND 2-D BOUNDARY VALUE ARRAYS ------------------------ COMMON/BDRY/ 1 TBE(KB),TBN(IM,KB),TBS(IM,KB),SBN(IM,KB),SBE(KB), 2 SBS(IM,KB),VABN(IM),VABS(IM),UGS(KB),VGS(IM,KB), 3 ELN(IM),ELS(IM),COVRHN(IM) DIMENSION DRHOX(IM,KB),XFLUX(IM,KB),YFLUX(IM,KB), 1 ADVX(IM,KB),CURV4(IM,KB) EQUIVALENCE (A,XFLUX),(C,YFLUX),(VH,CURV4) C DO 59 K=1,KB DO 59 I=1,IM UF(I,K)=0.E0 59 XFLUX(I,K)=0.E0 DO 60 K=1,KBM1 DO 60 I=1,IM 60 CURV4(I,K)=COR4(I) DO 61 I=1,IM 61 CURV42D(I)=0.E0 DO 70 K=1,KBM1 DO 70 I=1,IM 70 CURV42D(I)=CURV42D(I)+CURV4(I,K)*DZ(K) C******** HORIZONTAL ADVECTION ************************************ DO 100 K=1,KBM1 DO 100 I=2,IMM1 100 XFLUX(I,K)=.125E0*((DT(I+1)+DT(I))* 1 U(I+1,K)+(DT(I)+DT(I-1)) 2 *U(I,K))*(U(I+1,K)+U(I,K)) C****** HORIZONTAL DIFFUSION ************************************* DO 700 K=1,KBM1 DO 700 I=2,IMM1 XFLUX(I,K)=XFLUX(I,K) 1 -DT(I)*AAM(I,K)*2.E0*(UB(I+1,K)-UB(I,K))/DX(I) 700 CONTINUE C******** DO VERT. ADVECTION; ADD HORIZ. ADVECTION ******* CSB** DO 139 I=1,IM UF(I,1)=0.E0 139 UF(I,KB)=0.E0 DO 140 K=2,KBM1 DO 140 I=2,IM 140 UF(I,K)=.25E0*(W(I,K)+W(I-1,K))*(U(I,K)+U(I,K-1)) DO 145 K=1,KBM1 DO 145 I=1,IM 145 ADVX(I,K)=DZR(K)*(UF(I,K)-UF(I,K+1))*ARU(I) 1 +XFLUX(I,K)-XFLUX(I-1,K) C-------------------------------------------------------------------- C-------------------------------------------------------------------- C********* -FVD + GDEG/DX + BAROCLINIC TERM ********************** DO 150 K=1,KBM1 DO 150 I=2,IM 150 UF(I,K)=ADVX(I,K) 1 -ARU(I)*(CURV4(I,K)*DT(I)*2.*V(I,K) 2 +CURV4(I-1,K)*DT(I-1)*2.*V(I-1,K)) 3 +GRAV*.25E0*(DT(I)+DT(I-1)) 4 *(EGF(I)-EGF(I-1)+EGB(I)-EGB(I-1)) 6 +DRHOX(I,K) C******* STEP FORWARD IN TIME *********************************** DO 190 K=1,KBM1 DO 190 I=2,IM 190 UF(I,K)= 1 ((H(I)+ETB(I)+H(I-1)+ETB(I-1))*ARU(I)*UB(I,K) 2 -2.E0*DTI2*UF(I,K)) 3 /((H(I)+ETF(I)+H(I-1)+ETF(I-1))*ARU(I)) RETURN END C################################################################# SUBROUTINE ADVV(PGRAD,ADVY,DTI2,ADVECY) C REAL KM,KH,KQ,L PARAMETER (IM=41,KB=30,KS=20,IMM1=IM-1,KBM1=KB-1) PARAMETER (IMM2=IM-2,KBM2=KB-2) PARAMETER (IM2=IM*2,IMK=IM*KB,IMKM2=IM*KBM2) PARAMETER (IMKM1=IM*KBM1) COMMON/BLKCON/ 1 INT,IPRINT,DTE,DTI,TPRNU,UMOL,PR, 2 GRAV,TIME,RAMP C---------------- 1-D ARRAYS -------------------------------------- COMMON/BLK1D/ 1 DZR(KB),Z(KB),ZZ(KB),DZ(KB),DZZ(KB) C---------------- 2-D ARRAYS -------------------------------------- COMMON/BLK2D/H(IM),DX(IM),D(IM),DT(IM), 1 ART(IM),ARU(IM),ARV(IM),CBC(IM), 2 ALON(IM),ALAT(IM),ANG(IM), 3 DUM(IM),DVM(IM),FSM(IM),COR4(IM),CURV42D(IM), 4 WUSURF(IM),WVSURF(IM),WUBOT(IM),WVBOT(IM), 5 WTSURF(IM),WSSURF(IM),TPS(IM),AAM2D(IM), 6 UAF(IM),UA(IM),UAB(IM),VAF(IM),VA(IM), 7 VAB(IM),ELF(IM),EL(IM),ELB(IM), 8 ETF(IM),ET(IM),ETB(IM),FLUXUA(IM),FLUXVA(IM), 9 EGF(IM),EGB(IM), ELBB(IM) C---------------- 3-D ARRAYS -------------------------------------- COMMON/BLK3D/ 1 A(IM,KB),C(IM,KB),VH(IM,KB),VHP(IM,KB), 1 UF(IM,KB),VF(IM,KB), 2 KM(IM,KB),KH(IM,KB),KQ(IM,KB),L(IM,KB), 3 Q2(IM,KB),Q2B(IM,KB),AAM(IM,KB), 4 Q2L(IM,KB),Q2LB(IM,KB), 5 U(IM,KB),UB(IM,KB),W(IM,KB), 6 V(IM,KB),VB(IM,KB), 7 T(IM,KB),TB(IM,KB), 8 S(IM,KB),SB(IM,KB), 9 RHO(IM,KB),DTEF(IM,KB),RMEAN(IM,KB) C----------- 1 AND 2-D BOUNDARY VALUE ARRAYS ------------------------ COMMON/BDRY/ 1 TBE(KB),TBN(IM,KB),TBS(IM,KB),SBN(IM,KB),SBE(KB), 2 SBS(IM,KB),VABN(IM),VABS(IM),UGS(KB),VGS(IM,KB), 3 ELN(IM),ELS(IM),COVRHN(IM) DIMENSION XFLUX(IM,KB),YFLUX(IM,KB), 1 ADVY(IM,KB),CURV4(IM,KB) DIMENSION ADVECY(IM) EQUIVALENCE (A,XFLUX),(C,YFLUX),(VH,CURV4) C DO 299 K=1,KB DO 299 I=1,IM VF(I,K)=0.E0 299 XFLUX(I,K)=0.E0 C C********** HORIZONTAL ADVECTION ********************************* DO 300 K=1,KBM1 DO 300 I=2,IM 300 XFLUX(I,K)=.25E0*((DT(I)+DT(I-1))*U(I,K)) 2 *(V(I,K)+V(I-1,K)) DO I=1,IM ADVECY(I)=0 ENDDO DO K=1,KBM1 DO I=1,IMM1 ADVECY(I)=ADVECY(I)+(DZ(K)*(XFLUX(I+1,K)-XFLUX(I,K)))/DX(I) ENDDO ENDDO C******* HORIZONTAL DIFFUSION ************************************* DO 700 K=1,KBM1 DO 700 I=2,IM DTAAM=(DT(I)+DT(I-1)) 1 *(AAM(I,K)+AAM(I-1,K)) XFLUX(I,K)=XFLUX(I,K) 1 -.5E0*DTAAM* 3 (VB(I,K)-VB(I-1,K)) 4 /(DX(I)+DX(I-1)) 700 CONTINUE C C********** DO VERT. ADVECTION; ADD HORIZ. ADVECTION ************ DO 359 I=1,IM VF(I,1)=0. 359 VF(I,KB)=0.E0 DO 360 K=2,KBM1 DO 360 I=1,IM 360 VF(I,K)=.5E0*W(I,K)*(V(I,K)+V(I,K-1)) DO 400 K=1,KBM1 DO 400 I=1,IMM1 400 ADVY(I,K)=DZR(K)*(VF(I,K)-VF(I,K+1))*ARV(I) 1 +XFLUX(I+1,K)-XFLUX(I,K) C-------------------------------------------------------------------- C-------------------------------------------------------------------- C********* +FUD + GDEG/DX + BAROCIMNIC TERM ************************* DO 340 K=1,KBM1 DO 340 I=1,IMM1 340 VF(I,K)=ADVY(I,K) 1 +ARV(I)*(CURV4(I,K)*DT(I)*(U(I+1,K)+U(I,K)) 2 +CURV4(I,K)*DT(I)*(U(I+1,K)+U(I,K))) 4 +GRAV*ARV(I)*DT(I)*PGRAD C******* STEP FORWARD IN TIME *************************************** DO 390 K=1,KBM1 DO 390 I=1,IM 390 VF(I,K)= 1 ((H(I)+ETB(I))*ARV(I)*VB(I,K) 2 -DTI2*VF(I,K)) 3 /((H(I)+ETF(I))*ARV(I)) RETURN END C################################################################# SUBROUTINE BAROPG(DRHOX) C REAL KM,KH,KQ,L PARAMETER (IM=41,KB=30,KS=20,IMM1=IM-1,KBM1=KB-1) PARAMETER (IMM2=IM-2,KBM2=KB-2) PARAMETER (IM2=IM*2,IMK=IM*KB,IMKM2=IM*KBM2) PARAMETER (IMKM1=IM*KBM1) COMMON/BLKCON/ 1 INT,IPRINT,DTE,DTI,TPRNU,UMOL,PR, 2 GRAV,TIME,RAMP C---------------- 1-D ARRAYS -------------------------------------- COMMON/BLK1D/ 1 DZR(KB),Z(KB),ZZ(KB),DZ(KB),DZZ(KB) C---------------- 2-D ARRAYS -------------------------------------- COMMON/BLK2D/H(IM),DX(IM),D(IM),DT(IM), 1 ART(IM),ARU(IM),ARV(IM),CBC(IM), 2 ALON(IM),ALAT(IM),ANG(IM), 3 DUM(IM),DVM(IM),FSM(IM),COR4(IM),CURV42D(IM), 4 WUSURF(IM),WVSURF(IM),WUBOT(IM),WVBOT(IM), 5 WTSURF(IM),WSSURF(IM),TPS(IM),AAM2D(IM), 6 UAF(IM),UA(IM),UAB(IM),VAF(IM),VA(IM), 7 VAB(IM),ELF(IM),EL(IM),ELB(IM), 8 ETF(IM),ET(IM),ETB(IM),FLUXUA(IM),FLUXVA(IM), 9 EGF(IM),EGB(IM), ELBB(IM) C---------------- 3-D ARRAYS -------------------------------------- COMMON/BLK3D/ 1 A(IM,KB),C(IM,KB),VH(IM,KB),VHP(IM,KB), 1 UF(IM,KB),VF(IM,KB), 2 KM(IM,KB),KH(IM,KB),KQ(IM,KB),L(IM,KB), 3 Q2(IM,KB),Q2B(IM,KB),AAM(IM,KB), 4 Q2L(IM,KB),Q2LB(IM,KB), 5 U(IM,KB),UB(IM,KB),W(IM,KB), 6 V(IM,KB),VB(IM,KB), 7 T(IM,KB),TB(IM,KB), 8 S(IM,KB),SB(IM,KB), 9 RHO(IM,KB),DTEF(IM,KB),RMEAN(IM,KB) C----------- 1 AND 2-D BOUNDARY VALUE ARRAYS ------------------------ COMMON/BDRY/ 1 TBE(KB),TBN(IM,KB),TBS(IM,KB),SBN(IM,KB),SBE(KB), 2 SBS(IM,KB),VABN(IM),VABS(IM),UGS(KB),VGS(IM,KB), 3 ELN(IM),ELS(IM),COVRHN(IM) DIMENSION DRHOX(IM,KB),D1(IM),D2(IM) C C X COMPONENT OF BAROCLINIC PRESSURE GRADIENT C DO I=1,IM DO K=1,KB RHO(I,K)=RHO(I,K)-RMEAN(I,K) DRHOX(K,1)=0. ENDDO ENDDO DO 300 I=2,IM 300 DRHOX(I,1)=.5*GRAV*(0.-ZZ(1))*(DT(I)+DT(I-1)) 2 *(RHO(I,1)-RHO(I-1,1)) 3 -.5*GRAV*ZZ(1)*(DT(I)-DT(I-1)) 4 *(RHO(I,2)+RHO(I-1,2)-RHO(I,1)-RHO(I-1,1)) 5 *(0.-ZZ(1))/(ZZ(2)-ZZ(1)) DO 315 K=2,KBM1 DO 310 I=2,IM D1(I)= +GRAV*.25*(ZZ(K-1)-ZZ(K))*(DT(I)+DT(I-1)) 2 *(RHO(I,K)-RHO(I-1,K)+RHO(I,K-1)-RHO(I-1,K-1)) D2(I)=+GRAV*.25*(ZZ(K-1)+ZZ(K))*(DT(I)-DT(I-1)) 4 *(RHO(I,K)+RHO(I-1,K)-RHO(I,K-1)-RHO(I-1,K-1)) 310 DRHOX(I,K)=DRHOX(I,K-1)+D1(I)+D2(I) c WRITE(6,'('' K,ZZ,RHO ='',I5,9E11.2)') c 1 K,ZZ(K),(RHO(I,K) ,I=10,13) c WRITE(6,'('' K,ZZ,D1,D2 ='',I5,9E11.2)') c 1 K,ZZ(K),(D1(I),D2(I),I=10,13) 315 CONTINUE DO 320 K=1,KBM1 DO 320 I=2,IM 320 DRHOX(I,K)=DRHOX(I,K)*.5*(DT(I)+DT(I-1))*DUM(I) c CALL PRXZ(' DRHOX ',0.,DRHOX,IM,1,KB,0.,DT,ZZ) DO I=1,IM DO K=1,KB RHO(I,K)=RHO(I,K)+RMEAN(I,K) ENDDO ENDDO C RETURN END C################################################################# SUBROUTINE BCOND(IDX) C REAL KM,KH,KQ,L PARAMETER (IM=41,KB=30,KS=20,IMM1=IM-1,KBM1=KB-1) PARAMETER (IMM2=IM-2,KBM2=KB-2) PARAMETER (IM2=IM*2,IMK=IM*KB,IMKM2=IM*KBM2) PARAMETER (IMKM1=IM*KBM1) COMMON/BLKCON/ 1 INT,IPRINT,DTE,DTI,TPRNU,UMOL,PR, 2 GRAV,TIME,RAMP C---------------- 1-D ARRAYS -------------------------------------- COMMON/BLK1D/ 1 DZR(KB),Z(KB),ZZ(KB),DZ(KB),DZZ(KB) C---------------- 2-D ARRAYS -------------------------------------- COMMON/BLK2D/H(IM),DX(IM),D(IM),DT(IM), 1 ART(IM),ARU(IM),ARV(IM),CBC(IM), 2 ALON(IM),ALAT(IM),ANG(IM), 3 DUM(IM),DVM(IM),FSM(IM),COR4(IM),CURV42D(IM), 4 WUSURF(IM),WVSURF(IM),WUBOT(IM),WVBOT(IM), 5 WTSURF(IM),WSSURF(IM),TPS(IM),AAM2D(IM), 6 UAF(IM),UA(IM),UAB(IM),VAF(IM),VA(IM), 7 VAB(IM),ELF(IM),EL(IM),ELB(IM), 8 ETF(IM),ET(IM),ETB(IM),FLUXUA(IM),FLUXVA(IM), 9 EGF(IM),EGB(IM), ELBB(IM) C---------------- 3-D ARRAYS -------------------------------------- COMMON/BLK3D/ 1 A(IM,KB),C(IM,KB),VH(IM,KB),VHP(IM,KB), 1 UF(IM,KB),VF(IM,KB), 2 KM(IM,KB),KH(IM,KB),KQ(IM,KB),L(IM,KB), 3 Q2(IM,KB),Q2B(IM,KB),AAM(IM,KB), 4 Q2L(IM,KB),Q2LB(IM,KB), 5 U(IM,KB),UB(IM,KB),W(IM,KB), 6 V(IM,KB),VB(IM,KB), 7 T(IM,KB),TB(IM,KB), 8 S(IM,KB),SB(IM,KB), 9 RHO(IM,KB),DTEF(IM,KB),RMEAN(IM,KB) C----------- 1 AND 2-D BOUNDARY VALUE ARRAYS ------------------------ COMMON/BDRY/ 1 TBE(KB),TBN(IM,KB),TBS(IM,KB),SBN(IM,KB),SBE(KB), 2 SBS(IM,KB),VABN(IM),VABS(IM),UGS(KB),VGS(IM,KB), 3 ELN(IM),ELS(IM),COVRHN(IM) C DATA PI/3.14167E0/,GEE/9.807E0/ 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 BC BC C C U(1) EL(1) U(2) C C V(1) C C U(1) EL(1) U(2) C C 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 -------------------------------------------------------------------- C U(1) IS NOT USED C C GO TO (10,20,30,40,50,60), IDX C 10 CONTINUE C----------------------------------------------------------------------- C EXTERNAL ELEV. B.C.'S C----------------------------------------------------------------------- C*SB ELF(IM)=ELF(IMM1) RETURN C 20 CONTINUE C----------------------------------------------------------------------- C EXTERNAL VEL B.C.'S C----------------------------------------------------------------------- C ***** SEAWARD ******** C*SB UAF(IM)=RAMP+EL(IMM1) C*SB UMID=UA(IM) C*SB VAF(IM)=VA(IM)-DTI/(DX(IM)+DX(IMM1)) C*SB 1 *( (UMID+ABS(UMID))*(VA(IM)-VA(IMM1)) C*SB 2 +(UMID-ABS(UMID))*(0.E0 -VA(IM)) ) C*120 VAF(IM)=0.E0 DO 131 I=1,IM UAF(I)=UAF(I)*DUM(I) 131 VAF(I)=VAF(I)*DVM(I) RETURN C 30 CONTINUE C----------------------------------------------------------------------- C INTERNAL VEL B.C.'S C----------------------------------------------------------------------- C **** SEAWARD ******* C*SB DO 142 K=1,KBM1 C*SB DO 140 =2,MM1 C*SB GA=SQRT(H(IM)/2000.E0) C*140 UF(IM,K) C*SB 1 =GA*(.25E0*U(IMM1,K)+.5E0*U(IMM1,K)+.25E0*U(IMM1,K)) C*SB 2 +(1.E0-GA)*(.25E0*U(IM,K)+.5E0*U(IM,K)+.25E0*U(IM,K)) C*SB UF(IM,M,K)=0.E0 C*SB UF(IM,1,K)=0.E0 C*SB DO 141 =1,M C*SB UMID=.5E0*(U(IM,K)+U(IM,K)) C*SB VF(IM,K)=V(IM,K)-DTI/(DX(IM)+DX(IMM1)) C*SB 1 *( (UMID+ABS(UMID))*(V(IM,K)-V(IMM1,K)) C*SB 2 +(UMID-ABS(UMID))*(0.E0 -V(IM,K)) ) C*141 CONTINUE C*142 CONTINUE C ********************** DO 160 K=1,KBM1 DO 160 I=1,IM UF(I,K)=UF(I,K)*DUM(I) VF(I,K)=VF(I,K)*DVM(I) 160 CONTINUE RETURN C 40 CONTINUE C----------------------------------------------------------------------- C TEMP & SAL B.C.'S C----------------------------------------------------------------------- C *** SEAWARD ***** C*SB DO 220 K=1,KBM1 C*SB DO 220 =1,M C*SB UF(IM,K)=T(IM,K)-DTI/(DX(IM)+DX(IMM1)) C*SB 1 *( (U(IM,K)+ABS(U(IM,K)))*(T(IM,K)-T(IMM1,K)) C*SB 2 +(U(IM,K)-ABS(U(IM,K)))*(TBE(,K) -T(IM,K)) ) C*SB VF(IM,K)=S(IM,K)-DTI/(DX(IM)+DX(IMM1)) C*SB 1 *( (U(IM,K)+ABS(U(IM,K)))*(S(IM,K)-S(IMM1,K)) C*SB 2 +(U(IM,K)-ABS(U(IM,K)))*(SBE(,K) -S(IM,K)) ) C DO 240 K=1,KBM1 DO 240 I=1,IM UF(I,K)=UF(I,K)*FSM(I) VF(I,K)=VF(I,K)*FSM(I) 240 CONTINUE RETURN C C 50 CONTINUE C---------------VERTICAL VEL. B. C.'S -------------------------------- DO 250 K=1,KBM1 CSB** DO 250 I=1,IM W(I,K)=W(I,K)*FSM(I) 250 CONTINUE CSB** DO 251 K=1,KB W(1,K)=0 251 W(IM,K)=0. CSB** RETURN C 60 CONTINUE C---------------- Q2 AND Q2L B.C.'S ----------------------------------- DO 297 K=1,KB DO 297 I=1,IM UF(1,K)=1.E-10 VF(1,K)=1.E-10 UF(IM,K)=1.E-10 VF(IM,K)=1.E-10 UF(I,K)=UF(I,K)*FSM(I) VF(I,K)=VF(I,K)*FSM(I) 297 CONTINUE RETURN END c SUBROUTINE DIAGN c VTOT=0.E0 c TTOT=0.E0 c STOT=0.E0 c MKE=0.E0 c MPE=0.E0 c DO 8888 K=1,KBM1 c DO 8888 I=1,IM c DVTOT=DX(I)*DT(I)*DZ(K) c VTOT=VTOT+DVTOT c MKE=MKE+(UB(I,K)**2+VB(I,K)**2)*DVTOT c MPE=MPE+GRAV*ZZ(K)*DT(I)*(RHO(I,K)-RMEAN(I,K))*DVTOT c8888 CONTINUE c TTOT=TTOT/VTOT c STOT=STOT/VTOT c MKE=MKE/VTOT c MPE=MPE/VTOT c WRITE(6,'('' VTOT ='',E20.8,'' TAVE ='',E20.8, c 1 '' SAVE ='',E20.8,'' MKEAVE ='',E20.8,/, c 2 '' MPEAVE ='',E20.8 )') VTOT,TTOT,STOT,MKE,MPE c RETURN c END SUBROUTINE DENS C REAL KM,KH,KQ,L PARAMETER (IM=41,KB=30,KS=20,IMM1=IM-1,KBM1=KB-1) PARAMETER (IMM2=IM-2,KBM2=KB-2) PARAMETER (IM2=IM*2,IMK=IM*KB,IMKM2=IM*KBM2) PARAMETER (IMKM1=IM*KBM1) COMMON/BLKCON/ 1 INT,IPRINT,DTE,DTI,TPRNU,UMOL,PR, 2 GRAV,TIME,RAMP C---------------- 1-D ARRAYS -------------------------------------- COMMON/BLK1D/ 1 DZR(KB),Z(KB),ZZ(KB),DZ(KB),DZZ(KB) C---------------- 2-D ARRAYS -------------------------------------- COMMON/BLK2D/H(IM),DX(IM),D(IM),DT(IM), 1 ART(IM),ARU(IM),ARV(IM),CBC(IM), 2 ALON(IM),ALAT(IM),ANG(IM), 3 DUM(IM),DVM(IM),FSM(IM),COR4(IM),CURV42D(IM), 4 WUSURF(IM),WVSURF(IM),WUBOT(IM),WVBOT(IM), 5 WTSURF(IM),WSSURF(IM),TPS(IM),AAM2D(IM), 6 UAF(IM),UA(IM),UAB(IM),VAF(IM),VA(IM), 7 VAB(IM),ELF(IM),EL(IM),ELB(IM), 8 ETF(IM),ET(IM),ETB(IM),FLUXUA(IM),FLUXVA(IM), 9 EGF(IM),EGB(IM), ELBB(IM) C---------------- 3-D ARRAYS -------------------------------------- COMMON/BLK3D/ 1 A(IM,KB),C(IM,KB),VH(IM,KB),VHP(IM,KB), 1 UF(IM,KB),VF(IM,KB), 2 KM(IM,KB),KH(IM,KB),KQ(IM,KB),L(IM,KB), 3 Q2(IM,KB),Q2B(IM,KB),AAM(IM,KB), 4 Q2L(IM,KB),Q2LB(IM,KB), 5 U(IM,KB),UB(IM,KB),W(IM,KB), 6 V(IM,KB),VB(IM,KB), 7 T(IM,KB),TB(IM,KB), 8 S(IM,KB),SB(IM,KB), 9 RHO(IM,KB),DTEF(IM,KB),RMEAN(IM,KB) C----------- 1 AND 2-D BOUNDARY VALUE ARRAYS ------------------------ COMMON/BDRY/ 1 TBE(KB),TBN(IM,KB),TBS(IM,KB),SBN(IM,KB),SBE(KB), 2 SBS(IM,KB),VABN(IM),VABS(IM),UGS(KB),VGS(IM,KB), 3 ELN(IM),ELS(IM),COVRHN(IM) REAL RHOF(IM,KB),TF(IM,KB),SF(IM,KB) EQUIVALENCE (A,RHOF),(VH,TF),(UF,SF) C C THIS SUBROUTINE COMPUTES DENSITY-1.025 C T = POTENTIAL TEMPERATURE C DO 1 K=1,KBM1 DO 1 I=1,IM TF(I,K)=T(I,K)+10.E0 SF(I,K)=S(I,K)+35.E0 C Approximate pressure in units of bars c P(I,K)=-GRAV*1.025*ZZ(K)*DT(I)*0.01 C RHOF(I,K) = 999.842594 + 6.793952E-2*TF(I,K) $ - 9.095290E-3*TF(I,K)**2 + 1.001685E-4*TF(I,K)**3 $ - 1.120083E-6*TF(I,K)**4 + 6.536332E-9*TF(I,K)**5 C RHOF(I,K) = RHOF(I,K) + (0.824493 - 4.0899E-3*TF(I,K) $ + 7.6438E-5*TF(I,K)**2 - 8.2467E-7*TF(I,K)**3 $ + 5.3875E-9*TF(I,K)**4) * SF(I,K) $ + (-5.72466E-3 + 1.0227E-4*TF(I,K) $ - 1.6546E-6*TF(I,K)**2) * SF(I,K)**1.5 $ + 4.8314E-4 * SF(I,K)**2 C c CR(I,J,K)=1449.1+.0821*P(I,J,K)+4.55*TF(I,J,K)-.045*TR(I,J,K)**2 c $ +1.34*(SF(I,J,K)-35.) c RHOF(I,J,K)=RHOF(I,J,K) + 1.E5*P(I,J,K)/CR(I,J,K)**2 c $ *(1.-2.*P(I,J,K)/CR(I,J,K)**2) C RHO(I,K)=(RHOF(I,K)-1025.)*1.E-3*FSM(I) 1 CONTINUE C DO 2 K=1,KBM1 DO 2 I=1,IM RHO(I,K)=RHOF(I,K)-25. RHO(I,K)=RHO(I,K)*1.E-3*FSM(I) 2 CONTINUE DO 3 I=1,IM 3 RHO(I,KB)=0. RETURN END *DECK DEPTH SUBROUTINE DEPTH(Z,ZZ,DZ,DZZ,DZR,KB,KBM1) CSB IMPIMCIT HALF PRECISION (A-H,O-Z) C >>> DIMENSION Z(KB),ZZ(KB),DZ(KB),DZZ(KB),DZR(KBM1) KL1=7 KL2=KB-7 C *** C THIS SUBROUTINE ESTABIMSHES THE VERTICAL RESOLUTION WITH LOG C DISTRIBUTIONS AT THE TOP AND BOTTOM AND A IMNEAR DISTRIBUTION C BETWEEN. CHOOSE KL1 AND KL2. THE DEFAULT KL1 = .3*KB AND KL2 = KB-2 C YIELDS A LOG DISTRIBUTION AT THE TOP AND NONE AT THE BOTTOM. C *** BB=FLOAT(KL2-KL1)+4. CC=FLOAT(KL1)-2. DEL1=2./BB/EXP(.693147*FLOAT(KL1-2)) DEL2=2./BB/EXP(.693147*FLOAT(KB-KL2-1)) Z(1)=0. ZZ(1)=-DEL1/2. DO 30 K=2,KL1 Z(K)=-DEL1*EXP(.693147*FLOAT(K-2)) 30 ZZ(K)=-DEL1*EXP(.693147*(FLOAT(K)-1.5)) DO 40 K=KL1,KL2 Z(K)=-(FLOAT(K)-CC)/BB 40 ZZ(K)=-(FLOAT(K)-CC+0.5)/BB DO 50 K=KL2,KBM1 Z(K)=(1.0-DEL2*EXP(.693147*FLOAT(KB-K-1)))*(-1.) 50 ZZ(K)=(1.0-DEL2*EXP(.693147*(FLOAT(KB-K)-1.5)))*(-1.) Z(KB)=-1.0 ZZ(KBM1)=-1.*(1.-DEL2/2.) ZZ(KB)=-1.*(1.+DEL2/2.) C *** c DO 55 K=1,KB c DENOM=FLOAT(KB-1) c Z(K)=-FLOAT(K-1)/DENOM c ZZ(K)=Z(K)-0.5/DENOM c 55 CONTINUE C *** DO 60 K=1,KBM1 DZ(K)=Z(K)-Z(K+1) DZR(K)=1./DZ(K) 60 DZZ(K)=ZZ(K)-ZZ(K+1) WRITE(6,70) 70 FORMAT(/2X,'K',8X,'Z',9X,'ZZ',9X,'DZ',9X,'DZZ',/) DO 90 K=1,KB WRITE(6,80) K,Z(K),ZZ(K),DZ(K),DZZ(K) 80 FORMAT(' ',I5,4F11.5) 90 CONTINUE RETURN END C############################################################ SUBROUTINE PROFQ(DT2) C REAL KM,KH,KQ,L PARAMETER (IM=41,KB=30,KS=20,IMM1=IM-1,KBM1=KB-1) PARAMETER (IMM2=IM-2,KBM2=KB-2) PARAMETER (IM2=IM*2,IMK=IM*KB,IMKM2=IM*KBM2) PARAMETER (IMKM1=IM*KBM1) COMMON/BLKCON/ 1 INT,IPRINT,DTE,DTI,TPRNU,UMOL,PR, 2 GRAV,TIME,RAMP C---------------- 1-D ARRAYS -------------------------------------- COMMON/BLK1D/ 1 DZR(KB),Z(KB),ZZ(KB),DZ(KB),DZZ(KB) C---------------- 2-D ARRAYS -------------------------------------- COMMON/BLK2D/H(IM),DX(IM),D(IM),DT(IM), 1 ART(IM),ARU(IM),ARV(IM),CBC(IM), 2 ALON(IM),ALAT(IM),ANG(IM), 3 DUM(IM),DVM(IM),FSM(IM),COR4(IM),CURV42D(IM), 4 WUSURF(IM),WVSURF(IM),WUBOT(IM),WVBOT(IM), 5 WTSURF(IM),WSSURF(IM),TPS(IM),AAM2D(IM), 6 UAF(IM),UA(IM),UAB(IM),VAF(IM),VA(IM), 7 VAB(IM),ELF(IM),EL(IM),ELB(IM), 8 ETF(IM),ET(IM),ETB(IM),FLUXUA(IM),FLUXVA(IM), 9 EGF(IM),EGB(IM), ELBB(IM) C---------------- 3-D ARRAYS -------------------------------------- COMMON/BLK3D/ 1 A(IM,KB),C(IM,KB),VH(IM,KB),VHP(IM,KB), 1 UF(IM,KB),VF(IM,KB), 2 KM(IM,KB),KH(IM,KB),KQ(IM,KB),L(IM,KB), 3 Q2(IM,KB),Q2B(IM,KB),AAM(IM,KB), 4 Q2L(IM,KB),Q2LB(IM,KB), 5 U(IM,KB),UB(IM,KB),W(IM,KB), 6 V(IM,KB),VB(IM,KB), 7 T(IM,KB),TB(IM,KB), 8 S(IM,KB),SB(IM,KB), 9 RHO(IM,KB),DTEF(IM,KB),RMEAN(IM,KB) C----------- 1 AND 2-D BOUNDARY VALUE ARRAYS ------------------------ COMMON/BDRY/ 1 TBE(KB),TBN(IM,KB),TBS(IM,KB),SBN(IM,KB),SBE(KB), 2 SBS(IM,KB),VABN(IM),VABS(IM),UGS(KB),VGS(IM,KB), 3 ELN(IM),ELS(IM),COVRHN(IM) REAL KAPPA,KN C DIMENSION GH(IM,KB),SM(IM,KB),SH(IM,KB) DIMENSION PROD(IM,KB),KN(IM,KB),BOYGR(IM,KB) DIMENSION DH(IM),CC(IM,KB) C EQUIVALENCE (A,SM),(C,SH),(PROD,KN),(TPS,DH) C EQUIVALENCE (DTEF,CC,GH) DATA PROD/IMK*0./,BOYGR/IMK*0./,GH/IMK*0./,KM/IMK*0./ DATA KH/IMK*0./,SM/IMK*0./,SH/IMK*0./ 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 C DO 50 I=1,IM 50 DH(I)=H(I)+ETF(I) DO 100 K=2,KBM1 DO 100 I=1,IM A(I,K)=-DT2*(KQ(I,K+1)+KQ(I,K)+2.*UMOL)*.5 1 /(DZZ(K-1)*DZ(K)*DH(I)*DH(I)) C(I,K)=-DT2*(KQ(I,K-1)+KQ(I,K)+2.*UMOL)*.5 1 /(DZZ(K-1)*DZ(K-1)*DH(I)*DH(I)) 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 I=1,IM VH(I,1)=0. VHP(I,1)=SQRT( (.5*(WUSURF(I)+WUSURF(I+1)))**2 1 +(.5*(WVSURF(I)+WVSURF(I)))**2 )*CONST1 UF(I,KB)=SQRT( (.5*(WUBOT(I)+WUBOT(I+1)))**2 1 +(.5*(WVBOT(I)+WVBOT(I)))**2 )*CONST1 90 CONTINUE C--------------------------------------------------------------------- C RESTORE 10 DEG, 35 PPT AND 1.025 BIAS IN T, S AND RHO C--------------------------------------------------------------------- c DO 101 K=1,KBM1 c DO 101 I=1,IM c TP=T(I,K)+10. C Calculate pressure in units of decibars c P=-GEE*1.025*ZZ(K)*DH(I)*.1 c CC(I,K)=1449.1+.00821*P+4.55*TP -.045*TP**2 c $ +1.34*(S(I,K)- 0.) c CC(I,K)=CC(I,K)/SQRT((1.-.01642*P/CC(I,K)) c $ *(1.-0.40*P/CC(I,K)**2)) c 101 CONTINUE DO 102 K=2,KBM1 DO 102 I=2,IM Q2B(I,K)=ABS(Q2B(I,K)) Q2LB(I,K)=ABS(Q2LB(I,K)) BOYGR(I,K)=GEE*((RHO(I,K-1)-RHO(I,K))/(DZZ(K-1)*DH(I))) C 1 +GEE**2*2.*1.025/(CC(I,K-1)**2+CC(I,K)**2) 102 CONTINUE C------ CALC. T.K.E. PRODUCTION ----------------------------------- DO 120 K=2,KBM1 DO 120 I=2,IM PROD(I,K)=KM(I,K)*.25*SEF 1 *( (U(I,K)-U(I,K-1)+U(I+1,K)-U(I+1,K-1))**2 2 +(V(I,K)-V(I,K-1)+V(I,K)-V(I,K-1))**2 ) 3 /(DZZ(K-1)*DH(I))**2 120 PROD(I,K)=PROD(I,K)+KH(I,K)*BOYGR(I,K) DO 110 K=2,KBM1 DO 110 I=1,IM 110 DTEF(I,K)=Q2B(I,K)*SQRT(Q2B(I,K))/(B1*Q2LB(I,K)+SMALL) DO 140 K=2,KBM1 DO 140 I=1,IM VHP(I,K)=1./(A(I,K)+C(I,K)*(1.-VH(I,K-1)) 1 -(2.*DT2*DTEF(I,K)+1.) ) VH(I,K)=A(I,K)*VHP(I,K) VHP(I,K)=(-2.*DT2*PROD(I,K) 1 +C(I,K)*VHP(I,K-1)-UF(I,K))*VHP(I,K) 140 CONTINUE DO 150 K=1,KBM1 KI=KB-K DO 150 I=1,IM UF(I,KI)=VH(I,KI)*UF(I,KI+1)+VHP(I,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 I=1,IM VH(I,1)=0.0 155 VHP(I,1)=0.0 DO 160 K=2,KBM1 DO 160 I=1,IM DTEF(I,K) =DTEF(I,K)*(1.+E2*((1./ABS(Z(K)-Z(1))+ 1 1./ABS(Z(K)-Z(KB))) *L(I,K)/(DH(I)*KAPPA))**2) VHP(I,K)=1./(A(I,K)+C(I,K)*(1.-VH(I,K-1)) 1 -(DT2*DTEF(I,K)+1.)) VH(I,K)=A(I,K)*VHP(I,K) VHP(I,K)=(DT2*(-PROD(I,K) 1 *L(I,K)*E1)+C(I,K)*VHP(I,K-1)-VF(I,K))*VHP(I,K) 160 CONTINUE DO 170 K=1,KBM1 KI=KB-K DO 170 I=1,IM VF(I,KI)=VH(I,KI)*VF(I,KI+1)+VHP(I,KI) 170 CONTINUE DO 180 K=2,KBM1 DO 180 I=1,IM IF(UF(I,K).LE.SMALL. OR.VF(I,K).LE.SMALL) THEN UF(I,K)=SMALL VF(I,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,KBM1 DO 220 I=2,IMM1 L(I,K)=VF(I,K)/UF(I,K) 220 GH(I,K)=L(I,K)**2/UF(I,K)*BOYGR(I,K) DO 230 K=2,KBM1 DO 230 I=2,IMM1 GH(I,K)=MIN(GH(I,K),.028) SH(I,K)=COEF1/(1.-COEF2*GH(I,K)) SM(I,K)=COEF3+SH(I,K)*COEF4*GH(I,K) SM(I,K)=SM(I,K)/(1.-COEF5*GH(I,K)) 230 CONTINUE C theta=0.5 seems to be the least noisy theta=0.50 DO 280 K=2,KBM1 DO 280 I=2,IMM1 c Using UF instead of Q2 seems be less noisy KN(I,K)=L(I,K)*SQRT(ABS(UF(I,K))) KQ(I,K)=(1.-theta)*KN(I,K)*.41*SM(I,K)+theta*KQ(I,K) KM(I,K)=(1.-theta)*KN(I,K)*SM(I,K)+theta*KM(I,K) KH(I,K)=(1.-theta)*KN(I,K)*SH(I,K)+theta*KH(I,K) 280 CONTINUE C IF(INT.EQ.2592.OR.INT.EQ.15552) THEN C CALL PRXZ(' Q2 ',TIME,Q2,IM,1,KB,0.,DT,Z) C CALL PRXZ(' L ',TIME,L ,IM,1,KB,0.,DT,Z) C CALL PRXZ(' BOYGR ',TIME,BOYGR ,IM,1,KB,0.,DT,Z) C CALL PRXZ(' GH ',TIME, GH ,IM,1,KB,10.,DT,Z) C CALL PRXZ(' SM ',TIME, SM ,IM,1,KB,0.,DT,Z) C CALL PRXZ(' KM ',TIME, KM ,IM,1,KB,0.,DT,Z) C ENDIF RETURN END *DECK PROFT SUBROUTINE PROFT(F,WFSURF,FSURF,NBC,DT2) C REAL KM,KH,KQ,L PARAMETER (IM=41,KB=30,KS=20,IMM1=IM-1,KBM1=KB-1) PARAMETER (IMM2=IM-2,KBM2=KB-2) PARAMETER (IM2=IM*2,IMK=IM*KB,IMKM2=IM*KBM2) PARAMETER (IMKM1=IM*KBM1) COMMON/BLKCON/ 1 INT,IPRINT,DTE,DTI,TPRNU,UMOL,PR, 2 GRAV,TIME,RAMP C---------------- 1-D ARRAYS -------------------------------------- COMMON/BLK1D/ 1 DZR(KB),Z(KB),ZZ(KB),DZ(KB),DZZ(KB) C---------------- 2-D ARRAYS -------------------------------------- COMMON/BLK2D/H(IM),DX(IM),D(IM),DT(IM), 1 ART(IM),ARU(IM),ARV(IM),CBC(IM), 2 ALON(IM),ALAT(IM),ANG(IM), 3 DUM(IM),DVM(IM),FSM(IM),COR4(IM),CURV42D(IM), 4 WUSURF(IM),WVSURF(IM),WUBOT(IM),WVBOT(IM), 5 WTSURF(IM),WSSURF(IM),TPS(IM),AAM2D(IM), 6 UAF(IM),UA(IM),UAB(IM),VAF(IM),VA(IM), 7 VAB(IM),ELF(IM),EL(IM),ELB(IM), 8 ETF(IM),ET(IM),ETB(IM),FLUXUA(IM),FLUXVA(IM), 9 EGF(IM),EGB(IM), ELBB(IM) C---------------- 3-D ARRAYS -------------------------------------- COMMON/BLK3D/ 1 A(IM,KB),C(IM,KB),VH(IM,KB),VHP(IM,KB), 1 UF(IM,KB),VF(IM,KB), 2 KM(IM,KB),KH(IM,KB),KQ(IM,KB),L(IM,KB), 3 Q2(IM,KB),Q2B(IM,KB),AAM(IM,KB), 4 Q2L(IM,KB),Q2LB(IM,KB), 5 U(IM,KB),UB(IM,KB),W(IM,KB), 6 V(IM,KB),VB(IM,KB), 7 T(IM,KB),TB(IM,KB), 8 S(IM,KB),SB(IM,KB), 9 RHO(IM,KB),DTEF(IM,KB),RMEAN(IM,KB) C----------- 1 AND 2-D BOUNDARY VALUE ARRAYS ------------------------ COMMON/BDRY/ 1 TBE(KB),TBN(IM,KB),TBS(IM,KB),SBN(IM,KB),SBE(KB), 2 SBS(IM,KB),VABN(IM),VABS(IM),UGS(KB),VGS(IM,KB), 3 ELN(IM),ELS(IM),COVRHN(IM) C REAL AF(IM,KB),CF(IM,KB),VHF(IM,KB),VHPF(IM,KB) DIMENSION F(IM,KB),WFSURF(IM),FSURF(IM),DH(IM) EQUIVALENCE (A,AF),(RHO,CF),(VH,VHF) EQUIVALENCE (TPS,DH) UMOLPR=UMOL/PR C********************************************************************** C * C THE FOLLOWING SECTION SOLVES THE EQUATION * C DT2*(KH*F')'-F=-FB * C * C********************************************************************** CSB DO 90 I=1,IM DH(I)=H(I)+ETF(I) 90 CONTINUE DO 100 K=2,KBM1 CSB DO 100 I=1,IM AF(I,K-1)=-DT2*(KH(I,K)+UMOLPR)/(DZ(K-1)*DZZ(K-1)*DH(I) 1 *DH(I)) CF(I,K)=-DT2*(KH(I,K)+UMOLPR)/(DZ(K)*DZZ(K-1)*DH(I) 1 *DH(I)) 100 CONTINUE C----- SURFACE BC'S; WFSURF OR FSURF ----------------------------------- GO TO (50,51), NBC 50 CONTINUE CSB DO 501 I=1,IM VHF(I,1)=AF(I,1)/(AF(I,1)-1.) VHPF(I,1)=-DT2*WFSURF(I)/(-DZ(1)*DH(I))-F(I,1) 501 VHPF(I,1)=VHPF(I,1)/(AF(I,1)-1.) GO TO 52 51 CONTINUE CSB DO 511 I=1,IM VHF(I,1)=0. 511 VHPF(I,1)=FSURF(I) 52 CONTINUE C---------------------------------------------------------------------- DO 101 K=2,KBM2 CSB DO 101 I=1,IM VHPF(I,K)=1./(AF(I,K)+CF(I,K)*(1.-VHF(I,K-1))-1.) VHF(I,K)=AF(I,K)*VHPF(I,K) VHPF(I,K)=(CF(I,K)*VHPF(I,K-1)-REAL(F(I,K)))*VHPF(I,K) 101 CONTINUE CSB DO 102 I=1,IM 102 F(I,KBM1)=((CF(I,KBM1)*VHPF(I,KBM2)-F(I,KBM1)) 1 /(CF(I,KBM1)*(1.-VHF(I,KBM2))-1.)) DO 105 K=2,KBM1 KI=KB-K DO 105 I=1,IM F(I,KI)=(VHF(I,KI)*F(I,KI+1)+VHPF(I,KI)) 105 CONTINUE CSB RETURN END C############################################################### SUBROUTINE PROFU(DT2) C REAL KM,KH,KQ,L PARAMETER (IM=41,KB=30,KS=20,IMM1=IM-1,KBM1=KB-1) PARAMETER (IMM2=IM-2,KBM2=KB-2) PARAMETER (IM2=IM*2,IMK=IM*KB,IMKM2=IM*KBM2) PARAMETER (IMKM1=IM*KBM1) COMMON/BLKCON/ 1 INT,IPRINT,DTE,DTI,TPRNU,UMOL,PR, 2 GRAV,TIME,RAMP C---------------- 1-D ARRAYS -------------------------------------- COMMON/BLK1D/ 1 DZR(KB),Z(KB),ZZ(KB),DZ(KB),DZZ(KB) C---------------- 2-D ARRAYS -------------------------------------- COMMON/BLK2D/H(IM),DX(IM),D(IM),DT(IM), 1 ART(IM),ARU(IM),ARV(IM),CBC(IM), 2 ALON(IM),ALAT(IM),ANG(IM), 3 DUM(IM),DVM(IM),FSM(IM),COR4(IM),CURV42D(IM), 4 WUSURF(IM),WVSURF(IM),WUBOT(IM),WVBOT(IM), 5 WTSURF(IM),WSSURF(IM),TPS(IM),AAM2D(IM), 6 UAF(IM),UA(IM),UAB(IM),VAF(IM),VA(IM), 7 VAB(IM),ELF(IM),EL(IM),ELB(IM), 8 ETF(IM),ET(IM),ETB(IM),FLUXUA(IM),FLUXVA(IM), 9 EGF(IM),EGB(IM), ELBB(IM) C---------------- 3-D ARRAYS -------------------------------------- COMMON/BLK3D/ 1 A(IM,KB),C(IM,KB),VH(IM,KB),VHP(IM,KB), 1 UF(IM,KB),VF(IM,KB), 2 KM(IM,KB),KH(IM,KB),KQ(IM,KB),L(IM,KB), 3 Q2(IM,KB),Q2B(IM,KB),AAM(IM,KB), 4 Q2L(IM,KB),Q2LB(IM,KB), 5 U(IM,KB),UB(IM,KB),W(IM,KB), 6 V(IM,KB),VB(IM,KB), 7 T(IM,KB),TB(IM,KB), 8 S(IM,KB),SB(IM,KB), 9 RHO(IM,KB),DTEF(IM,KB),RMEAN(IM,KB) C----------- 1 AND 2-D BOUNDARY VALUE ARRAYS ------------------------ COMMON/BDRY/ 1 TBE(KB),TBN(IM,KB),TBS(IM,KB),SBN(IM,KB),SBE(KB), 2 SBS(IM,KB),VABN(IM),VABS(IM),UGS(KB),VGS(IM,KB), 3 ELN(IM),ELS(IM),COVRHN(IM) DIMENSION DH(IM) DATA UMOLB/1.36E-6 / C*********************************************************************** C * C THE FOLLOWING SECTION SOLVES THE EQUATION * C DT2*(KM*U')'-U=-UB * C * C*********************************************************************** CSB DO 84 I=1,IM 84 DH(I)=1.0 DO 85 I=2,IM 85 DH(I)=.5E0*(H(I)+ETF(I)+H(I-1)+ETF(I-1)) DO 90 K=1,KB DO 90 I=2,IM 90 C(I,K)=(KM(I,K)+KM(I-1,K))*.5E0 DO 100 K=2,KBM1 CSB DO 100 I=1,IM A(I,K-1)=-DT2*(C(I,K)+UMOL )/(DZ(K-1)*DZZ(K-1)*DH(I) 1 *DH(I)) C(I,K)=-DT2*(C(I,K)+UMOL )/(DZ(K)*DZZ(K-1)*DH(I) 1 *DH(I)) 100 CONTINUE CSB DO 1011 I=1,IM VH(I,1)=A(I,1)/(A(I,1)-1.E0) VHP(I,1)=(-DT2*WUSURF(I)/(-DZ(1)*DH(I))-UF(I,1)) 1 /(A(I,1)-1.E0) 1011 CONTINUE DO 101 K=2,KBM2 DO 101 I=2,IM VHP(I,K)=1.E0/(A(I,K)+C(I,K)*(1.-VH(I,K-1))-1.) VH(I,K)=A(I,K)*VHP(I,K) VHP(I,K)=(C(I,K)*VHP(I,K-1)-UF(I,K))*VHP(I,K) 101 CONTINUE DO 102 I=2,IM TPS(I)=0.5E0*(CBC(I)+CBC(I-1)) 1 *SQRT(UB(I,KBM1)**2+(.25E0*(VB(I,KBM1) 2 +VB(I,KBM1)+VB(I-1,KBM1)+VB(I-1,KBM1)))**2) UF(I,KBM1)=UF(I,KBM2) UF(I,KBM1)=(C(I,KBM1)*VHP(I,KBM2)-UF(I,KBM1))/(TPS(I) 1 *DT2/(-DZ(KBM1)*DH(I))-1.E0-(VH(I,KBM2)-1.E0)*C(I,KBM1)) 102 UF(I,KBM1)=UF(I,KBM1)*DUM(I) DO 103 K=2,KBM1 KI=KB-K DO 103 I=1,IM UF(I,KI)=(VH(I,KI)*UF(I,KI+1)+VHP(I,KI))*DUM(I) 103 CONTINUE CSB DO 104 I=1,IM 104 WUBOT(I)=-TPS(I)*UF(I,KBM1) c104 WUBOT(I)=0.0 RETURN END *DECK PROFV SUBROUTINE PROFV(DT2) C REAL KM,KH,KQ,L PARAMETER (IM=41,KB=30,KS=20,IMM1=IM-1,KBM1=KB-1) PARAMETER (IMM2=IM-2,KBM2=KB-2) PARAMETER (IM2=IM*2,IMK=IM*KB,IMKM2=IM*KBM2) PARAMETER (IMKM1=IM*KBM1) COMMON/BLKCON/ 1 INT,IPRINT,DTE,DTI,TPRNU,UMOL,PR, 2 GRAV,TIME,RAMP C---------------- 1-D ARRAYS -------------------------------------- COMMON/BLK1D/ 1 DZR(KB),Z(KB),ZZ(KB),DZ(KB),DZZ(KB) C---------------- 2-D ARRAYS -------------------------------------- COMMON/BLK2D/H(IM),DX(IM),D(IM),DT(IM), 1 ART(IM),ARU(IM),ARV(IM),CBC(IM), 2 ALON(IM),ALAT(IM),ANG(IM), 3 DUM(IM),DVM(IM),FSM(IM),COR4(IM),CURV42D(IM), 4 WUSURF(IM),WVSURF(IM),WUBOT(IM),WVBOT(IM), 5 WTSURF(IM),WSSURF(IM),TPS(IM),AAM2D(IM), 6 UAF(IM),UA(IM),UAB(IM),VAF(IM),VA(IM), 7 VAB(IM),ELF(IM),EL(IM),ELB(IM), 8 ETF(IM),ET(IM),ETB(IM),FLUXUA(IM),FLUXVA(IM), 9 EGF(IM),EGB(IM), ELBB(IM) C---------------- 3-D ARRAYS -------------------------------------- COMMON/BLK3D/ 1 A(IM,KB),C(IM,KB),VH(IM,KB),VHP(IM,KB), 1 UF(IM,KB),VF(IM,KB), 2 KM(IM,KB),KH(IM,KB),KQ(IM,KB),L(IM,KB), 3 Q2(IM,KB),Q2B(IM,KB),AAM(IM,KB), 4 Q2L(IM,KB),Q2LB(IM,KB), 5 U(IM,KB),UB(IM,KB),W(IM,KB), 6 V(IM,KB),VB(IM,KB), 7 T(IM,KB),TB(IM,KB), 8 S(IM,KB),SB(IM,KB), 9 RHO(IM,KB),DTEF(IM,KB),RMEAN(IM,KB) C----------- 1 AND 2-D BOUNDARY VALUE ARRAYS ------------------------ COMMON/BDRY/ 1 TBE(KB),TBN(IM,KB),TBS(IM,KB),SBN(IM,KB),SBE(KB), 2 SBS(IM,KB),VABN(IM),VABS(IM),UGS(KB),VGS(IM,KB), 3 ELN(IM),ELS(IM),COVRHN(IM) DIMENSION DH(IM) DATA UMOLB/1.36E-6/ C*********************************************************************** C * C THE FOLLOWING SECTION SOLVES THE EQUATION * C DT2*(KM*U')'-U=-UB * C * C*********************************************************************** CSB DO 84 I=1,IM 84 DH(I)=1. DO 85 I=1,IM 85 DH(I)=H(I)+ETF(I) DO 90 K=1,KB DO 90 I=1,IM 90 C(I,K)=KM(I,K) DO 100 K=2,KBM1 CSB DO 100 I=1,IM A(I,K-1)=-DT2*(C(I,K)+UMOL )/(DZ(K-1)*DZZ(K-1)*DH(I) 1 *DH(I)) C(I,K)=-DT2*(C(I,K)+UMOL )/(DZ(K)*DZZ(K-1)*DH(I) 1 *DH(I)) 100 CONTINUE CSB DO 1001 I=1,IM VH(I,1)=A(I,1)/(A(I,1)-1.E0) VHP(I,1)=(-DT2*WVSURF(I)/(-DZ(1)*DH(I))-VF(I,1)) 1 /(A(I,1)-1.E0) 1001 CONTINUE DO 101 K=2,KBM2 CSB DO 101 I=1,IM VHP(I,K)=1.E0/(A(I,K)+C(I,K)*(1.E0-VH(I,K-1))-1.E0) VH(I,K)=A(I,K)*VHP(I,K) VHP(I,K)=(C(I,K)*VHP(I,K-1)-VF(I,K))*VHP(I,K) 101 CONTINUE DO 102 I=1,IMM1 TPS(I)=CBC(I) 1 *SQRT((.25E0*(UB(I,KBM1)+UB(I+1,KBM1) 2 +UB(I,KBM1)+UB(I+1,KBM1)))**2+VB(I,KBM1)**2) VF(I,KBM1)=VF(I,KBM2) VF(I,KBM1)=(C(I,KBM1)*VHP(I,KBM2)-VF(I,KBM1))/(TPS(I) 1 *DT2/(-DZ(KBM1)*DH(I))-1.E0-(VH(I,KBM2)-1.E0)*C(I,KBM1)) 102 VF(I,KBM1)=VF(I,KBM1)*DVM(I) DO 103 K=2,KBM1 KI=KB-K DO 103 I=1,IM VF(I,KI)=(VH(I,KI)*VF(I,KI+1)+VHP(I,KI))*DVM(I) 103 CONTINUE CSB DO 104 I=1,IM 104 WVBOT(I)=-TPS(I)*VF(I,KBM1) c104 WVBOT(I)=0.0 RETURN END C############################################################# SUBROUTINE PRXZ(LABEL,TIME,A,IM,ISKP,KB,SCALA,DT,ZZ,LAB) DIMENSION A(IM,KB),NUM(100),LINE(100),IDT(100), 1 ZZ(KB),DT(IM) CHARACTER LABEL*(*) C C THIS SUBROUTINE WRITES VERTICAL SECTIONS OF A 3-D FIELD C C TIME=TIME IN DAYS C A= ARRAY(IM,KB) TO BE PRINTED C ISPL=PRINT SKIP FOR I C SCALE=DIVISOR FOR VALUES OF A C LAB=1 include labels C DATA ZERO /1.E-10/ C SCALE=SCALA IF(SCALE.GT.ZERO) GO TO 160 AMX=ZERO DO 150 K=1,KB DO 150 I=1,IM,ISKP AMX=MAX(ABS(A(I,K)),AMX) 150 CONTINUE SCALE=10.**(INT(LOG10(AMX)+1.E2)-103) 160 CONTINUE SCALEI=1./SCALE IF(LAB.EQ.1) THEN WRITE(6,9980) LABEL 9980 FORMAT(' ',A30) WRITE(6,9981) TIME,SCALE 9981 FORMAT(' TIME = ',F9.3,' DAYS MULTIPLY ALL VALUES BY',1PE10.3) ENDIF IB=1 IE=IB+23*ISKP 50 CONTINUE IF(IE.GT.IM) IE=IM DO 100 I=IB,IE,ISKP IDT(I)=DT(I) 100 NUM(I)=I IF(LAB.EQ.1) THEN 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,/,' ZZ ') ENDIF DO 200 K=1,KB DO 190 I=IB,IE,ISKP 190 LINE(I)=INT(SCALEI*A(I,K)) WRITE(6,9966) K,ZZ(K),(LINE(I),I=IB,IE,ISKP) 9966 FORMAT(1X,I2,2X,F6.3,24I5) 200 CONTINUE IF(LAB.EQ.1) WRITE(6,'(//)') IF(IE.EQ.IM) GO TO 10 IB=IB+24*ISKP IE=IB+23*ISKP GO TO 50 10 CONTINUE RETURN END C############################################################# SUBROUTINE PRXZSH(LABEL,TIME,A,IM,ISKP,KB,SCALA,DT,ZZ,LAB) DIMENSION A(IM,KB),NUM(100),LINE(100),IDT(100), 1 ZZ(KB),DT(IM) CHARACTER LABEL*(*) C C THIS SUBROUTINE WRITES VERTICAL SECTIONS OF A 3-D FIELD C C TIME=TIME IN DAYS C A= ARRAY(IM,KB) TO BE PRINTED C ISPL=PRINT SKIP FOR I C SCALE=DIVISOR FOR VALUES OF A C DATA ZERO /1.E-10/ C IL=21 IF(KB.GT.1)THEN KF=KB-15 ELSE KF=1 ENDIF C SCALE=SCALA IF(SCALE.GT.ZERO) GO TO 160 AMX=ZERO DO 150 K=KF,KB DO 150 I=1,IL,ISKP AMX=MAX(ABS(A(I,K)),AMX) 150 CONTINUE SCALE=10.**(INT(LOG10(AMX)+1.E2)-103) 160 CONTINUE SCALEI=1./SCALE IF(LAB.EQ.1) THEN WRITE(6,9980) LABEL 9980 FORMAT(' ',A30) WRITE(6,9981) TIME,SCALE 9981 FORMAT(' TIME = ',F9.3,' DAYS MULTIPLY ALL VALUES BY',1PE10.3) ENDIF IB=1 IE=IB+23*ISKP 50 CONTINUE IF(IE.GT.IL) IE=IL DO 100 I=IB,IE,ISKP IDT(I)=DT(I) 100 NUM(I)=I IF(LAB.EQ.1) THEN 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,/,' ZZ ') ENDIF DO 200 K=KF,KB DO 190 I=IB,IE,ISKP 190 LINE(I)=INT(SCALEI*A(I,K)) WRITE(6,9966) K,ZZ(K),(LINE(I),I=IB,IE,ISKP) 9966 FORMAT(1X,I2,2X,F8.5,24I5) 200 CONTINUE IF(LAB.EQ.1) WRITE(6,'(//)') IF(IE.EQ.IL) GO TO 10 IB=IB+24*ISKP IE=IB+23*ISKP GO TO 50 10 CONTINUE RETURN END C######################################################### SUBROUTINE SHPFIL(UR,UI,IM,KB,NORD) C C THIS IS A SHAPIRO FILTER WHICH REMOVES SMALL SCALE C NOISE (TWO-DELTA-X NOISE IS COMPLETELY ELIMINATED) C WHILE REDUCING FILTER EFFECT ON LARGER SCALES. C PARAMETER(II=41) DIMENSION UI(IM,KB),UR(IM,KB) COMPLEX FLUX(II),VCM(II),FLCON(8) DATA FLCON/ (.5,0.), (-.5,0.), (0.,.5), (0.,-.5), 1 (.354,.354), (-.354,-.354), (-.354,.354), (.354,-.354)/ C DO 1000 K=1,KB C DO 70 I=1,IM 70 VCM(I)=CMPLX(UR(I,K),UI(I,K)) C DO 500 LOOP=1,NORD C C TIME=FLOAT(LOOP)+.0001 C WRITE(6,'(1X,I5,2F10.3)') LOOP,FLCON(LOOP) C DO 100 I=1,IM-1 100 FLUX(I)= FLCON(LOOP)*(VCM(I+1)-VCM(I)) DO 120 I=2,IM-1 120 VCM(I)=VCM(I) 1 +0.5*(FLUX(I)-FLUX(I-1)) 500 CONTINUE C DO 10 I=1,IM UR(I,K)=REAL(VCM(I)) 10 UI(I,K)=AIMAG(VCM(I)) 1000 CONTINUE RETURN END C *DECK VERTVL SUBROUTINE VERTVL(DTI2) C REAL KM,KH,KQ,L PARAMETER (IM=41,KB=30,KS=20,IMM1=IM-1,KBM1=KB-1) PARAMETER (IMM2=IM-2,KBM2=KB-2) PARAMETER (IM2=IM*2,IMK=IM*KB,IMKM2=IM*KBM2) PARAMETER (IMKM1=IM*KBM1) COMMON/BLKCON/ 1 INT,IPRINT,DTE,DTI,TPRNU,UMOL,PR, 2 GRAV,TIME,RAMP C---------------- 1-D ARRAYS -------------------------------------- COMMON/BLK1D/ 1 DZR(KB),Z(KB),ZZ(KB),DZ(KB),DZZ(KB) C---------------- 2-D ARRAYS -------------------------------------- COMMON/BLK2D/H(IM),DX(IM),D(IM),DT(IM), 1 ART(IM),ARU(IM),ARV(IM),CBC(IM), 2 ALON(IM),ALAT(IM),ANG(IM), 3 DUM(IM),DVM(IM),FSM(IM),COR4(IM),CURV42D(IM), 4 WUSURF(IM),WVSURF(IM),WUBOT(IM),WVBOT(IM), 5 WTSURF(IM),WSSURF(IM),TPS(IM),AAM2D(IM), 6 UAF(IM),UA(IM),UAB(IM),VAF(IM),VA(IM), 7 VAB(IM),ELF(IM),EL(IM),ELB(IM), 8 ETF(IM),ET(IM),ETB(IM),FLUXUA(IM),FLUXVA(IM), 9 EGF(IM),EGB(IM), ELBB(IM) C---------------- 3-D ARRAYS -------------------------------------- COMMON/BLK3D/ 1 A(IM,KB),C(IM,KB),VH(IM,KB),VHP(IM,KB), 1 UF(IM,KB),VF(IM,KB), 2 KM(IM,KB),KH(IM,KB),KQ(IM,KB),L(IM,KB), 3 Q2(IM,KB),Q2B(IM,KB),AAM(IM,KB), 4 Q2L(IM,KB),Q2LB(IM,KB), 5 U(IM,KB),UB(IM,KB),W(IM,KB), 6 V(IM,KB),VB(IM,KB), 7 T(IM,KB),TB(IM,KB), 8 S(IM,KB),SB(IM,KB), 9 RHO(IM,KB),DTEF(IM,KB),RMEAN(IM,KB) C----------- 1 AND 2-D BOUNDARY VALUE ARRAYS ------------------------ COMMON/BDRY/ 1 TBE(KB),TBN(IM,KB),TBS(IM,KB),SBN(IM,KB),SBE(KB), 2 SBS(IM,KB),VABN(IM),VABS(IM),UGS(KB),VGS(IM,KB), 3 ELN(IM),ELS(IM),COVRHN(IM) DIMENSION XFLUX(IM,KB),YFLUX(IM,KB) EQUIVALENCE (XFLUX,A),(YFLUX,C) C C C CALCULATE NEW VERTICAL VELOCITY C C REESTABIMSH BOUNDARY CONDITIONS DO 100 K=1,KBM1 DO 100 I=2,IM 100 XFLUX(I,K) 1 =.5E0*(DT(I)+DT(I-1))*U(I,K) DO 120 K=1,KBM1 DO 120 I=1,IM 120 YFLUX(I,K) 1 =DX(I)*DT(I)*V(I,K) CSB DO 125 I=1,IM 125 W(I,1)=0.E0 DO 710 K=1,KBM1 DO 710 I=1,IMM1 710 W(I,K+1)=W(I,K) 1 +DZ(K)*((XFLUX(I+1,K)-XFLUX(I,K)) 2 /DX(I) 3 +(ETF(I)-ETB(I))/DTI2 ) RETURN END C############################################################# SUBROUTINE ZTOSIG(ZS,TB,ZZ,H,T,IM,KS,KB) PARAMETER (KBB=30,KSS=20) DIMENSION ZS(KS),TB(IM,KS),ZZ(KB),H(IM),T(IM,KB), 1 TIN(KSS),TOUT(KBB),ZZH(KBB) DO 40 I=1,IM IF(H(I).LT.1.00001)GO TO 40 DO 45 K=1,KS TIN(K)=TB(I,K) IF(TIN(K).LT.0.00001.AND.K.NE.1)TIN(K)=TIN(K-1) IF(ZS(K).LT.H(I)) KBOT=K C IF(ZND(K).LT.-H(I).AND.K.NE.1)TIN(K)=TIN(K-1) 45 CONTINUE C DO 50 K=1,KB 50 ZZH(K)=-ZZ(K)*H(I) C CALL SPLINC(ZS,TIN,KS ,2.e30,2.e30,ZZH,TOUT,KB) C IPR=IM/2 IF(I.EQ.IPR) THEN WRITE(6,'(//,'' Data interpolated from z to sigma grid at I ='' 1 ,I5)') IPR WRITE(6,'('' H ='',F10.1)') H(IPR) WRITE(6,'(1X,I5,4F10.4)') (K,ZS(K),TIN(K),ZZH(K),TOUT(K),K=1,KS) WRITE(6,'(1X,I5,20X,2F10.4)') (K, ZZH(K),TOUT(K),K=KS+1,KB) ENDIF C DO K=1,KB T(I,K)=TOUT(K) T(I,KB)=0. ENDDO 40 CONTINUE C RETURN END C C SUBROUTINE SPLINC(X,Y,N,YP1,YPN,XNEW,YNEW,M) PARAMETER (NMAX=100) DIMENSION X(N),Y(N),Y2(NMAX),U(NMAX),XNEW(M),YNEW(M) 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 11 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 11 CONTINUE 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 12 K=N-1,1,-1 Y2(K)=Y2(K)*Y2(K+1)+U(K) 12 CONTINUE C DO 20 I =1,M CALL SPLINT(X,Y,Y2,N,XNEW(I),YNEW(I)) 20 CONTINUE RETURN END SUBROUTINE SPLINT(XA,YA,Y2A,N,X,Y) DIMENSION XA(N),YA(N),Y2A(N) 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) 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 SUBROUTINE SLPMIN(H,IM,FSM) DIMENSION H(IM),FSM(IM) DIMENSION SL(100) DO 3 LOOP=1,10 SLMIN=0.2 DO 1 I=2,IM-1 IF(FSM(I).EQ.0..OR.FSM(I+1).EQ.0.) GOTO 1 SL(I)=ABS(H(I+1)-H(I))/(H(I)+H(I+1)) IF(SL(I).LT.SLMIN) GOTO 1 DH=0.5*(SL(I)-SLMIN)*(H(I)+H(I+1)) SN=-1. IF(H(I+1).GT.H(I)) SN=1. H(I+1)=H(I+1)-SN*DH H(I)=H(I)+SN*DH 1 CONTINUE DO 2 I=IM-1,2,-1 IF(FSM(I).EQ.0..OR.FSM(I+1).EQ.0.) GOTO 2 SL(I)=ABS(H(I+1)-H(I))/(H(I)+H(I+1)) IF(SL(I).LT.SLMIN) GOTO 2 DH=0.5*(SL(I)-SLMIN)*(H(I)+H(I+1)) SN=-1. IF(H(I+1).GT.H(I)) SN=1. H(I+1)=H(I+1)-SN*DH H(I)=H(I)+SN*DH 2 CONTINUE 3 CONTINUE CALL PRXZ(' SL ',0.,SL ,IM,1,1 ,0., H,ZZ) RETURN END