olnm - compute loading coefficients from global (netCDF) maps


USAGE:
olnm @ instruction-file [:<target>]
PURPOSE:
Compute loading coefficients from the global maps that are produced with gotsalm, gotsalsm or gotsalzm.
These netCDF-grids may cover the whole globe exactly (the standard assumption);
or a spherical quadrangle; or a grid that has been augmented with east and west margins
to avoid the hassle with modulus operations.

Every complex field is read from, and eventually written to, two files, amplitude and phase.
The file names carry `amp´ respectively `pha´ before their dot+extension ending. 

Also a range of other tasks can be executed on one or two complex*8 grids (see Commands). Range of grids
may include netcdf ocean tide
Our netcdf-files distinguish between "wet" and "dry" nodes, the latter signified by z-values 999.0
(not a very robust concept for the future - anyway...). Upon reading, these markers can be made to
set flags in an array (FLZ), 1 or 0, respectively. By default, we assume an all-through wet grid.
For those, flag checking during processing would render a curse; and indeed, a value of 17 in
FLZ(1,1) (a curse in Swedish ;-) signals that dry cells will be absent. Flag checking can be activated
(presently only in namelist; a command option may be coded into olnm.f later).

Spherical quadrangle: the circular connection must be de-activated (Namelist, zqopt='-C').
Augmented grid: olnm does not allow this feature; worse still: augmented grids with 1/16 x 1/16
resolution are too big for input. With lower resolution, the margins must be cut away
(Namelist, specify mcrop). 

Program source: ~/Oload/afor/p/m/olnm.f
INSTRUCTION FILE:
  1.  Namelist &forms
  2.  Open-file block
  3.  [site list specified here unless file unit 11 is opened]
  4.  Commands
1. Namelist &forms
Gray lines: currently insignificant / not used
Blue lines: can be modified in commands TIDE


Parameter name
Type
Explanation
Default

otmodel
char*48__
Name of the tide model
Replaces '%%%' in the name of the input grid files
(parameter otpath)
'Schwiderski'

mcode
char*8
Tide model code; a string that goes into the comment of a BLQ output file
'schw-'

tide
char*4
Darwin tide symbol
'M2'

otpath
char*64
Path to and file name's initial string of the input grid files  
'/home/hgs/Oload/
%%%/####_'

quse_flags
logical
distinguish wet and dry nodes; skip the latter while processing
.false.

fext




qbin




nsites
integer
Upper limit of number of sites that will be processed (eventhough more are being read in). Specify -1 to allow 10,000 
-1

comp(2)
char*4
'RADI' (radial) and 'PTAN' (tangential displacement potential)
Replaces '####' in input file names (parameter otpath)
Consider 'EAST' and 'NRTH' when using the PROFILE command.
'radi','ptan'

ocomp(2)
char*4
'RADI' and 'TANG', appears in comment lines in BLQ files
'RADI','TANG'

fmtin




phacut
real*4
For output of amplitude and phase, where to cut the phase circle in degrees (usually 0.0 or 180.0)
0.0

zqopt
char*8
Options in addition to ZQUINT_OPT('+G'), see zquint.f


nmsg
integer
Limit of messages from geocoord.f routines
9999

sclout
real*4
Factor to multiply amplitudes with. 1.0 means units of meter
1000.

blq_site_fmt
char*64
See subroutine put_ol_rec in getolrec.f
An empty string will revert to the standard format.
' '

opt_putolrec
char*16
See subroutine put_ol_rec in getolrec.f
Options to
put_ol_rec. An empty string defers to the subroutine's defaults.
' '

qfft_grad
logical
If .true., uses the Fourier method for gradient calculation.
.false.

monitor_grid_at(2)
integer
Output complex arrays, profile samples, before and after forward and backward transforms. Subroutine outsp is used, and the tags are CSA, CSB, CSC, and CSD, respectively. The east component is written to unit 91, the north to 92. (outsp source file is ~/sas/p/ufiosps.f )
-1,-1

monitor_iun
integer
With the GRAD_ commands, a value > 0 will initiate output of profile samples for east to this file unit, north to value+1
-1


smm_opt

char*32

Options for PTAN smoothing


lfssm
integer
Filter off-centre length
-1

issm_from
issm_to
integer
Column index range where to smooth
-1


landcount_tolerance

integer

For command COASTCUT

3

scissors_length
real*4
For command COASTCUT, unit = degrees
1.0


amp_putolrec





ltide




xcyl




mcrop
integer
Crop this number of columns from the east and west margins of the grid file.


kstop
integer
Debugging purpose:
1 : Program stops after reading the netCDF file
2 : ... after grid cropping
More stops could be placed in the main program.
99

 
2. Open-file block

