TALK=F;RUN(1,1) DISPLAY This Kelvin-Helmholz simulation can be run in the following modes: 'pure k-e' or 'pure MFM' ---------- ---------- 1. Effective viscosity from (a) k**2 / e or (b) mixl * mnsq 2. Micro-mixing rate from (a) e / k or (b) mnsq / mixl 3. Micro-mixing length from (e) k**1.5 / e or (b) source = const * mnsq Eight combinations are obviously possible. One or other of the 'pure' cases can be set by way of a single boolean variable. This Q1 comprises that created by Dr Hornby for the k-epsilon model, with addition of MFM features. It requires the use of a special ground (see appendix D). ENDDIS boolean(pureke,puremfm) pureke=t puremfm=.not.pureke GROUP 1. Run title and other preliminaries ** TANH PROFILES HAVE MAXIMUM INFLEXION AT THE C/L * SHEAR MODEL REFERENCE * 1 LINEAR SHEAR * 2 TANH SHEAR PROFILE (THICKNESS PARAMETER GTSH) REAL(GYF) GYF=20.0 INTEGER (JSH) JSH=2 IG(1)=JSH * TEMPERATURE MODEL REFERENCE (THICKNESS PARAMETER GTTS) * 1 LINEAR PROFILE * 2 TANH PROFILE INTEGER (JTEM) JTEM=2 IG(2)=JTEM * SALINITY MODEL REFERENCE (THICKNESS PARAMETER GTTS) * 1 LINEAR PROFILE * 2 TANH PROFILE INTEGER (JSAL) JSAL=2 IG(3)=JSAL * DIMENSIONS REAL(GX1,GY1) * X DOMAIN GX1=0.16 * Y DOMAIN GY1=0.1 BOOLEAN(JTURB,JRES) INTEGER(JTMDL) * JTURB SOLVE JTURB=T IF(JTURB)THEN * TURBULENCE MODEL * 1 LES * 2 K-E JTMDL=2 IF(JTMDL.EQ.1)THEN TEXT(LES:KH instability) ELSE TEXT(K-EP:KH instability) ENDIF ENDIF * JRES RESTART FLAG JRES=F REAL(GCPW,GRHOW,KCOND1,GDIFH,GDIFS) REAL(GHS,GHB,GHREF,GTREF,GSS,GSB,GSREF,GPRTH,GPRTS) REAL(GTS,GTB) REAL(GBT,GBS,G) REAL(GUIN,GSU) * WATER SPECIFIC HEAT GCPW=4180. * TEMPERATURE GTS=12.0 GTB=10.0 GHS=GCPW*GTS GHB=GCPW*GTB GHREF=0.5*(GHB+GHS) GTREF=0.5*(GTB+GTS) * SALINITY GSS=3.0 GSB=218.0 GSREF=0.5*(GSB+GSS) * CONSTANTS FOR DENSITY * Joint Panel on Oceanographic Tables and Standards (UNESCO 1981) * Formulae simplified by assuming zero pressure * Valid for depths ~1km REAL(A1,A2,A3,A4,A5) REAL(A6,A7,A8,A9,A10) REAL(A11,A12,A13,A14,A15) A1=999.842594 A2=6.793952E-2 A3=-9.095290E-3 A4=1.001685E-4 A5=-1.120083E-6 A6=6.536332E-9 A7=0.824493 A8=-4.0899E-3 A9=7.6438E-5 A10=-8.2467E-7 A11=5.3875E-9 A12=-5.72466E-3 A13=1.0227E-4 A14=-1.6546E-6 A15=4.8314E-4 REAL(G1,G2,G3,GRHOS,GRHOB) G1=A1+A2*GTREF+A3*GTREF**2+A4*GTREF**3+A5*GTREF**4+A6*GTREF**5 G2=GSREF*(A7+A8*GTREF+A9*GTREF**2+A10*GTREF**3+A11*GTREF**4) G3=GSREF**1.5*(A12+A13*GTREF+A14*GTREF**2) GRHOW=G1+G2+G3+A15*GSREF**2 G1=A1+A2*GTB+A3*GTB**2+A4*GTB**3+A5*GTB**4+A6*GTB**5 G2=GSB*(A7+A8*GTB+A9*GTB**2+A10*GTB**3+A11*GTB**4) G3=GSB**1.5*(A12+A13*GTB+A14*GTB**2) GRHOB=G1+G2+G3+A15*GSB**2 G1=A1+A2*GTS+A3*GTS**2+A4*GTS**3+A5*GTS**4+A6*GTS**5 G2=GSS*(A7+A8*GTS+A9*GTS**2+A10*GTS**3+A11*GTS**4) G3=GSS**1.5*(A12+A13*GTS+A14*GTS**2) GRHOS=G1+G2+G3+A15*GSS**2 GRHOB GRHOW GRHOS G1=A2+2*A3*GTREF+3*A4*GTREF**2+4*A5*GTREF**3+5*A6*GTREF**4 G2=GSREF*(A8+2*A9*GTREF+3*A10*GTREF**2+4*A11*GTREF**3) G3=GSREF**1.5*(A13+2*A14*GTREF) GBT=(G1+G2+G3)/GRHOW GBT G1=A7+A8*GTREF+A9*GTREF**2+A10*GTREF**3+A11*GTREF**4 G2=1.5*GSREF**0.5*(A12+A13*GTREF+A14*GTREF**2) G3=2*GSREF*A15 GBS=(G1+G2+G3)/GRHOW GBS * USE FOR TESTS grhow=1000.0 gbt=-1e-3 gbs=2e-3 grhow gbt gbs * WATER THERMAL CONDUCTIVITY KCOND1=0.6 * WATER MOMENTUM DIFFISIVITY ENUL=1E-6 * HEAT DIFFUSIVITY GDIFH=KCOND1/(GRHOW*GCPW) * SALT DIFFUSIVITY GDIFS=1E-9 * TURBULENT PRANDTL NUMBER OF HEAT GPRTH=0.7 * TURBULENT PRANDTL NUMBER OF SALT GPRTS=0.7 * GRAVITY G=9.81 RG(3)=GRHOW RG(4)=GCPW RG(5)=GHREF RG(6)=GBT RG(7)=GHS RG(8)=GHB RG(9)=G RG(10)=GBS RG(11)=GSS RG(12)=GSB RG(14)=GSREF GROUP 2. Transience; time-step specification REAL(GTIME,JNT) STEADY=F NREGT=2 IREGT=1; GRDPWR(T,1,10.0,1.0) IREGT=2; GRDPWR(T,1,10.0,1.0) JNT=600 GTIME=JNT*5.0/2000.0 GRDPWR(T,JNT,GTIME,1.0) GROUP 3. X-direction grid specification CARTES=T REAL(GDX1) REAL(JNX1) JNX1=100 JNX1=20 GDX1=GX1/JNX1 RG(16)=GDX1 NREGX=1 IREGX=1; GRDPWR(X,JNX1,GX1,1.0) XULAST XCYCLE=T GROUP 4. Y-direction grid specification REAL(GDY1,GPOW) INTEGER(JNY1,JNY0) JNY1=140 JNY1=30 GDY1=GY1/JNY1 GPOW=1.2 RG(2)=GPOW JNY0=JNY1/2 G1=GY1/2 NREGY=2 IREGY=1; GRDPWR(Y,JNY0,G1,-GPOW) IREGY=2; GRDPWR(Y,JNY0,G1,GPOW) YVLAST GROUP 5. Z-direction grid specification NREGZ=1 IREGZ=1; GRDPWR(Z,1,1.0,1.0) ZWLAST GROUP 6. Body-fitted coordinates or grid distortion GROUP 7. Variables stored, solved & named SOLVE(P1,U1,V1,H1) SOLVE(C1,C2) NAME(16)=SALT * Y in SOLUTN argument list denotes: * 1-stored 2-solved 3-whole-field * 4-point-by-point 5-explicit 6-harmonic averaging SOLUTN(P1 ,Y,Y,Y,N,N,Y) * Stored variables list STORE(TMP1,DEN1) IF(JTURB)THEN IF(JTMDL.EQ.1)THEN TURMOD(SGSMOD) * SPECIFY LENGTH SCALE USING GRID (EL1B=0.0) * OTHERWISE SET EL1B=0.41 (VON KARMAN'S CONSTANT) EL1B=0.0 * SMAGORINSKY CONSTANT EL1A=0.17 STORE(ENUT,WDIS) ELSE TURMOD(KEMODL) KELIN=3 STORE(ENUT) ENDIF ENDIF GROUP 8. Terms (in differential equations) & devices * Y in TERMS argument list denotes: * 1-built-in source 2-convection 3-diffusion 4-transient * 5-first phase variable 6-interphase transport TERMS(H1,N,P,P,P,P,P) ! cut out built-in sources, i.e. kinetic ! heating GROUP 9. Properties of the medium (or media) PRNDTL(H1)=ENUL/GDIFH IF(JTURB)THEN PRT(H1)=GPRTH ENDIF PRNDTL(SALT)=ENUL/GDIFS IF(JTURB)THEN PRT(SALT)=GPRTS ENDIF RHO1=GRND1 *USING ENTHALPY RELATIVE TO GHREF AND SALINITY RELATIVE TO GSREF *RHO=RHO1A+RHO1B*NAME(IBUOYB)+RHO1C*NAME(IBUOYC) IBUOYA=0 *ENTHALPY IBUOYB=14 *SALT(C1) IBUOYC=16 RHO1A=GRHOW RHO1B=GBT*GRHOW/GCPW RHO1C=GBS*GRHOW *CORRECT FOR USE OF ENTHALPY RELATIVE TO GHREF TMP1=GRND2 TMP1A=GHREF/GCPW TMP1B=1.0/GCPW * SHEAR TANH PROFILE THICKNESS REAL(GTSH) IF(JSH.EQ.1)THEN GTSH=YVLAST ELSE GTSH=YVLAST/GYF ENDIF RG(15)=GTSH * TEMPERATURE/SALINITY TANH PROFILE THICKNESS REAL(GTTS) IF(JSH.EQ.1)THEN GTTS=YVLAST ELSE GTTS=YVLAST/GYF ENDIF RG(17)=GTTS * BUOYANCY FREQUENCY REAL(GBF) GBF=G*(GRHOB-GRHOS)/(GRHOW*GTTS) GBF=SQRT(GBF) GBF *RESIDUAL FLOW GUIN=0.3 RG(1)=GUIN * SHEAR FREQUENCY REAL(GSH) GSH=2*GUIN/GTSH GSH *CFL REAL(GCFL,GDX,GDT) GDT=GTIME/JNT GDX=GX1/JNX1 GCFL=GUIN*GDT/GDX GCFL GROUP 10. Inter-phase-transfer processes and properties GROUP 11. Initialization of variable or porosity fields INIADD=F REAL (GENUT,GEP,GKE) GENUT=1E-10 GENUT=1E-5 GEP=1E-10 G1=(GENUT*GEP)/0.09 GKE=SQRT(G1) GKE GEP FIINIT(KE)=GKE ;FIINIT(EP)= GEP *INITIAL VELOCITY PERTURBATION FACTOR REAL(GVPER) GVPER=0.0025/GUIN RG(18)=GVPER IF(JRES)THEN RESTRT(ALL) ELSE PATCH(SOL0,INIVAL,1,NX,1,NY,1,NZ,1,1) INIT(SOL0,U1,0.000E+00,GRND) INIT(SOL0,V1,0.000E+00,GRND) INIT(SOL0,H1,0.000E+00,GRND) INIT(SOL0,SALT,0.000E+00,GRND) ENDIF fiinit(enut)= 1.e-4 ! these are not active fiinit(ke) = 1.e-3 fiinit(enut)= 1.e-3 GROUP 12. Patchwise adjustment of terms (in differential equations) GROUP 13. Boundary conditions and special sources * NUMERICAL SCHEME (DEFAULT TO UPWIND)...... * (KOREN BEST HIGHER ORDER SCHEME) SCHEME(KOREN,ALL) ! these are not active SCHEME(SUPBEE,ALL) SCHEME(CUS,ALL) SCHEME(SMART,ALL) * PRESSURE REF PATCH(PREF,CELL,1,1,1,1,1,1,1,#NREGT) COVAL(PREF,P1,FIXP,0.0) * CONCENTRATION SOURCE PATCH(TRACE,CELL,1,NX,NY/2,NY/2+1,1,NZ,1,#NREGT) COVAL(TRACE,C2,FIXVAL,1.0) * wall patches for wall distance IF(JTURB)THEN IF(JTMDL.EQ.1)THEN PATCH(WALL1,SWALL,1,NX,1,1,1,NZ,1,#NREGT) * LOG LAW COVAL(WALL1,U1,GRND2,0.0) PATCH(WALL2,NWALL,1,NX,NY,NY,1,NZ,1,#NREGT) * LOG LAW (NOTE A/C TAKEN OF INFLOW) COVAL(WALL2,U1,GRND2,0.0) ENDIF ENDIF * BUOYANCY SOURCE FOR V1..................... * 3.2 setting (full density) REAL(GALPH) * Tilt angle GALPH=(2.0/180.0)*3.141593 BUOYA=-G*SIN(GALPH);BUOYB=-G*COS(GALPH);BUOYC=0.0;BUOYD=GRHOW PATCH(BUOY,PHASEM,1,NX,1,NY,1,NZ,1,#NREGT) COVAL(BUOY,U1,FIXFLU,GRND2) COVAL(BUOY,V1,FIXFLU,GRND2) * BUOYANCY SOURCE FOR K IF(JTMDL.EQ.2)THEN PATCH(KEBUOY,PHASEM,1,NX,1,NY,1,NZ,1,#NREGT) COVAL(KEBUOY,KE,GRND4,GRND4) * C3E ZERO, 0.2 FOR STABLE STRATIFICATION * C3E=1.0 FOR UNSTABLE STRATIFICATION COVAL(KEBUOY,EP,GRND4,GRND4) SPEDAT(SET,KEBUOY,GCEB,R,0.2) ENDIF GROUP 14. Downstream pressure for PARAB=.TRUE. GROUP 15. Termination of sweeps LSWEEP=5 SELREF=T RESFAC=1E-3 TSTSWP=-1 GROUP 16. Termination of iterations GROUP 17. Under-relaxation devices real(GUMX,gdtf) GUMX=GUIN GDTF=0.2*XULAST/(NX*GUMX) RELAX(P1,LINRLX,0.7) RELAX(U1,FALSDT,GDTF) RELAX(V1,FALSDT,GDTF) RELAX(W1,FALSDT,GDTF) RELAX(H1,FALSDT,GDTF) RELAX(SALT,FALSDT,GDTF) IF(JTMDL.EQ.2)THEN RELAX(KE,LINRLX,0.7) RELAX(EP,LINRLX,0.7) RELAX(ENUT,LINRLX,0.5) ENDIF GROUP 18. Limits on variables or increments to them GROUP 19. Data communicated by satellite to GROUND GROUP 20. Preliminary print-out GROUP 21. Print-out of variables GROUP 22. Spot-value print-out IXMON=NX/2;IYMON=NY/2;IZMON=1 GROUP 23. Field print-out and plot control CSG1=A; ! deactivation of this line creates a 3D parphi file IDISPA=50 ! which allows the devlopment of the flow to be ! viewed in a single photon session * Y in OUTPUT argument list denotes: * 1-field 2-correction-eq. monitor 3-selective dumping * 4-whole-field residual 5-spot-value table 6-residual table IF(JTURB)THEN OUTPUT(ENUT ,Y,Y,Y,Y,Y,N) ENDIF OUTPUT(SALT ,Y,Y,Y,Y,Y,N) GROUP 24. Dumps for restarts !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! MFM section starts here INTEGER(NFLUIDS,NFLR,NFLF) CHAR(MFMMOD) ! declare NFLF=10; NFLR=1 ! choose NFLR for one-D population NFLUIDS=NFLF*NFLR REAL(VISCON,LENCON,CONMIX) ! declare necessary constants MFMMOD=MNSQ; ! calculate and use RMS of population LENCON=0.5; VISCON=20.0; CONMIX=5.0 ! model & constants ! see also enut = $ ! $MNSQ patch CONMIX=1.0 ! a non-optimised value GROUP 7. Variables stored, solved & named SOLVE(MIXL) ! solve and STORE(LEN1,RATE,MNSQ,AVEF) ! store needed quantities INTEGER(III) ! activate fluids settings III=NFLUIDS+1 DO II = 1,NFLUIDS ! DO loop order chosen so that III = III-1 ! F1 is solved first SOLVE(F:iii:) VARMAX(F:iii:)=1.0; VARMIN(F:iii:)=0.0; RELAX(F:iii:,LINRLX,0.5) ENDDO GROUP 11. Initialization of variable or porosity fields FIINIT(LEN1)= 0.00001*YVLAST ! it is important to initialise to FIINIT(MIXL)= 0.00001*YVLAST ! reasonable values PATCH(LOWER,INIVAL,1,NX,1,NY/2,1,1,1,1) ! the lower fluid is fluid 1 COVAL(LOWER,F1,0.0,1.0) PATCH(UPPER,INIVAL,1,NX,NY/2+1,NY,1,1,1,1) ! the upper fluid is fluid 10 COVAL(UPPER,F:NFLUIDS:,0.0,1.0) GROUP 13. Boundary conditions and special sources if(puremfm) then the following statements dictate that there is a source of MIXL, per unit volume, equal to LENCON*MNSQ PATCH($MNSQ,PHASEM,1,NX,1,NY,1,NZ,1,1) ! note use of 'dollar patch' COVAL($MNSQ, MIXL, LENCON*1.0E-5, 1.E5)! source of mixl ENUT=GRND10 else ! pure ke (STORED of MIXL IS len1 ) ! get mixl from k-e's length (STORED of RATE IS EPKE ) ! get rate frim ep/ke endif GROUP 18. Limits on variables or increments to them VARMAX(MNSQ)=1.0; VARMIN(MNSQ)=0.0 VARMAX(AVEF)=1.0; VARMIN(AVEF)=0.0 ! physical-realism demands GROUP 19. Data communicated by satellite to GROUND SPEDAT(MFM,MFMMOD, C,:MFMMOD:) SPEDAT(MFM,NFLUIDS,I,:NFLUIDS:) SPEDAT(MFM,NFLR, I,:NFLR:) SPEDAT(MFM,NFLF, I,:NFLF:) SPEDAT(MFM,CONMIX, R,:CONMIX:) SPEDAT(MFM,VISCON, R,:VISCON:) SPEDAT(MFM,POPMIN, R,-0.3) SPEDAT(MFM,POPMAX, R,+0.3) GROUP 23. Field print-out and plot control CHAR(NAMPROF) ! line-printer plot output DO II=1,NFLUIDS NAMPROF=PROF:II: PATCH(:NAMPROF:,PROFIL,nx/2,nx/2,1,NY,1,NZ,1,1) COVAL(:NAMPROF:,F:II:,0.0,0.0) ENDDO lsweep=3 STOP