CCM Test: Salinity propagation in Oresund Channel.
  **************************************************************
  DISPLAY
  ----------------------------------------------------------
  This case concerns simulation of a stratified flow through
  Oresund Channel. Originally, problem had been set by Urban
  Svensson (CFE). The library case is its simplified version.
 
  The  main  objective  of  this  case is to demonstrate new
  features of CCM specially  developed for simulation of  so
  called shallow water flows, i.e flows with character  size
  in  one  direction  much  less  than  these in two others.
  'Shallow' direction is Z axis. (See MBFGE entry in PHENC).
 
  There are two options:
   LSG10 when TRUE activates simulation of shallow flows  by
         truncated 3D Navier-Stokes equations;
   LSG8  when TRUE  activates same as LSG10  treatment; plus
         pressure-correction equation is solved 2D, with WC1
         defined by integration of continuity equation.
  ----------------------------------------------------------
  ENDDIS
l(pause
  **************************************************************
BOOLEAN(LTURB,LSHWT,SHWT3Z)
CCM= T;  LTURB= T;  SHWT3Z= T;  LSHWT= F
  **************************************************************
   IGRID : 1 =>  7x13xNZ grid
           2 => 14x27xNZ grid
           3 => 29x55xNZ grid
  **************************************************************
INTEGER(NTIM,IGRID); IGRID= 2
REAL(REFDEN,ALFA,SALT0,SALTIN,PRSALT,DTS)
REFDEN= 1000.; ALFA= 8.E-4; SALT0= 10.; SALTIN= 30.; PRSALT= 2.
NTIM  = 10;   DTS = 600.
  **************************************************************
  PHOTON USE
   p ; ; ; ; ;
 
   set vec av off
   msg Velocity Vectors:
   vi z
   vec k 1 sh
   pause
   cl
   vec k m sh
  ENDUSE
  **************************************************************
    GROUP 1. Run title and other preliminaries
IF(CCM) THEN
+IF(LSHWT) THEN
+ IF(LTURB) THEN
+  TEXT( CCM: Bay flow (2D,K-E)
+ ELSE
+  TEXT( CCM: Bay flow (2D,Lam)
+ ENDIF
+ELSE
+ IF(LTURB) THEN
+  TEXT( CCM: Bay flow (3D,K-E)
+ ELSE
+  TEXT( CCM: Bay flow (3D,Lam)
+ ENDIF
+ENDIF
ELSE
+IF(LTURB) THEN
+ TEXT(PHOE: Bay flow (3D,K-E)
+ELSE
+ TEXT(PHOE: Bay flow (3D,Lam)
+ENDIF
ENDIF
TITLE
    GROUP 6. Transience; time-step specification
STEADY= F;  GRDPWR(T,NTIM,NTIM*DTS,1.0)
SPEDAT(SET,GXMONI,TRANSIENT,L,F)
    GROUP 6. Body-fitted coordinates or grid distortion
BFC= T;  NONORT= .NOT.CCM; READCO(GRID5)
    GROUP 7. Variables stored, solved & named
NAME(H1)= SALT; SOLVE(P1,V1,U1,W1,SALT);
STORE(VPOR,DUDX,DUDY,DUDZ,DVDX,DVDY,DVDZ,DWDX,DWDY,DWDZ)
SOLUTN(P1,Y,Y,Y,N,N,N)
IF(CCM) THEN
+ NAME(C1)= UC1;  NAME(C2)= VC1;  NAME(C3)= WC1
+  SOLVE(UC1,VC1,WC1);      SOLUTN(UC1,Y,Y,Y,P,P,P)
+ SOLUTN( VC1,Y,Y,Y,P,P,P); SOLUTN(WC1,Y,Y,Y,P,P,P)
+ SOLUTN(SALT,Y,Y,Y,P,P,P)
+ TERMS(U1,N,N,N,N,N,N); TERMS(UC1,N,Y,Y,P,P,P)
+ TERMS(V1,N,N,N,N,N,N); TERMS(VC1,N,Y,Y,P,P,P)
+ TERMS(W1,N,N,N,N,N,N); TERMS(WC1,N,Y,Y,P,P,P)
+ FIINIT(UC1)= 1.E-6;  FIINIT(VC1)= 1.E-6;  FIINIT(WC1)= 0.0
ENDIF
IF(LTURB) THEN
+ TURMOD(KEMODL);  STORE(ENUT,GEN1,RHO1,LEN1)
+ SOLUTN(KE,Y,Y,Y,P,P,P);  SOLUTN(EP,Y,Y,Y,P,P,P)
+ FIINIT(KE)= 1.E-4;  FIINIT(EP)= 1.E-6
+ ENUL= 1.E-6;  FIINIT(ENUT)= 1.E-3
ELSE
+ ENUL= 1.E-6
ENDIF
    GROUP 8. Terms (in differential equations) & devices
TERMS(SALT,N,Y,Y,P,P,P)
    GROUP 9. Properties of the medium (or media)
RHO1= LINSCAL; RHO1A= REFDEN; RHO1B= ALFA*REFDEN
PRNDTL(SALT)= 1000.;  PRT(SALT)= PRSALT
    GROUP 11. Initialization of variable or porosity fields
INIADD= F;  FIINIT(SALT)= SALT0
    *** Model VPOR setting for 7x13xNZ grid:
IF(IGRID.EQ.1) THEN
    ***** Saltholm:
+ CONPOR(0.0,CELL,4,4, 5, 7,1,NZ)
    ***** Helsingor:
+ CONPOR(0.0,CELL,1,2,13,13,1,NZ)
+ CONPOR(0.0,CELL,1,1,12,12,1,NZ)
ENDIF
    *** Model VPOR setting for 14x27xNZ grid:
IF(IGRID.EQ.2) THEN
    ***** Saltholm:
+ CONPOR(0.0,CELL,7,7,10,14,1,NZ)
    ***** Helsingor:
+ CONPOR(0.0,CELL,1,4,26,27,1,NZ)
+ CONPOR(0.0,CELL,1,3,25,25,1,NZ)
+ CONPOR(0.0,CELL,1,2,24,24,1,NZ)
+ CONPOR(0.0,CELL,1,1,23,23,1,NZ)
    ***** Hven:
+ CONPOR(0.0,CELL,8,8,21,21,1,NZ)
ENDIF
    *** VPOR setting for main grid 29x55x16:
IF(IGRID.EQ.3) THEN
    ***** Falsterbo:
+ CONPOR(0.0,CELL,25,29,5,5,1,NZ)
+ CONPOR(0.0,CELL,25,25,6,6,1,NZ)
    ***** Klagshamn:
+ CONPOR(0.0,CELL,29,29,9,9,1,NZ)
    ***** Saltholm:
+ CONPOR(0.0,CELL,13,13,19,28,1,NZ)
+ CONPOR(0.0,CELL,14,14,19,20,1,NZ)
    ***** Hven:
+ CONPOR(0.0,CELL,16,17,42,42,1,NZ)
    ***** Helsingborg:
+ CONPOR(0.0,CELL,28,29,47,47,1,NZ)
+ CONPOR(0.0,CELL,29,29,48,48,1,NZ)
    ***** Helsingor:
+ CONPOR(0.0,CELL,1, 1,47,47,1,NZ)
+ CONPOR(0.0,CELL,1, 3,48,48,1,NZ)
+ CONPOR(0.0,CELL,1, 7,49,49,1,NZ)
+ CONPOR(0.0,CELL,1,10,50,51,1,NZ)
+ CONPOR(0.0,CELL,1, 9,52,55,1,NZ)
ENDIF
    GROUP 13. Boundary conditions and special sources
    ***** INLET of salted water is defined from 1 to NZ/2
IF(IGRID.EQ.1) THEN
    ** North Boundary:
+PATCH( NORR,NORTH,3,NX,NY,NY,NZ/2+1,  NZ,1,LSTEP)
+PATCH(NORRL,NORTH,3,NX,NY,NY,     1,NZ/2,1,LSTEP)
    ** South Boundary:
+PATCH(  SYD,SOUTH,3,NX, 1, 1,1,NZ,1,LSTEP)
ENDIF
IF(IGRID.EQ.2) THEN
    ** North Boundary:
+PATCH( NORR,NORTH,5,NX,NY,NY,NZ/2+1,  NZ,1,LSTEP)
+PATCH(NORRL,NORTH,5,NX,NY,NY,     1,NZ/2,1,LSTEP)
    ** South Boundary:
+PATCH(  SYD,SOUTH,5,NX, 1, 1,1,NZ,1,LSTEP)
ENDIF
IF(IGRID.EQ.3) THEN
    ** North Boundary:
+PATCH( NORR,NORTH,10,NX,NY,NY,NZ/2+1,  NZ,1,LSTEP)
+PATCH(NORRL,NORTH,15,24,NY,NY,     1,NZ/2,1,LSTEP)
    ** South Boundary:
+PATCH(  SYD,SOUTH, 9,NX, 1, 1,1,NZ,1,LSTEP)
ENDIF
 COVAL(NORR, P1,10.*FIXP,0.); COVAL(NORR,SALT,ONLYMS,SAME)
 COVAL(NORRL,P1,1.E-20,1.E20*REFDEN*0.3)
 COVAL(NORRL,SALT,ONLYMS,SALTIN)
 COVAL(SYD,P1,1000.,1000.); COVAL(SYD,SALT,ONLYMS,SALT0)
    ** Bottom Friction:
PATCH(BOTFRI,LWALL,1,NX,1,NY,1,1,1,LSTEP)
IF(LTURB) THEN
+ COVAL(   SYD,KE,ONLYMS,1.E-4); COVAL(   SYD,EP,ONLYMS,1.E-6)
+ COVAL( NORRL,KE,ONLYMS,1.E-4); COVAL( NORRL,EP,ONLYMS,1.E-6)
+ COVAL(  NORR,KE,ONLYMS,1.E-4); COVAL(  NORR,EP,ONLYMS,1.E-6)
+ COVAL(BOTFRI,KE, LOGLAW,LOGLAW); COVAL(BOTFRI,EP, LOGLAW,LOGLAW)
+IF(CCM) THEN
+ COVAL(BOTFRI,UC1,LOGLAW,0.0); COVAL(BOTFRI,VC1,LOGLAW,0.0)
+ COVAL(BOTFRI,WC1,LOGLAW,0.0)
+ELSE
+ COVAL(BOTFRI,U1,LOGLAW,0.0); COVAL(BOTFRI,V1,LOGLAW,0.0)
+ COVAL(BOTFRI,W1,LOGLAW,0.0)
+ENDIF
   ** Turbulence/Buoyancy interaction:
+ PATCH(KEBUOY,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP)
+  COVAL(KEBUOY,KE,GRND4,GRND4); COVAL(KEBUOY,EP,GRND4,GRND4)
ELSE
+IF(CCM) THEN
+ COVAL(BOTFRI,UC1,1.0,0.0); COVAL(BOTFRI,VC1,1.0,0.0)
+ COVAL(BOTFRI,WC1,1.0,0.0)
+ELSE
+ COVAL(BOTFRI,U1,1.0,0.0); COVAL(BOTFRI,V1,1.0,0.0)
+ENDIF
ENDIF
    ** Buoyancy forces:
BUOYC= -9.81
IF(CCM) THEN
+ BUOYE= REFDEN*(1.+ALFA*SALT0)
ELSE
+ PATCH(BUOY,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP)
+ BUOYD= REFDEN*(1.+ALFA*SALT0); COVAL(BUOY,W1,FIXFLU,DENSDIFF)
+ COVAL(BUOY,U1,FIXFLU,DENSDIFF);COVAL(BUOY,V1,FIXFLU,DENSDIFF)
ENDIF
    ** Coriolis forces:
CORIOL= 1.E-4
    GROUP 15. Termination of sweeps
LSWEEP= 50;  TSTSWP= -1
    GROUP 16. Termination of iterations
SELREF= T;  RESFAC= 1.E-6
    GROUP 17. Under-relaxation devices
RELAX(P1,LINRLX,0.5); RELAX(SALT,FALSDT,1000.)
IF(LTURB) THEN
+ KELIN= 3; RELAX(KE,LINRLX,0.25); RELAX(EP,LINRLX,0.25)
+ LITER(KE)= 10; LITER(EP)= 10; OUTPUT(ENUT,Y,N,N,Y,Y,Y)
ENDIF
    GROUP 18. Limits on variables or increments to them
VARMIN(SALT)= SALT0;  VARMAX(SALT)= SALTIN
VARMAX(ENUT)= 1.0
    GROUP 19. Data communicated by satellite to GROUND
IF(CCM) THEN
+ LSG4= T; LSG8= LSHWT; LSG10= SHWT3Z
+ RELAX(UC1,FALSDT,1000.); RELAX(VC1,FALSDT,1000.)
+ RELAX(WC1,FALSDT,   1.)
+ LITER(P1)= 20; LITER(UC1)= 5; LITER(VC1)= 5; LITER(WC1)= 5
ELSE
+ RELAX(U1,FALSDT,1000.); RELAX(V1,FALSDT,1000.)
+ RELAX(W1,FALSDT,   1.)
ENDIF
    GROUP 22. Spot-value print-out
IXMON= 7;  IYMON= 5;  IZMON= NZ/2+1
    GROUP 23. Field print-out and plot control
IDISPA= 10; CSG1=URB