*DECK MAIN PROGRAM MAIN C********************************************************************** C C THIS IS A 1-D VERSION OF THE MELLOR YAMADA, LEVEL 2.5 MODEL C AS APPLIED TO OCEAN BOUNDAY LAYERS. THE MOST RECENT DISCRIPTION C OF THE MODEL WILL BE FOUND IN REV. GEOPHYS. AND SPACE PHYS., C VOL.20,PP.851-875,1982.(THE MODEL CAN BE CONVERTED TO AN ATMOS- C PHERIC B.L. MODEL WITH DUE REGARD FOR THERMODYNAMICS AND THE DIRECT- C ION OF THE GRAVITY VECTOR.) C ESSENTIALLY 'THE' MODEL RESIDES IN SUBROUTINE PROFQ WHICH PROVID- C ES MIXING COEFFICIENTS, KM FOR VELOCITY AND KH FOR TEMPERATURE AND C SALINTY. SUBROUTINES PROFU AND PROFV SOLVE FOR THE U AND V VELOCITY C WHEREAS PROFT SOLVES FOR TEMPERATURE AND SALINITY (OR ANY OTHER C SCALER). ALL OF THESE "PROF" ROUTINES USE THE INTEGRATION METHOD C DESCRIBED ON P.198-201 0F RICHTMEYER AND MORTON, INTERSCIENCE PUB.,1967. C THE MODEL IS HERE CONFIGURED TO MATCH THE VELOCITY, TEMPERATURE C AND SALINITY OF UNDERLYING WATER. HOWEVER, SEE PROFU,PROFV AND C PROFQ FOR SIMPLE MODIFICATIONS TO INCLUDE BOTTOM BOUNDARY LAYERS C WHICH, WITH SUFFICIENT VERTICAL RESOLUTION, SHOULD PRODUCE A LOG C DISTRIBUTION NEAR THE BOTTOM. C THE NUMERICS ARE ADAPTED FROM THE 3-D MODEL OF BLUMBERG AND C MELLOR WHICH IS IMPLICIT IN THE VERTICAL DIFFUSION TERMS AND LEAP C FROG IN THE ADVECTIVE AND HORIZONTAL DIFFUSION TERMS. SINCE THE C LATTER ARE EXCLUDED HERE, THE CODE MAY BE CONVERTED TO A TWO LEVEL C TIME STEP ALGORITHM FAIRLY EASILY. C THERE ARE NUMERICAL PARAMETERS SUCH AS CHOICE OF GRID C AND PHYSICAL PARAMETERS SUCH AS THE GRAVITY CONSTANT TO BE CHOSEN C BY THE USER. THE MODEL PARAMETERS IN PROFQ MAY BE FIDDLED WITH, C OF COURSE, BUT WE PREFER TO THINK OF THEM AS FIXED. THE PARAMETER, C UMOL, IS "BACKROUND VISCOSITY" PRESUMABLY RELATED TO INTERNAL WAVES C OR WHATEVER. EVIDENCE IN THE OCEAN PLACES IT THE RANGE, 1.E-5 TO C 1.E-4 MKS; IT IS, HOWEVER, AN IGNORANCE FACTOR TO BE CHOSEN BY THE C USER. C G. MELLOR, MARCH 1983 C C This code is not supported in the sense of pom.f. It has been C updated since 1983, but I forget when and how. C G.M. May 1994 C*********************************************************************** PARAMETER (KB=40) COMMON/BLKT/TF(KB),T(KB),TB(KB),SF(KB),S(KB),SB(KB) COMMON/BLKU/UF(KB),U(KB),UB(KB),VF(KB),V(KB),VB(KB) COMMON/BLKQ/Q2F(KB),Q2(KB),Q2B(KB),Q2LF(KB),Q2L(KB),Q2LB(KB) COMMON/BLKR/RHO(KB),WT(KB) COMMON/BLKWU/WUSURF,WVSURF,WUBOT,WVBOT COMMON/BLKWT/WTSURF,WSSURF COMMON/BLKZ/Z(KB),ZZ(KB),DZ(KB),DZZ(KB),UMOL COMMON/BLKH/H,D,DT COMMON/BLKA/A(KB),C(KB),VH(KB),VHP(KB) COMMON/BLKK/KM(KB),KH(KB),KQ(KB),L(KB) COMMON/MODEL/GM(KB),GH(KB),SM(KB),SH(KB),KN(KB), 1 SPROD(KB),BPROD(KB),PROD(KB),DTEF(KB) REAL KAPPA,KM,KH,KQ,L DIMENSION UG(KB),VG(KB) DATA KAPPA/0.40/,EPS/1.E-10/,PI/3.14159/ C DATA ISTART/1/,IEND/481/,IPRINT/96/ C C********************************************************************* C INSERT USER CHOICE OF PROBLEM PARAMETERS. C SUBROUTINE DEPTH ESTABLISHES THE VERTICAL GRID. Z VARIES FROM C -1 TO 0 SO THAT H*Z IS THE DIMENSIONAL GRID. H*ZZ IS A STAGGERED C GRID ON WHICH ARE LOCATED U,V,T AND S. SEE SUBROUTINE DEPTH FOR C MORE DETAILS. C ********************************************************************* CALL DEPTH(Z,ZZ,DZ,DZZ,KB) H=80. TIME=0. DT1=15.*60. COR=0. COR=2.*PI/86400. GRAV=9.806 SMOTH=0.10 UMOL=1.0E-5 C****INSERT USER CHOICE OF INITIAL CONDITIONS************************* DO 5 K=1,KB UB(K)=0. C UB(K)=2.*EXP(-((-ZZ(K)-.5)/.20)**2) U(K)=UB(K) UF(K)=U(K) C VB(K)=-.50*(1.+ZZ(K)) VB(K)=0. V(K)=VB(K) VF(K)=V(K) UG(K)=U(K) VG(K)=V(K) WT(K)=0. TB(K)=24.+.05*ZZ(K)*H T(K)=TB(K) TF(K)=T(K) SB(K)=35. S(K)=SB(K) SF(K)=S(K) Q2B(K)=1.E-5 Q2(K)=Q2B(K) Q2LB(K)=1.E-5 Q2L(K)=Q2LB(K) Q2LF(K)=Q2L(K) L(K)=Q2L(K)/Q2(K) KM(K)=0. KH(K)=0. KQ(K)=0. 5 CONTINUE CALL DENS(T,S,ZZ,H,RHO,KB) WUBOT=0. WVBOT=0. C*********************************************************************** DT2=2.*DT1 DT4=4.*DT1 DAYI=DT1/86400. C C WRITE(6,15) ISTART,IEND,DT1 15 FORMAT(//' MODEL START UP ISTART =',I4,' IEND =',I4, 1 ' DT =',F5.0,' S',/) C C C*********************************************************************** C* * C* BEGIN TIME MARCH * C* * C*********************************************************************** DO 9000 INT=ISTART,IEND C PER=2*PI/COR PER=3600. TIME=TIME+DT1/PER C C*********************************************************************** C INSERT USER CHOICE OF GEOSTROHIC VELOCITY AND SURFACE FORCEING C AS FUNCTIONS OF TIME. C QSURF IS HEAT ABSORBED AT SURFACE (POSITIVE IS WARMING) IN W M-2. C QSWRAD IS SHORT WAVE RADIATION (CAN ONLY BE POSITIVE) IN W M-2 WHICH C PENETRATES THE SURFACE AND ATTENUATES ACCORDING TO THE ALGORIITHM C IN PROFT. C WUSURF=-.0004 WVSURF=0. QSURF=600.*.484 QSURF=0. QSWRAD=0. WTSURF=-QSURF/(1.026*4.186E6) SWRAD =-QSWRAD/(1.026*4.186E6) WSSURF=0. C*********************************************************************** C CALL PROFQ(DT2) DO 325 K=1,KB Q2(K)=Q2(K)+.5*SMOTH*(Q2F(K)+Q2B(K)-2.*Q2(K)) Q2B(K)=Q2(K) Q2(K)=Q2F(K) Q2L(K)=Q2L(K)+.5*SMOTH*(Q2LF(K)+Q2LB(K)-2.*Q2L(K)) Q2LB(K)=Q2L(K) 325 Q2L(K)=Q2LF(K) C C CALL PROFT(TF,TB,WTSURF,SWRAD,TSURF,1,DT2) CALL PROFT(SF,SB,WSSURF,0.,SSURF,1,DT2) DO 345 K=1,KB T(K)=T(K)+.5*SMOTH*(TF(K)+TB(K)-2.*T(K)) TB(K)=T(K) T(K)=TF(K) S(K)=S(K)+.5*SMOTH*(SF(K)+SB(K)-2.*S(K)) SB(K)=S(K) S(K)=SF(K) 345 CONTINUE C CALL DENS(T,S,ZZ,H,RHO,KB) C DO 380 K=1,KB-1 UF(K)=UB(K)+DT2*COR*(V(K)-VG(K)) 380 VF(K)=VB(K)-DT2*COR*(U(K)-UG(K)) CALL PROFU(DT2) CALL PROFV(DT2) DO 382 K=1,KB U(K)=U(K)+.5*SMOTH*(UF(K)+UB(K)-2.*U(K)) V(K)=V(K)+.5*SMOTH*(VF(K)+VB(K)-2.*V(K)) UB(K)=U(K) U(K)=UF(K) VB(K)=V(K) 382 V(K)=VF(K) C**** Begin diagnostic and print section IF(MOD(INT,IPRINT).NE.0) GO TO 9000 CALL MLDPTH(ZZ,T,KB-1,ZZMLD) ZZDMLD=-ZZMLD*H SPROD(1)=(ABS(WUSURF)+EPS)**1.5/KAPPA BPROD(1)=0. SPROD(KB)=(ABS(WUBOT)+EPS)**1.5/KAPPA BPROD(KB)=0. WRITE(6,500) TIME,ZZDMLD 500 FORMAT(/,' TIME =',F7.2,' HOURS MLD =',F7.2) WRITE(6,501) 501 FORMAT(1X,' DEPTH U V T S DEPTH 1 Q2 L KM KH RF SM SH 2 ') DO 550 K=1,KB ZZD=ZZ(K)*H ZD=Z(K)*H/10. RF=BPROD(K)/(SPROD(K)+1.E-20) WRITE(6,502) ZZD,U(K),V(K),T(K),S(K),ZD,Q2(K),L(K),KM(K), 1 KH(K),RF,SM(K),SH(K) 502 FORMAT(1X,F7.1,4(1PE10.2),F7.1,7(1PE10.2)) 550 CONTINUE WRITE(6,504) 0.,WUSURF,WVSURF,WTSURF,WSSURF 504 FORMAT(/,1X,F6.0,4(1PE10.2),' = SURFACE FLUXES') WRITE(6,505) H,WUBOT ,WVBOT,0.,0. 505 FORMAT(1X,F6.0,4(1PE10.2),' = BOTTOM FLUXES') 9000 CONTINUE STOP END *DECK DENS SUBROUTINE DENS(t,s,zz,dt,rho,kb) implicit double precision (a-h,o-z) real t,s,zz,dt,grav,rho DIMENSION t(KB),s(KB),rho(KB),zz(KB) DATA GRAV/9.806/ C C ......... THIS SUBROUTINE COMPUTES DENSITY-1............. C T = POTENTIAL TEMPERATURE C Mellor,G.L., 1991: An equation of state for numerical models of C Oceans and Estuaries. J.Atmos.and Oceanic Tech. 8, 609-611 C DO 1 K=1,KB-1 TR=T(K) SR=S(K) C Here, the (approximate) pressure is in units of decibars. P=-GRAV*1.025*ZZ(K)*DT*0.1 C RHOR = 999.842594 + 6.793952E-2*TR $ - 9.095290E-3*TR**2 + 1.001685E-4*TR**3 $ - 1.120083E-6*TR**4 + 6.536332E-9*TR**5 C RHOR = RHOR + (0.824493 - 4.0899E-3*TR $ + 7.6438E-5*TR**2 - 8.2467E-7*TR**3 $ + 5.3875E-9*TR**4) * SR $ + (-5.72466E-3 + 1.0227E-4*TR $ - 1.6546E-6*TR**2) *(ABS(SR))**1.5 $ + 4.8314E-4 * SR**2 C For shallow water the pressure dependency can be neglected C in which case it should also be omitted in PROFQ CR =1449.1+.00821*P+4.55*TR-.045*TR**2+1.34*(SR-35.) RHOR=RHOR + 1.E4*P/CR**2*(1.-0.2*P/CR**2) C RHO(K)=(RHOR-1000.)*1.E-3 1 CONTINUE C return end *DECK DEPTH SUBROUTINE DEPTH(Z,ZZ,DZ,DZZ,KB) DIMENSION Z(KB),ZZ(KB),DZ(KB),DZZ(KB) KL1=.3*KB KL2=KB-2 C*********************************************************************** C THIS SUBROUTINE ESTABLISHES THE VERTICAL RESOLUTION WITH LOG C DISTRIBUTIONS AT THE TOP AND BOTTOM AND A LINEAR 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*********************************************************************** C KL1=.2*KB 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 3 K=2,KL1 Z(K)=-DEL1*EXP(.693147*FLOAT(K-2)) 3 ZZ(K)=-DEL1*EXP(.693147*(FLOAT(K)-1.5)) DO 4 K=KL1,KL2 Z(K)=-(FLOAT(K)-CC)/BB 4 ZZ(K)=-(FLOAT(K)-CC+0.5)/BB DO 5 K=KL2,KB-1 Z(K)=(1.0-DEL2*EXP(.693147*FLOAT(KB-K-1)))*(-1.) 5 ZZ(K)=(1.0-DEL2*EXP(.693147*(FLOAT(KB-K)-1.5)))*(-1.) Z(KB)=-1.0 ZZ(KB-1)=-1.*(1.-DEL2/2.) ZZ(KB)=-1.*(1.+DEL2/2.) C DO 10 K=1,KB 10 Z(K)=-FLOAT(K-1)/FLOAT(KB-1) DO 11 K=1,KB-1 11 ZZ(K)=.5*(Z(K)+Z(K+1)) C DO 6 K=1,KB-1 DZ(K)=Z(K)-Z(K+1) 6 DZZ(K)=ZZ(K)-ZZ(K+1) RETURN END * DECK MLDPTH SUBROUTINE MLDPTH(ZZ,T,KB,ZZMLD) DIMENSION ZZ(KB),T(KB) DO 100 K=1,KB-1 IF(T(1).GT.(T(K)+.2)) GO TO 100 ZZMLD=ZZ(K)-(T(K)+.2-T(1))*(ZZ(K)-ZZ(K+1))/(T(K)-T(K+1)) 100 CONTINUE RETURN END *DECK PROFQ SUBROUTINE PROFQ(DT2) PARAMETER (KB=40) COMMON/BLKT/TF(KB),T(KB),TB(KB),SF(KB),S(KB),SB(KB) COMMON/BLKU/UF(KB),U(KB),UB(KB),VF(KB),V(KB),VB(KB) COMMON/BLKQ/Q2F(KB),Q2(KB),Q2B(KB),Q2LF(KB),Q2L(KB),Q2LB(KB) COMMON/BLKR/RHO(KB),WT(KB) COMMON/BLKWU/WUSURF,WVSURF,WUBOT,WVBOT COMMON/BLKWT/WTSURF,WSSURF COMMON/BLKZ/Z(KB),ZZ(KB),DZ(KB),DZZ(KB),UMOL COMMON/BLKH/H,D,DT COMMON/BLKA/A(KB),C(KB),VH(KB),VHP(KB) COMMON/BLKK/KM(KB),KH(KB),KQ(KB),L(KB) COMMON/MODEL/GM(KB),GH(KB),SM(KB),SH(KB),KN(KB), 1 SPROD(KB),BPROD(KB),PROD(KB),DTEF(KB) DIMENSION BOYGR(KB),CC(KB) REAL KM,KH,KQ,L REAL KAPPA,KN 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 SM/40*0.39/,SH/40*0.49/,GM/40*0.154/,GH/40*0.154/ DATA KAPPA/0.40/,SQ/0.2/,CIWC/1.0/ DATA GEE/9.806/ DATA SMALL/1.E-10/ DATA EPS/1.E-30/ C DH=H DO 100 K=2,KB-1 A(K)=-DT2*(KQ(K+1)+KQ(K) )*.5/(DZZ(K-1)*DZ(K) 1 *DH*DH) C(K)=-DT2*(KQ(K-1)+KQ(K) )*.5/(DZZ(K-1)*DZ(K-1) 1 *DH*DH) 100 CONTINUE C*********************************************************************** C * C THE FOLLOWING SECTION SOLVES THE EQUATION * C DT2*(KQ*Q2')' - Q2*(2.*DT2*DTEF+1.) = -Q2B * C * C*********************************************************************** CONST1=16.6**.6666667*CIWC C Boundary Conditions VH (1)=0.0 VHP(1)=SQRT(WUSURF**2+WVSURF**2)*CONST1 Q2F(KB)=.5*SQRT((WUBOT+WUBOT)**2 2 + (WVBOT+WVBOT)**2)*CONST1 C DO 101 K=1,KB-1 C Calculate pressure in units of decibars P=-GEE*1.025*ZZ(K)*DH*.1 CC(K)=1449.1+.00821*P+4.55*T(K) -.045*T(K)**2 $ +1.34*(S(K)-35.) CC(K)=CC(K)/SQRT((1.-.01642*P/CC(K)) $ *(1.-0.40*P/CC(K)**2)) 101 CONTINUE DO 102 K=2,KB-1 BOYGR(K)=GEE*(RHO(K-1)-RHO(K))/(DZZ(K-1)*DH) 1 +GEE**2*2.*1.025/(CC(K-1)**2+CC(K)**2) DTEF(K)=Q2B(K)*SQRT(Q2B(K))/(B1*Q2LB(K)+EPS) SPROD(K)=.25*KM(K)*((U(K)+U(K)-U(K-1)-U(K-1))**2 2 +(V(K)+V(K)-V(K-1)-V(K-1))**2) 3 /(DZZ(K-1)*DH)**2*CIWC**2 BPROD(K)=KH(K)*BOYGR(K) 102 PROD(K)=SPROD(K)+BPROD(K) C Sweep downward DO 103 K=2,KB-1 VHP(K)=1./(A(K)+C(K)*(1.-VH(K-1)) 1 -(2.*DT2*DTEF(K)+1.) ) VH(K)=A(K)*VHP(K) VHP(K)=(-2.*DT2*PROD(K) 1 +C(K)*VHP(K-1)-Q2B(K))*VHP(K) 103 CONTINUE C Sweep upward DO 104 K=1,KB-1 KI=KB-K Q2F(KI)=VH(KI)*Q2F(KI+1)+VHP(KI) 104 CONTINUE C*********************************************************************** C * C THE FOLLOWING SECTION SOLVES THE EQUATION * C DT2(KQ*Q2L')' - Q2L*(DT2*DTEF+1.) = -Q2LB * C * C*********************************************************************** C Boundary Conditions VH(1)=0. VHP(1)=0. Q2LF(KB)=0. C Sweep downward DO 110 K=2,KB-1 DTEF(K) =DTEF(K)*(1.+E2*((1./ABS(Z(K)-Z(1))+ 1 1./ABS(Z(K)-Z(KB))) *L(K)/(DH*KAPPA))**2) VHP(K)=1./(A(K)+C(K)*(1.-VH(K-1)) 1 -(DT2*DTEF(K)+1.)) VH(K)=A(K)*VHP(K) VHP(K)=(DT2*(-(SPROD(K)+E3*BPROD(K)) 1 *L(K)*E1)+C(K)*VHP(K-1)-Q2LB(K))*VHP(K) 110 CONTINUE C Sweep upward DO 111 K=1,KB-1 KI=KB-K Q2LF(KI)=VH(KI)*Q2LF(KI+1)+VHP(KI) 111 CONTINUE DO 112 K=2,KB-1 IF(Q2F(K).GT.SMALL.AND.Q2LF(K).GT.SMALL) GO TO 112 Q2F(K)=SMALL Q2LF(K)=SMALL 112 CONTINUE c DO 113 K=2,KB-1 c L(K)=Q2LF(K)/Q2F(K) C IF(L(K)**2*BOYGR(K).LT.-0.281*Q2F(K)) C 1 Q2LF(K)=Q2F(K)*SQRT(-0.281*Q2F(K)/BOYGR(K)) c113 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 APPROAHES 0.0288 L(1)=0. L(KB)=0. GH(1)=0. GH(KB)=0. DO 220 K=2,KB-1 L(K)=Q2LF(K)/Q2F(K) GH(K)=L(K)**2/Q2F(K)*BOYGR(K) 220 CONTINUE DO 230 K=1,KB GH(K)=MIN(GH(K),.028) SH(K)=COEF1/(1.-COEF2*GH(K)) SM(K)=COEF3+SH(K)*COEF4*GH(K) SM(K)=SM(K)/(1.-COEF5*GH(K)) 230 CONTINUE C DO 280 K=1,KB KN(K)=L(K)*SQRT(ABS(Q2(K))) KQ(K)=(KN(K)*.41*SH(K)+KQ(K))*.5 KM(K)=(KN(K)*SM(K)+KM(K))*.5 KH(K)=(KN(K)*SH(K)+KH(K))*.5 280 CONTINUE RETURN END *DECK PROFT SUBROUTINE PROFT(FF,FB,WFSURF,SWRAD,FSURF,NBC,DT2) PARAMETER (KB=40) COMMON/BLKT/TF(KB),T(KB),TB(KB),SF(KB),S(KB),SB(KB) COMMON/BLKU/UF(KB),U(KB),UB(KB),VF(KB),V(KB),VB(KB) COMMON/BLKQ/Q2F(KB),Q2(KB),Q2B(KB),Q2LF(KB),Q2L(KB),Q2LB(KB) COMMON/BLKR/RHO(KB),WT(KB) COMMON/BLKWU/WUSURF,WVSURF,WUBOT,WVBOT COMMON/BLKWT/WTSURF,WSSURF COMMON/BLKZ/Z(KB),ZZ(KB),DZ(KB),DZZ(KB),UMOL COMMON/BLKH/H,D,DT COMMON/BLKA/A(KB),C(KB),VH(KB),VHP(KB) COMMON/BLKK/KM(KB),KH(KB),KQ(KB),L(KB) REAL KM,KH,KQ,L REAL KAPPA DIMENSION FF(KB),FB(KB) DIMENSION TR(5),EXTC(5),RAD(KB) DATA KAPPA/0.4/,PR/1./ UMOLPR=UMOL 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----------------------------------------------------------------------- C NBC=1: SURF. B.C. IS WFSURF. NO RADIATIVE PENETRATION. C NBC=2; SURF. B.C. IS WFSURF+SWRAD*(1.-TR). WITH SWRAD*TR PENETRATION. C NBC=3; SURF. B.C. IS TSURF. NO SW RADIATIVE PENETRATION. C C NOTE THAT WTSURF (=WFSURF) AND SWRAD ARE NEG. VALUES WHEN WATER COLUMN IS C WARMING. C----------------------------------------------------------------------- DH=H DO 40 K=2,KB-1 A(K-1)=-DT2*(KH(K)+UMOLPR)/(DZ(K-1)*DZZ(K-1)*DH 1 *DH) C(K)=-DT2*(KH(K)+UMOLPR)/(DZ(K)*DZZ(K-1)*DH 1 *DH) 40 CONTINUE GO TO (50,51,52), NBC 50 CONTINUE C Boundary conditions and radiation flux calculation VH(1)=A(1)/(A(1)-1.) VHP(1)=-DT2*WFSURF/(-DZ(1)*DH)-FB(1) VHP(1)=VHP(1)/(A(1)-1.) DO 60 K=1,KB 60 RAD(K)=0. GO TO 53 C 51 CONTINUE VH(1)=A(1)/(A(1)-1.) VHP(1)=-DT2*(WFSURF+SWRAD*(1.-TR(NTP))) 1 /(-DZ(1)*DH)-FB(1) VHP(1)=VHP(1)/(A(1)-1.) DO 61 K=1,KB 61 RAD(K)=SWRAD*TR(NTP)*EXP(EXTC(NTP)*Z(K)*DH) GO TO 53 C 52 CONTINUE VH(1)=0. VHP(1)=FSURF DO 62 K=1,KB 62 RAD(K)=0. 53 CONTINUE C Calculate contribution due to radiation flux divergence DO 100 K=1,KB-1 100 FB(K)=FB(K)-DT2*(RAD(K)-RAD(K+1))/(DZ(K)*DH) C*********************************************************************** C THE FOLLOWING SECTION SOLVES THE EQUATION * C DT2*(KH*FF')' -FF = FB C*********************************************************************** C IF(MOD(INT,IPRINT).EQ.0) THEN C CALL PRXZ('RAD ',TIME,RAD,IM,2,JM,85,0,0,KB,0.,DT,Z) C ENDIF C Sweep downwards DO 101 K=2,KB-2 VHP(K)=1./(A(K)+C(K)*(1.-VH(K-1))-1.) VH(K)=A(K)*VHP(K) VHP(K)=(C(K)*VHP(K-1)-FB(K))*VHP(K) 101 CONTINUE C FF(KB-1)=FB(KB-1) C********************************************************************* C INSTEAD OF MATCHING SOLUTION TO A LOWER LAYER VALUE OF F(KB-1), C ONE MAY IMPOSE AN ADIABATIC BOTTOM B.C. BY REMOVING C'S FROM C COL. 1 OF THE NEXT TWO LINES. C********************************************************************* C FF(KB-1)=(C(KB-1)*VHP(KB-2)-FB(KB-1)) C 1 /(C(KB-1)*(1.-VH(KB-2))-1.) 99 CONTINUE C Sweep upwards DO 102 K=2,KB-1 KI=KB-K FF(KI)=VH(KI)*FF(KI+1)+VHP(KI) 102 CONTINUE RETURN END *DECK PROFU SUBROUTINE PROFU(DT2) PARAMETER (KB=40) COMMON/BLKT/TF(KB),T(KB),TB(KB),SF(KB),S(KB),SB(KB) COMMON/BLKU/UF(KB),U(KB),UB(KB),VF(KB),V(KB),VB(KB) COMMON/BLKQ/Q2F(KB),Q2(KB),Q2B(KB),Q2LF(KB),Q2L(KB),Q2LB(KB) COMMON/BLKR/RHO(KB),WT(KB) COMMON/BLKWU/WUSURF,WVSURF,WUBOT,WVBOT COMMON/BLKWT/WTSURF,WSSURF COMMON/BLKZ/Z(KB),ZZ(KB),DZ(KB),DZZ(KB),UMOL COMMON/BLKH/H,D,DT COMMON/BLKA/A(KB),C(KB),VH(KB),VHP(KB) COMMON/BLKK/KM(KB),KH(KB),KQ(KB),L(KB) REAL KM,KH,KQ,L C C*********************************************************************** C * C THE FOLLOWING SECTION SOLVES THE EQUATION * C DT2*(KM*U')' - U= -UB * C * C*********************************************************************** 85 DH=H DO 100 K=2,KB-1 A(K-1)=-DT2*(KM(K)+UMOL )/(DZ(K-1)*DZZ(K-1)*DH 1 *DH) C(K)=-DT2*(KM(K)+UMOL )/(DZ(K)*DZZ(K-1)*DH 1 *DH) 100 CONTINUE VH(1)=A(1)/(A(1)-1.) VHP(1)=(-DT2*WUSURF/(-DZ(1)*DH)-UF(1)) 1 /(A(1)-1.) DO 101 K=2,KB-2 VHP(K)=1./(A(K)+C(K)*(1.-VH(K-1))-1.) VH(K)=A(K)*VHP(K) VHP(K)=(C(K)*VHP(K-1)-UF(K))*VHP(K) 101 CONTINUE CBC=AMAX1(.0025,.16/ALOG((ZZ(KB-1)-Z(KB))*DH/.01)**2) CBC=CBC*SQRT(UB(KB-1)**2+(.25*(VB(KB-1) 1 +VB(KB-1)+VB(KB-1)+VB(KB-1)))**2) C******************************************************************** C TO RESTORE BOTTOM B.L. DELETE NEXT LINE C********************************************************************* CBC=0. UF(KB-1)=(C(KB-1)*VHP(KB-2)-UF(KB-1))/(CBC 1 *DT2/(-DZ(KB-1)*DH)-1.-(VH(KB-2)-1.)*C(KB-1)) DO 103 K=2,KB-1 KI=KB-K UF(KI)=VH(KI)*UF(KI+1)+VHP(KI) 103 CONTINUE 92 WUBOT=-CBC*UF(KB-1) RETURN END *DECK PROFV SUBROUTINE PROFV(DT2) PARAMETER (KB=40) COMMON/BLKT/TF(KB),T(KB),TB(KB),SF(KB),S(KB),SB(KB) COMMON/BLKU/UF(KB),U(KB),UB(KB),VF(KB),V(KB),VB(KB) COMMON/BLKQ/Q2F(KB),Q2(KB),Q2B(KB),Q2LF(KB),Q2L(KB),Q2LB(KB) COMMON/BLKR/RHO(KB),WT(KB) COMMON/BLKWU/WUSURF,WVSURF,WUBOT,WVBOT COMMON/BLKWT/WTSURF,WSSURF COMMON/BLKZ/Z(KB),ZZ(KB),DZ(KB),DZZ(KB),UMOL COMMON/BLKH/H,D,DT COMMON/BLKA/A(KB),C(KB),VH(KB),VHP(KB) COMMON/BLKK/KM(KB),KH(KB),KQ(KB),L(KB) REAL KM,KH,KQ,L C C*********************************************************************** C * C THE FOLLOWING SECTION SOLVES THE EQUATION * C DT2*(KM*V')' -V= -VB * C * C*********************************************************************** DH=H DO 100 K=2,KB-1 A(K-1)=-DT2*(KM(K)+UMOL )/(DZ(K-1)*DZZ(K-1)*DH 1 *DH) C(K)=-DT2*(KM(K)+UMOL )/(DZ(K)*DZZ(K-1)*DH 1 *DH) 100 CONTINUE VH(1)=A(1)/(A(1)-1.) VHP(1)=(-DT2*WVSURF/(-DZ(1)*DH)-VF(1)) 1 /(A(1)-1.) 98 CONTINUE DO 101 K=2,KB-2 VHP(K)=1./(A(K)+C(K)*(1.-VH(K-1))-1.) VH(K)=A(K)*VHP(K) VHP(K)=(C(K)*VHP(K-1)-VF(K))*VHP(K) 101 CONTINUE 104 CBC=AMAX1(.0025,.16/ALOG((ZZ(KB-1)-Z(KB))*DH/.01)**2) CBC=CBC*SQRT((.25*(UB(KB-1)+UB(KB-1) 1 +UB(KB-1)+UB(KB-1)))**2+VB(KB-1)**2) C******************************************************************** C TO RESTORE BOTTOM B.L. DELETE NEXT LINE C********************************************************************* CBC=0. VF(KB-1)=(C(KB-1)*VHP(KB-2)-VF(KB-1))/(CBC 1 *DT2/(-DZ(KB-1)*DH)-1.-(VH(KB-2)-1.)*C(KB-1)) DO 103 K=2,KB-1 KI=KB-K VF(KI)=VH(KI)*VF(KI+1)+VHP(KI) 103 CONTINUE 92 WVBOT=-CBC*VF(KB-1) RETURN END C******************************************************************** C END OF FILE C********************************************************************* C*********************************************************************