------------------------------------------------------------------------------------------------------
File unit - data type -
Purpose etc.
------------------------------------------------------------------------------------------------------
    11    - ascii     - optional; if open, the sites list is read from this file.
    31    - ascii     - Output file; loading coeffcients (BLQ, HARPOS, simple or default)
    40    - ascii     - Input; if output format is BLQ, a file with the header comments, else ignored.
    89    - ascii     - Output of smoothing with option +RMS, produces an xyz-file ready to process
                        with xyz2grd .
  91 92   - binary    - output; see
monitor_grid_at (namelist) or PROFILE (command). 
------------------------------------------------------------------------------------------------------
Other file units are involved in specific command context; the unit numbers are most often selectable.

4. Commands
Environment variables are recognized,  ${VARIABLE}  or  ${VARIABLE:[default]}
To withdraw a variable, setenv VARIABLE  is not sufficient; you must  unsetenv VARIABLE
  --------- Syntax: -----------------------------------------------------------------------------------

<empty> comment
; comment
COMMAND options -
<empty> options continued

OUT format directive    - ascii, results for the sites in BLQ or another format, see below.

OUTD path/              - the gradient fields of PTAN are written to this path in netCDF.
                         
File names: EAST_amp.<tide> EAST_pha.<tide>  NRTH_amp.<tide>
NRTH_amp.<tide
                         
Sites will (ought) not be processed.

     The output statements above must appear before the executive commands following:

[ NOOP | DO | DIFF | GRAD_N | GRAD_E | LCOMB | PROFILE ] TIDE Darwin-symbol directives   


                        - NOOP is a placeholder instead of a command.

                        - If a command is not followed with a TIDE declaration, it will act on the
                          field(s) currently in storage. Thus, the first command of a series of
                          commands must either be preceded by a command-less TIDE declaration or
                          or must be followed by a TIDE declaration on the same line.
                          A TIDE line can be continued on the next line if it ends with `-´ or `/`

                        - with DO , will execute the loading calculation for the sites.
                          Symbol = tide's Darwin symbol; must follow the string TIDE with exactly one
                          blankspace inbetween. 

                          Not true - or?  TIDE commands without preceding DO can be repeated to e.g. set parameters
                          one-by-one. (Rather: use continuation character `/´ or `-´)

                          A concluding DO will execute the calculation.
                          If no parameter or directive occurs in subsequent TIDE lines, the
                          previous settings are kept.
                          The TIDE command effectuates reading of the netcdf-data sets.

                        - with DIFF computes the complex difference between two fields.
                          Example, first tide is stored in field 2, second in 1, in- and output components = RADI :

                          OUTD   ./HAMTIDEq-coarse/diff_hr_
                          NOOP TIDE M2  P=/home/hgs/Oload/HAMTIDEq/RADI_        M=${OTM}/- G=${OTMGRID:[F008]} C=${OTM} -
                                =>2 I=RADI
                          DIFF TIDE M2  P=/home/hgs/Oload/HAMTIDEq-coarse/RADI_ M=${OTM}/- G=${OTMGRID:[F008]} C=${OTM} -
                                =>1 N=RADI I=RADI
                          STOP
                          (note the continuation character `-´ at line endings)

                        - with GRAD_x will call one of the two gradient subroutines and output the resulting
                          grids to the directory that was specified in the OUTD command. See namelist qfft_grad.
                          Example:

                          GRAD_N TIDE M2  P=/home/hgs/Oload/%%%/####_ M=${OTM}/- G=${OTMGRID:[F014]} C=${OTM}
                          GRAD_E
                          which could have been split (and continued)
                          TIDE M2  P=/home/hgs/Oload/%%%/####_ M=${OTM}/- G=${OTMGRID:[F014]} C=${OTM} 
                          GRAD_N
                          GRAD_E
                          TIDE S2
                          GRAD_N
                          GRAD_E
                          ...

                        - with PROFILE , scans two profiles, one along a latitude ring, one along a meridian.
                          Binary output using outsp (~/sas/p/ufiosps.f); check the log-file for details,
                          use  splist filename X  for an overview of file's content.
                          This command has options:
                          +P        - process the second grid (usually the PTAN component).
                          U=iun     - for writing to a--preferably--opened file.
                          DIFF=S[X] - compute simple  differences, with X in the cross directions of the
                                      profiles.
                          DIFF=D[X] -
compute central differences (default is: no differencing).
                          S=scale   - a scale factor, reset after profile is processed, default=1.d0
                          X,Y=x,y   - longitude,latitude index numbers (default is monitor-grid_at from namelist).
                          E,N=
x,y   - longitude,latitude in degrees.

                        - with LCOMB computes linear combination
Z(,,3) = Z(,,1) + c·Z(,,2)
                          where c can be given a value using option SF=c (default = 1.0)


STOP                      Last line of input.


