Programs and utilities for
generation of an ocean tide model.
Catalog of modules in library
liboteq.a at the end of this document.
Utilities: c.f. OTEU.DOC,
chapters I and II.
-----------------------
º MAIN PROGRAM CREAM º
-----------------------
CREate A Model
CREAM requires an instruction
file CREAM.INS to be placed into a subdirectory,
the ModelSubDirectory. The file
is read via unit 4. The MSDir will hold all
data files related with a
specific area as modelling will proceed.
Job phases: The model creation
job consists of the following phases:
-----------
Phase Task
1a Take the >Area
Outline<, spherical quadrangle bounded by integer
latitudes/longitudes and project
it to the plane. The system
determines a plane that encloses
this region completely. Grid boxes
back-projecting outside the Area
Outline are normally flagged
"Land" (except under
>Tiled< >Entire< area.
System assembles a first
'Z'-flags array, the >Raw Flag array< from
a temporary depth array on the
'Z'-grid.
1b Assembling the
>Depth Array< on a temporary 'M'-grid.
Optional
interactive
inspection of the depth arrays.
1c Optional output
of Raw arrays: Goto 7 and stop
1d Alternative to
1a+b: Input of previously saved Raw Arrays.
2a User supplied
update of land/sea ('Z'-grid) and depth array
('M'-grid) using UPDFLM (input from file unit 4).
Four depth values are affected
for each flag value. If this is
undesired, depth may be updated
in phase 3 instead.
2b Graphic display of 'Z'-flags together with GMT
coastlines
using procedure GRA_GLCOAST (glcopls.f).
Opportunity to update flag array
interactively, graphic screen.
Update parameters are output to
the protocol file (unit 7); they
can be used in subsequent CREAM
jobs, file unit 4, under (2a).
3 User
supplied update (file unit 4) of depth array,
exact location on 'M'-grid.
4 Assembling
passive boundary condition flags. User supplied data
is processed (file unit 4 / CONsole) indicating the style of
connectivity of land / sea in
diagonal directions.
The user is prompted at
locations of connectivity ambiguities
eventually remaining.
The interaction is via a graphic
interface. The
conflicts
are shown
on a zoom-in map.
There are two passes where the user
can enter the type of
connection. In the first pass, the
passive boundary conditions have
been devised up to the
conflict point but not beyond.
Thus, the grid may look a
little difficult to interpret.
In the second pass the
situation is more clear-cut.
Decision can be deferred to the
second pass by answering D
(Defer) or A (defer All). Otherwise
the user may enter Y or N
depending upon whether the pair of
cells is to be sea-connected or
not.
The answers are collected in the
MSDIR/prflm.prt file.
From there they can be copied
and pasted into the
CONNCT section of
MSDIR/cream.ins. If that resolves all conflicts,
CREAM will run quicker next time
(unless the user updates
the grid under 2b).
5 Assembly of
active boundaries (user supplied data on file unit 4).
This means: Injection of unresolved flags for Active Boundaries
and Out-Of-Area nodes, 'Z'- and
'M'-grids.
6 Display of
depth array using procedure GraSol;
interactive update of depth array using graphic screen:
Enter: Answer at GraSol-prompt
C x
where x > 0
Exit: Set STOP (enter S-key) at
GraSol-prompt.
7 Output
The
printed output
The file connected to unit 7 is
used for PRFLM, printing interpreted flag
arrays. 'M'-flags representing
land bridges are indicated by \ or /,
channels by u or v, for NW-SE or
NE-SW connections, respectively.
Passive boundary symbols are
intuitive.
Example:
Side remark:
Example is Denmark. You realize
that an opening exists between mainland and
Fyn, whereas Store Baelt is
closed ! There is a lot of manual updating needed.
Baltic sea is to be closed off as
the area is too close to the right border
of the model.
--------------------
Side remark 2
There are opportunities to update
the depth array produced under (1b) at later
stages: Use ...
(1) interactive feature in TTEQ
/ OTEMT1. C.f. OTET.doc, OTES12.doc
(2) main program UPDEPTH.exp
(4) rerun Prep.2 phase,
OTEM92u.exp
(3) subroutine POKER, oteu013.f
File MSDir/cream.ins Structure
and codes
=========================================================================
Namelist &MODEL parameters:
_________________________________________________________________________
XLB,XUB,YLB,YUB (real*4):
Lower and Upper Bounds of the area:
X in
longitude, Y in latitude. Values in degrees.
SCALE
(real*4): Grid spacing in
kilometers.
SHORE, TOPMAX, HRESX, HRESY,
CFRSEA (real*4): c.f. routine ASSDEP
TOPODSC (char*64): Name of the
topo-bathy-descriptor file (default
/home/hgs/TOPO/topo.dsc) e.g.
/home/hgs/ETOPO1/topo.dsc
Qdbg_TOPO
.true.: extensive printing for each call of CTOPO
QASS
(logical): C.f. above, job
phases.
.true.: All job phases
.false.: Input "raw" arrays and do job phases 2...7
QSTOP
(logical): Stop after reading the
topo files (and QLOOK_H_).
QLOOK_H_Z
(log.) Interactive inspection of the depth array on
the
'Z'-grid.
QLOOK_H_M
(log.) Interactive inspection of the depth array on
the
'M'-grid.
QEntire_Area
(log.) The entire area is used (default =
.false.)
By default the area declared by XLB,XUB... determines the
exploitable parts of arrays. Due
to the projection of that
spherical area, dead land areas
will result along the edges
of the plane grids.
QTiled_Area
(log.) (default=.false.) under QEntire_Area=.true.
Outside the area outline, some 1x1 degree grid cells on the
sphere project partly outside
the plane grid. Use the Tiled-
area option to preserve only
"entire tiles"; grid boxes in
broken tiles will be flagged
"land". Yet, the amount of
"dead" land is much less than
with QEntire_Area=.false. The
area edges will be rugged as
they will line up with integer
meridians and latitude circles.
N.B.: This option must be backed
up with an (excessive)
pile of "Keep_Global_Land" alt.
"Keep_Global_Ocean" areas
in order to prohibit double
account of tide loading effects.
The Keep_Global_.. COMMON area
is presently too small for
this.
QShow_Coast
(log.) (default=.true.)
Show the coastline when the flag array is display for the
first time.
TEMPDIR
(char*64) (default=' ')
The directory to which GMT coastline dumps are made. An empty
string defaults to using
/home/hgs/GLCOAST. This catalog is
not open to strangers; they need
to create their own temporary
directory. There is a scratch
commons though. Use
TEMPDIR='/scratch/hgs/glcoast-dump/'
QSKIP_SHOW_CPROBLEM (log.)
(default=.salse.)
Don't show / don't prompt for decisions concerning
connection problems during pass
1 of Passive-boundary
detection at the stage of M-grid
assembly.
Instead defer decision to to
pass 2.
=========================================================================
INSTRUCTION INPUT
_________________________________________________________________________
Max 80 characters of comment.
_________________________________________________________________________
Namelist &MODEL
_________________________________________________________________________
File open block for input of
"raw" arrays (QASS=.false.), respectively
file open block for output of
"raw" arrays or final arrays (QASS=.true.)
file unit data type input/output
31 - 'Z'-flags and /CAREA/ I O
32 - 'M'-flags and depth I O
33 - Active boundary locations
(if phase 5 is nonempty) O
_________________________________________________________________________
File open block for output
(processed if QASS=.false., skipped if .true.)
_________________________________________________________________________
Saving raw flag arrays and exit:
Instruction: SAVE RAW
The system continues with phase 7
alt.:
Instruction: GRAPHIC
The system continues with phase 6
alt.:
Update sets for the Z-flag array
(transparent reading):
Instruction (1) UPDFLM
Instruction (2) N, skip-flag. N =
number of Update-instructions
until next instruction of type
(2).
Skip-flag Meaning
'T' Take
'S' Skip (=ignore)
'E' End of update sets.
Instructions (3..): Update
information to UPDFLM:
I1,J1,I2,J2,Flagvalue,Depth,'comment'
c.f. keyword UPDFLM below
_________________________________________________________________________
Updating the depth array at exact
location on 'M'-grid:
Instruction (1) UPDPTH
Instruction (2) N, skip-flag. N =
number of update-instructions
until next instruction of type
(2)
Instructions (3..): Update
information:
I,J,Depth,'comment'
_________________________________________________________________________
Connectivity information for
pairs of corners
Instruction (1) CONNCT
Instruction (2..) Y i,j - yes,
connect sea
N i,j - no, land bridge
(n) T - terminate CONNCT
_________________________________________________________________________
Injecting active boundaries:
Instruction (1) BOUNDZ
Instructions (2..) Boundary
parameter sets
c.f. keyword ABOUND below
_________________________________________________________________________
These instructions are difficult
(impossible) to anticipate. Most efficient
processing is:
(1) to execute with QASS=.true.
and SAVE RAW
(2) to execute with QASS=.false.,
retrieving the raw arrays, and using
trivial instruction sets under
UPDFLM, UPDPTH, CONNCT and ABOUND
(3) to iterate (2), building up
the instruction sets successively.
Hints are found in the protocol
file MSDIR/prflm.prt
Example: Med-Sea.
24 km grid interval, longitude
range 10 W... 37 E, latitude range 29 N ... 47 N
10x10 grid cell refinement to
interpolate depth
66% sea area required for a box
to obtain "Sea" status.
Do all job phases at once.
____________________________
START OF EXAMPLE ________________________________
Prep.step.1 of model OT12:
MEDITERRANEAN (new)
&MODEL
XLB=-10., XUB=37., YLB=29.,
YUB=47., SCALE=24.
HRESX=10., HRESY=10.
SHORE=-0.5,
TOPMAX=9.,
CFRSEA=-0.667
QASS=.true.,
qstop=.false.
&END
1 * Output files
case
QASS=.TRUE., input files case QASS=.FALSE.
7 B prflm.prt Protocol file
31 B MDL12/FLZRAW.DAT
32 B MDL12/FLMHRAW.DAT
33 B MDL12/ABOUND.DAT
Q
1 * This block is
bypassed due
to QASS=.TRUE.
31 B MDL12/FLZ.DAT
32 B MDL12/FLMH.DAT
33 B MDL12/ABOUND.DAT
Q
UPDFLM
9, 'T'
1,21, 3,26,0,-9999.,
'Mauretanian coast'
67,30, 75,33,0,-9999., 'Algerian
lowland'
126,10,183,16,0,-9999., 'Libyan
lowlands etc.'
182,23,183,33,0,-9999.,
'Israelian lowland'
183,55,183,55,0,-9999., 'Turkeian
lowland'
144,62,180,94,0,-9999., 'Black
Sea'
24,46, 24,46,1, 100., 'Strait of
Gibraltar'
24,45, 25,45,1, 200., 'Strait of
Gibraltar'
23,44, 23,44,1, 100., 'Strait of
Gibraltar'
25, 'T'
103,48,103,48,0,-9999., 'Messina,
Land'
104,47,104,47,0,-9999., 'Messina,
Land'
128,50,130,50,0,-9999., 'Korinth'
140,46,141,46,0,-9999., 'Kyklades'
138,48,138,48,0,-9999., 'Kyklades'
143,49,143,49,0,-9999., 'Kyklades'
92,78, 92,79,1, 3., 'Comacchio'
97,83, 97,83,1, 3., 'Trieste'
32,40, 32,40,1, 25., 'Melilla'
128,44,128,44,1, 75., 'Kalamai
Peloponnes'
182,60,128,60,1, 15.,
'Thessaloniki'
142,56,142,56,1, 15., 'Izmir'
143,37,143,37,0,-9999., 'Kreta'
141,36,141,36,0,-9999., 'Kreta'
135,38,135,38,0,-9999., 'Kreta'
131,43,131,43,0,-9999.,
'Peloponnes'
174,43,174,43,0,-9999., 'Cyprus'
149,45,149,45,0,-9999., 'Rhodos'
150,46,150,46,0,-9999., 'Rhodos'
149,42,149,42,0,-9999., 'Rhodos'
123,48,123,48,0,-9999.,
'Zakinthos Ionian'
123,50,123,50,0,-9999., 'Levkas
Ionian'
54,55, 54,56,0,-9999.,
'Ibiza+Formentera'
61,58, 61,58,0,-9999., 'Menorca'
87,43, 87,43,0,-9999., 'C.Bon,
Tunisia'
1, 'T'
101,38,101,38,0,-9999., 'Malta'
0, 'E'
UPDPTH
0, 'E'
CONNCT
N 86,28
N 87,43
N 150,43
N 174,43
N 150,46
Y 104,48
N 105,48
Y 123,49
Y 123,50
N 124,51
N 146,51
Y 133,53
N 144,53
Y 142,56
Y 143,56
N 141,60
Y 99,80
Y 100,80
T
ABOUND
W 20,64,17 "GIBR"
O 1,65,49,95
T
______________________________
END OF EXAMPLE ________________________________
ASSDEP: PARAMETERS and some
comments
============================
SUBROUTINE ASSDEP (H,FLM,TYPE,M,N,DEPTH,SHORE
&, TOPLND,RESX,RESY,CFRSEA,iuw)
Model construction * Make Depth
and Flag arrays
ASSDEP creates depth array H and
flag array FLM by calling the
EXTERNAL topography-function
DEPTH from north to south.
Border points (outside the area
determined by AOUTLI) will obtain
flag=100,000. Otherwise,
flag=0 for topography above
SHORE,
flag=1 below SHORE.
TYPE:
(char*1) An 'M' or a 'Z' type depth matrix (and accompanying
flag matrix) can be produced.
H(M,N): The
bathymetry grid (negative topography) is returned in H(M,N).
SHORE:
Topography values higher than this value (in meters)
are
taken as "land"; this affects the flag array and optionally
the H(i,j) values on "land".
A reasonable value is always
less than zero.
TOPLND: The value
which is to limit H(i,j) in case "land".
TOPLND
is of type topography (not bathymetry).
TOPLND
> 1E6 or < -1E6: ASSDEP will return H as the full
topographic
or bathymetric field.
RESX &
Resolution for obtaining mean depth values (in km).
RESY
If values of less or equal 0. or greater than SCALE (= the
grid
spacing) are used, default = SCALE will be taken,
which
means that the depth value
at the box centre is used.
Make sure
that RESX/Y is an
integer divisor of SCALE.
(SCALE is
specified to
subroutine AOUTLI.)
CFRSEA: Critical
fraction of sea content in a cell; discriminates
whether a
cell is counted "sea" or "land". Applies only
if -1E6
< TOPLAND < 1E6,
i.e. in "not-full topo/bathy" mode.
Specify a
value greater than 1
in order to obtain sea-land
discrimination on the basis of
average topography.
If CFRSEA
< 0 is specified,
the absolute of CFRSEA value is taken;
however,
average water depth is
only computed for the water
fraction
of the cell.
If
CFRSEA=0 is specified, one of
100 cells = sea is enough
to declare
the whole cell = sea;
pure land will be preserved
as land
however.
The
definition of SHORE is
applied after CFRSEA.
IUW
Protocol print log.unit
For use in preparation of OTEQ
depth and flag arrays, it is recommended to
(1) use CFRSEA = -.5 ... -.7
(2) SHORE = -10.
with the effect that average
depth is computed for the water-covered
fraction; if this depth is too
shallow, the cell is flagged "land".
Call Entire_Area (m,n) OTES811.f
prior to ASSDEP
to resolve depth array also
outside projected model area.
Then, ITOZ will require a wide
Grace limit !
UPDFLM: SPECIFYING UPDATE
PARAMETERS FOR FLAG ARRAYS
============================================
i1,j1,i2,j2,Flagvalue,Depth,'comment'
Set flag value in array from
south-west corner i1,j1 to north-east i2,j2
and set bathymetry to specified
depth value. The comment is printed on
the protocol.
Changing bathymetry only: Specify
Flagvalue=-9999
Changing flag value only: Specify
Depth=-9999.9
N.B.: CREAM.f calls UPDFLM on the
'Z'-flags and 'M'-depth arrays. This
implies that the four neighbour
nodes of bathymetry are subject
to change at each processing of
a flag node. In this mode, the
convention is that UPDFLM can
only increase depth.
If changes of bathymetry are
required at exact places, section
UPDPTH must be used.
Legal flag values:
Y/N Y/N indicates applies
to 'Z'-flags and 'M'-flags, respectively.
0
Land
Y Y
1
Sea
Y Y
10*b
Passive boundary
N Y
100,000
O_o_area
unresolved Y Y
-n
" storage ref.
Y Y
200,000
Act.bound
unresolved Y Y
200,000+10*b
" with pass.bound.cond. N Y
200,000-n " storage ref.
Y Y
n Storage reference, < 9000
b Pass.bound.cond. number
ABOUND: SPECIFYING ACTIVE
BOUNDARIES
============================
1. Summary
-----------
Processing a boundary = resolving
storage references at flags ("A")
P i1,j1,i2,j2 Protect a strip
when processing the next boundary.
The protected strip is limited
by the straight
line which connects (i1,j1) and
(i2,j2).
Boundary type N and S: i1 <
i2 is required,
" " E " W: j1
< j2
x i1,i2,j "NAME" where x = { N S
}: define north or south boundary from
i
= i1 to i2 at j = const.
x j1,j2,i "NAME" where x = { E W
}: define west or east boundary from
j
= j1 to j2 at i = const.
"NAME" is a unique name for each
boundary, enclosed in
double quotes. If no name is
given, a default name
"00##", ## = running number as a
string, is supplied.
The name can be used in TTEQ to
tweek
the boundary values (search for
_UPDAB).
New feature since Feb 2011.
Example
N 3,85,67 "NRTH"
goes from 3 to 85 eastward at j-position 67
E 3,85,67 "WEST"
goes from 3 to 85 northward at i-position 67
O
i1,j1,i2,j2 Change "Sea" into
"Out-of-area" in the rectangle
(i1,j1)
= Southwest corner, (i2,j2) = Northeast.
D n n >
0: Delete the last n
active boundaries
n <
0:
Step
up boundary index by n
T
Terminate
reading.
2. Codes and cases
------------------
(C1) Simple case.
Specify
x i,j,l
x = boundary orientation, i,j =
range from, to; l = location.
Boundary orientation = { E N S W
} (East, North ...)
Examples:
N 12,40,62 "NRTH"
meaning: Northern boundary, from
I=12 to 40 (eastward) at north 62.
Range and location are array
indices.
W 3,70,11 "WEST"
meaning: Western boundary, from
J=3 to 70 (northward) at east 11.
The Skagerak-Kattegat model
(KATT/cream.ins) runs such a simple case;
ABOUND
W 100, 180, 10 "WEST"
T
(C2) Difficult geometry.
ABOUND will turn "Open
Sea"-flags into "Out-of-area"-flags in the
region west of an active
W-boundary, east of an E-boundary etc.
Example:
........XXX,,,...XXXX
27
..........XXX,...XXXX
26
............XXX...XXX
25
..............aX...XX
24
..............A.....X
23
..............A.....X
22
Type N act.bound.
AAAAAAAAAAAAAAA...... 21
.....................
20
|
|
|
30
40
50
When the boundary is resolved,
the sea cells (".") to the north of it are
changed to "O", which affects
the model area at the locations marked ","
above.
In order to protect a model area
against the flag change, introduce
"protection polygons" before
defining the boundary:
P ib1,jb1,ie1,je1
P ib2,jb2,ie2,je2
... .... (max 10 polygone
sections can be handled).
A polygon section protects the
array cells behind it, seen from the
active boundary's vantage. It is
a line connecting start point (ibn,jbn)
and end point (ien,jen).
In the example above,
P 39,27,45,23
would be appropriate.
(C3) Left-over "Sea" cells to be
turned into "Out-of-area" cells.
ABOUND does not declare cells
"Out-of-area" if they lie north or south
of the end points of an 'E' or
'W' boundary, or - analogously - if they
lie west or east of the end
points of an 'N' or 'S' boundary. These cells
must be changed into "Out-of
area" by explicit instructions:
O i1,j1,i2,j2
Change "Sea" into "Out-of-area" in the rectangle
(i1,j1)
= Southwest corner, (i2,j2) = Northeast.
(C4) An active boundary can be
introduced temporarily.
Protection polygones can be
attached to a temporary boundary, i.e. a
boundary that is deleted again
after processing.
This facilitates the conversion
of "Sea" to "O-o-area" cells in those
cases were many O-instructions
would be needed.
3. Split jobs. First job
terminates after ABOUND.
Next job continues with ABOUND
and ABFLZ.
--------------------------------------------------
CREAM does not split these
phases. If reprogrammed, observe that...
The instructions for (C2) and
(C3) are passed to ABFLZ (entry in MKFLZ)
through common /CABOUP/.
The boundary location data (1) is
stored in common /CABOUN/. It can be
in-/output using GETAB / OUTAB.
The instruction D -1 can be used
after a set of P-instructions to
store protection polygones
without processing of any active boundary.
3. Examples
------------
(E1) The third boundary is to be
protected.
Instructions to ABOUND with
processing:
N 1,30,60 "N001"
N 50,55,65 "N002"
P 50,25,51,30
P 51,30,45,40
E 12,40,45 "EAST"
T
(E2) The third boundary is to be
protected.
Instructions to ABOUND without
processing:
Assume that three boundaries
have been declared earlier. They are already
represented in the flag arrays.
ABFLZ is to be rerun.
D 9999 Don't process any boundary
D -2 Assign next instruction(s)
to the 3'rd boundary
P 50,25,51,30 Introduce
protection polygones
P 51,30,45,40
D -1 Keep information for three
boundaries
T
Contents of library:
liboteq.a, context cream
(This is not really up-to-date)
Module:
otes41
Bathymetry,
internal routines
_depth_
fct. Bathymetry at a
random point
_extcdp_ fct. entry Extend buffer for
depth
Module:
otes71
Passive
boundaries, more or less obsolete.
_boundt_
subr. Bound.cond.numbers from
'M'-flags
_testsu_
" Used only
by BOUNDT
_mirr_
"
"
" " "
_rot_
"
"
" " "
Module: otes72
_boundz_iuo entry Output unit
for e.g. copying CONNCT requests.
_boundz_
subr. Bound.cond.numbers from
'Z'-flags; preferred.
_connct_
" Declare
connection status for land-land diagonals
Module:
otes81
Model
assembly, user's calls
_assdep_
Make
the bathymetry array
_aoutli_
subr. Define model area
_updflm_
" Update flag and depth array
_mkflz_
" Make 'Z'-flags from 'M'-flags
_abflz_
" Mark active boundary in 'Z'-flag
array
_abound_
" Mark active boundary and
determine 'O_o_area'
Module: otes811
_bentar_ block data /centar/
_entire_area_ subr. Depth resolved in
entire area (default: model outline)
_tiled_area_ entry Cut area borders
at integer degrees of sphere.
_enq_entar_ subr. Enquire
Module:
otes42
Bathymetry
file routines
_rtopo_
subr. Transfer data to a
spherical grid area.
_rtopo_e_ entry
...with alternate exit for EOF
_ctopo_
entry Continue transfer
_ctopo_e_ entry
... ... EOF
Module: glcopls
_gra_glcoast_ subr. Update flag array on
graph.screen. Disspla coasts
_cursor_upd_ subr. Cursor move and
prompt
.bye
Graphic
displays
2b GRA_GLCOAST
Above: This is the display of the
raw Z-flag array. After pressing U for cursor-Upd, ...
... a few changes have been made:
Yellow squares have been inserted to flip status to "Land" by
left-clicking a box and pressing L,
light-blue boxes to flip status
to "Sea". Finally, a rectangle has been started (yellow frame) to
introduce "Sea" into
southern Öresund by first
left-clicking on a box and pressing B to mark one corner, then
left-clicking on a box to mark
the diametric corner. Next S
should be pressed to indicate "Sea" as the new status. Then, the
Return-key will bring
the user back to the first menu bar, where A for accept can be pressed.
Or U to continue with updating (Sjællands Odde e.g.
should be turned into a chain of
land-boxes for instance).
There is no way to regret single
edits. The work done is not lost since the edits are written to the
protocol file
(prtflm.prt) from where they can
be copied-and-pasted into the instruction file (cream.ins). There is
nothing to get distroyed
either, since the raw flag array
is kept unchanged; thus, you can start all over.
Back to stages, 2b.
4 Diagonal dams and
straights
Defer = don't decide now.
Defer1 = this one problem.
The program can be configured
such that the first pass through the array is not shown. In the second
pass, the ambiguous connections
should all be resolved. To do this iteratively (no necessity, but still
there's this option) verifiable
decisions can be copied-and-pasted from the PRFLM.PRT file into the
cream.ins file. If the last
iteration of cream (successive runs of cream) is successful, PREP2 can
be
started (otem92).
Back to Phase 4