----- Tide directives:
-----------------------------------------------------------------------------------
P=path                  - Path to the grid file.
                          May contain a string `####´ which will be replaced with the component's name
                          (namelist variables COMP(1) and COMP(2) ) 
                          May contain a string `%%%´ which will be replaced by the ...

M=model-code[/-]          e.g. the name of a subdirectory containing grid files. Overrides mcode in namelist.
                          The `/-´ is necessary to avoid auto-names _amp.tide and _pha.tide
                          Example:
                            COMP='RADI','PTAN'   in namelist and 
                          TIDE M2  P=/home/hgs/Oload/%%%/####_ M=FES2014q/-
                         
will read the files
/home/hgs/Oload/FES2014q/RADI_amp.m2  _pha.m2
                                              /home/hgs/Oload/FES2014q/PTAN_amp.m2  _pha.m2
                         
, and the next TIDE command
                          TIDE S2
                         
will read the S2 files of kin.

G=grid-code             - obsolete

DO STORE n              - For BLQ type output, the column number of the tide (there should be a
                          nice default, but there isn't.
                          Issue has been handed over to the Mañana department). 


----- Format directives:
----------------------------------------------------------------------------------
                          For loading results, specify one of

BLQ

HARPOS

SIMPLE                  - writes to 31
(amp in mm, pha in degrees):
                          SITE lon lat TIDE radial east north

<empty>
                 - writes to 31 (amp in mm, pha in degrees):
                          SITE MONU TIDE RADI     radial
                          SITE MONU TIDE TANG     east    north

PTAN smoothing:
If a scalar field is noisy with crackles, they might get amplified when e.g. a gradient is calculated.
As a remedy, smoothing with a short low-pass filter can be applied. This is controlled through the namelist. Smoothing may be invoked in three different contexts, i.e. the command in action:
  1.  DO - (loading effects) smoothing is applied to the PTAN field but only in a narrow band of columns that depend on the site's longitudes.
  2.  GRAD_d - smoothing is applied on the PTAN field, which must be read as the second component (namelist parameter COMP(2)). Smoothing will not be applied if the starting column is less than zero (namelist parameter ISMM_FROM; if ISMM_FROM is set to zero, the full range of columns 1...m is smoothed). An upper limit can be set with namelist parameter ISMM_TO. If the filter length is less than zero (namelist parameter LFSMM), smoothing is not carried out.
  3.  PROFILE - The field that the PROFILE-command acts upon (+P given or omitted) is smoothed under thesame conditions as in the GRAD_d commands. 
The responsible subroutines are in Oload/afor/p/zsmoothms.f
 
SMM_OPT     
String of options for smoothing (namelist parameter, char*32):
+P       
 - print.
+DBG     
 - more extensive print, for debugging.
+F{R|T|C}  - one of three filter types: simple average ("Rectangle"), triangular weights, or
   raised-cosine. Off-centre length is determined by namelist parameter LFSMM .
+RMS[>t
 - Output on unit 89 those samples where the effect of smoothing exceeds threshold t.
   The output is formatted ready for reading with xyz2grd .
   Since PTAN is in units of meter-times-radian, the effect measure is converted to
   a more intuitive millimeter-times-degree. Default t = 0.1
+RMS[>t]=>u 
 - Write to file unit u. The previous unit, if open, will be closed. This prevents
   writing over a file that was opened through the open-block. Default u = 89
 
If the parameters SMM_OPT, ISMM_FROM, ISFMM_TO and LFSMM need to be changed on the fly, use command
SMOOTH [SMM_OPT=opt ] [ISMM_FROM_TO=i,j ] [LFSMM=len]
(note: the assigments must be coded with parameter names in upper-case and without intervening
blanks around `=´ and `,´)
For switching the feature off, issue  SMOOTH LFSMM=-1   or  SMOOTH OPT=OFF
OFF
(or any string not containing +F ) disables the smoothing filter in subsequent calls, causing
an early return from the subroutine.

For repeating the smoothing operation in another PROFILE or GRAD_d command on the same field,
use
SMOOTH RESET
e.g.

OUTD   ./FES2014qdg2/
TIDE M2  P=/home/hgs/Oload/%%%/####_ M=${OTM}/- G=${OTMGRID:[F014]} C=${OTM}
SMOOTH OPT=+RMS>0.1=>87
PROFILE U=91 +P E,N=-38.50,-4.00
SMO
OTH RESET
SMOOTH +RMS>0.01=>88
GRAD_N
SMOOTH +RMS=>89
GRAD_E
PROFILE U=92 +P E,N=-38.50,-4.00

The second PROFILE re-applies the smoothing parameters; however, the record of already smoothed
columns is cleared.

Coastline scissoring
... the purpose being production of grids of partial loading, removing wet cells within a specified distance from dry points and clean-cutting to obtain straight edges of the same measure. The controlling parameter (namelist) is scissors_length in degrees (henceforth L). The detection of land in a cell of  L x L with mostly wet nodes for flagging that tile for state change may allow a small number of dry points, namelist parameter landcount_tolerance.
The parameters can also be specified on the command line.

The nodes that change state from wet to dry can be written to an ASCII file in GMT xyz-style; a grd-file can then be created. 
   

The command sequence is
OUTD   dir_out/
COASTCUT [TIDE ...] [ =>1 ] [ L=scissors_length ][ T=landcount_tolerance ][ N=tag ][ MASK>iun ]
     Requires namelist specification of
       comp='....','SKIP'
     in order to input only one field, a field that does not carry the component in the file name.

     If the masking data is to be output, a file should be opened in the open-file block.

     EXAMPLE:
Prepare environment:
setenv SCISSORS 0.5
if ( x$1 != x ) setenv SCISSORS $1
setenv CUTL `calc -f "%2.2i" "16*$SCISSORS"`
Instruction file:
CC>
 &forms
 kstop=99
 mcrop=0
 otmodel='FES2014b'
 quse_flags=.true.
 nsites=-99
 comp='....','SKIP'
 &end
   C
71 B tmp/cc$CUTL.xyz
   Q

!
! Computing the difference of Radial between integrated and non-integrated Greens functions
!
OUTD   ./FES2014s/
COASTCUT TIDE M2  P=FES2014b M=FES2014b/- G=${OTMGRID:[F008]} C=FES2014b -
      =>1 L=${SCISSORS:[1.0]} N=CC${CUTL:[16]} MASK>71
STOP

     Note the continuation character at the end of the COASTCUT command line.
     The gray code is presently insignificant.   


EXAMPLE, A LOADING JOB FOR PRECISE SITE POSITIONS:
setenv OTM FES2014q
olnm @ olnm-OTM.ins :PENNA

setenv OTM <model> where model is the subdirectory name of a
                   set of loading maps, e.g. HAMTIDEq

OBS We have produced TPXO.7.1e/ with longitude margins of 4 (&forms below: mcrop=4)
OBS                  TPXO.7.1/

PENNA>
 &forms
 kstop=99
 mcrop=0
 otmodel='FES2014p'
 fmtin='(t6,a9,a5,t27,3f15.4)'
 zqopt='+C'
 nmsg=-1
 comp='radi','ptan'
 ocomp='RADI','TANG'
 opt_putolrec=' '
 &end
40 O ~/Oload/rray/header.blq
31 B ${OTM}/olnm_penna.rsl
11 * O sites.lonlat
   * instead of the site table below, active the above by removing '* '
   Q
 3 sites
-/ +P /
-/ EN FMT (a8,1x,a6,t21,2f10.0)
WTZR     WTZR          12.8774   49.1450
ALBH     ALBH         236.5126   48.3898
ALRT     ALRT         297.6596   82.4943
OUT SIMPLE
DO TIDE M2  P=/home/hgs/Oload/%%%/####_ M=${OTM}/- G=${OTMGRID:[F014]} C=${OTM} STORE 1
DO TIDE S2  STORE 2
DO TIDE N2  STORE 3
DO TIDE K2  STORE 4
DO TIDE K1  STORE 5
DO TIDE O1  STORE 6
DO TIDE P1  STORE 7
DO TIDE Q1  STORE 8
DO TIDE MF  STORE 9
DO TIDE MM  STORE 10
DO TIDE SSA STORE 11
STOP

EXAMPLE, A LINEAR-COMBINATION JOB
setenv OTM FES2014b-STW105-c-s1ni
setenv COMP RADI
setenv TIDE N2
exit

U>
 &forms
 kstop=99
 mcrop=0
 otmodel='FES2104q'
 nsites=-99
 comp='RADI','RADI'
 &end

40 O ~/Oload/rray/header.blq
31 B tmp/olnm.rsl
11 * O sites.lonlat
   Q
1 sites
-/ +P /
-/ EN FMT (a8,1x,a6,t21,2f10.0)
NEWC     NEWC         358.3834   54.9791
;
; Add coarse result and Machiel's corrections together
;
OUT BLQ

OUTD   ./FES2014b-s/${COMP}_

NOOP   TIDE ${TIDE}  P=FES2014b-STW105-c-s1ni/${COMP}_  M=${OTM}/- G=${OTMGRID:[F016]} C=${OTM} -
       =>1 I=${COMP}

LCOMB  TIDE ${TIDE}  P=FES2014b-c/${COMP}_ SF=0.001  +R M=${OTM}/- G=${OTMGRID:[F016]} C=${OTM} -
       =>2 I=${COMP}
STOP

Exec:
source olnm-add.ins
olnm @ olnm-add.ins :P

Note, +R is a highly special option; the two grids differed in the east-length by 1, whence the grid must be advanced by 1 every line in
subroutine ZGRID_LINCOMB (zgriddiffs.f)
.bye