Home

The CASCO4b User's Guide

1. 6 J 6f b t 0hr t 72 hr _ 4 J _ 4 2 2pa b c d a 2la Cc Cc oa E 0H Z Z oO 2 7 o 2 F ul ul 4 J 4 1 6 L L L a 6 L i T T T T 6 liz T T a L t 216 hr E aus Mr 5 5 2ta Cc Cc 2 OF E S o 2f ul ul _4 L L L L 6 a L L L 1 21 41 61 1 21 41 61 Boundary Node Number Boundary Node Number Figure 15 Inverse solution boundary conditions dots for Mode A with velocity forcing Also shown are the truth BC s line Initial conditions are zero for both truth and inverse solutions The a b c d represent the location of the corners of the mesh as shown in Fig 9 31 8 8 6 6 z 3 8 O gt a 2 2 2 0 0 0 5 10 15 0 5 10 15 1 1 5 3 0 5 3 l o o D 2 0 6 6 6 6 6 60 6 0 6 60 0 0 6 60 0 3 cU 3 0 5 LU _0 5 co D 1 0 0 5 10 15 0 5 10 15 lteration Number Iteration Number Figure 16 Cost function plotted as a function of iteration number upper left The components of the cost function are also shown The horizontal lines represent the cost contributions of the true BC s and the prior estimates of velocity and elevation noise Since there is no noise in the sampling and a perfect inverse exists the velocity contribution is driven toward zero The regularization results in abandoning the true BC s in favor of smoothing decreasing the BC cost below the truth at the expense of a small nonz
2. obs 1 1 Nobs E noise W 1 6 5 2 Velocity Velocity measurements are assumed to be sampled at a point in z y t and averaged over a vertical interval 21 22 Modeled velocity is sampled in the same way To adjust for discrepancies between model and observed bathymetries the vertical averaging is done in a bathymetry normalized vertical coordinate Velocity errors are assumed to be homogeneous and independent Ws is configured to compute the weighted mean squared velocity error F Ws d L T obs where V2 is the expected value of the i squared error the velocity noise variance and Noss is the number of velocity observations Number of ADCP points Accordingly We is diagonal 1 1 Noss Vai noise Ws U 6 5 3 Boundary Conditions W is configured to compute a weighted sum of mean sguared BC size slope and ten dency The weight factors are wo 101 and wg respectively This is written as pue nage F wen e ce ag Jn 58 Note that this measure is normalized with respect to the boundary length ds and the time span f dt It represents a weighted average of the boundary signals not a summation or integration The ideal value of wo would be 1 the expected size 2 of a boundary elevation for wy 1 expected boundary slope 2 and for wy 1 expected boundary tendency 2 It is useful to relate w to the expected size of geostrophic boundary inflows V 90
3. Mode A Elevation driven inverse error elevation error mm velocity error 10 cm s Or we 7 Or f V V N y Y Y ON 7 AS fs Wr Be SP AT de AC YN UN D 4 1 y E T gt s NR 1 1 4 t 9 9 A ra e 50 F 50 Ra Ge ae OU g n SL i J E E i x 100 Y 100 N 4 7 gt haa NON NV NN i 150F 4 150 NY ai Ni soe yA RE 09 RE KR YY As SR SI 4 Ts e _200 L LR Yh 8 200 X gt lt eri J 0 50 100 0 50 100 x km x km Figure 20 Mode A inverse solution error with elevation forcing The difference between the truth and the inverse solution is plotted Notice that the elevation error is plotted in mm compared with cm for the solution in Fig 19 38 Boundary Elevations Boundary Elevations t 72 hr t 0hr O N A O L O N A O T N 1 Elevation cm ob 4 4 SE 6 J 6 J Elevation cm t 216hr Elevation cm Elevation cm 1 21 41 61 1 21 41 61 Boundary Node Number Boundary Node Number Figure 21 Inverse solution boundary conditions at four points in time Also shown are the steady state truth BC s Initial conditions are zero for both truth and inverse solution The a b c d represent the location of the corners of the mesh as shown in Fig 9 39 39 1 2 5 R z 0 5 B 2 3 1 5 E 00000000000000000000O o
4. Ws where Ns is the number of velocity observations u and v pairs and is the identity matrix Similarily the elevation errors are assumed to be uncorrelated and equal The resulting covariance matrix W is diagonal and requires a single input the elevation noise level E noise a W xe 3 noise where N is the number of elevation observations For CASCO4 the assumption about uncorrelated errors is part of the core code It does not allow for a more general i e nondiagonal covariance structure However the assumption about the equal error levels are made at the input level and can be changed if the user wants to assign different error levels to different observations The user would be required to provide the code to bring the required data into the model and fill the covariance arrays W and W Problems in coastal ocean are always underdetermined with more unknowns than independent observations As such regularization or penalty terms W are required to control the solution The regularization used here contains terms which penalize the amplitude slope and tendency time derivative of the elevation boundary conditions 9 Error Fields Vins Expected value of the velocity error m s Epms Expected value of the elevation error m Descent Algorithm method 1 steepest descent 2 conjugate gradient Regularization Wo WI lt p gt W1 WI lt dp ds gt g PPV W2 W2 lt dp ot gt A
5. f Os Thus w scales with V g i e roughly 10 9 V 2 It is useful to scale V to the observed velocity data variance Vi adjusted for the noise level V2 bs moise V 6 5 4 Construction of W The construction of the boundary covariance proceeds as follows p s t is represented in the separate linear bases s and wj t y pirdils by t pik represents the value of the boundary condition at spatial node i and time level k The spatial basis s is identical to the FEM basis for z y along the boundary The temporal basis x t represents linear interpolation among egually spaced time points k but these points in boundary condition time are separated by a larger time step At than that used for t in the simulation This represents an extra form of temporal smoothing Expanding p in its bases we have the quadratic form for p W p dopa ask wo f r yds f vd w yi 0o ds f wedt we sds EU It is useful to recognize standard FEM matrices M the boundary integral mass matrix and IC the negative of the FEM boundary Laplacian 00 dei ds Ki Os Os Mi Lo ds 59 Analogous matrices and W comprise integrals of the temporal bases Ob Ku or dt Nu fuor dt All integrals in W are evaluated using the Trapezoidal Rule nodal quadrature With linear bases and assuming uniform At we have the following properties e M and N are diagonal K and
6. ins mesh name test6 bel file test6 3sides bel start time dd mm yyyy s 06 04 1995 0 echo file testIl echo Latitude 43 5 N stop time dd mm yyyy s 20 04 1995 0 time step 200 s THETA 0 75 TO 2 0 x 1074 s7 NNV 11 NLBS 0 kmin 0 0005 m s Nin 0 0100 m s EPSN 0 5 Table 5 Hydrodynamics model parameters for the Mode B test case The only difference from Mode A test case is the boundary element bel file The complete shell script for running the inversion is given here The script file name is casco _ModeB_basecase_both bin csh script to run casco to compute the optimal elevation bcs and then run saco to compute the final estimate of the linear solution Input files set ins_file inputs testli ins saco input file set inp_vel_error inputs inversel modeB land uxv set inp_elev_error inputs inversel modeB land uxe set samp_nodes_list inputs sample_node list set samp_dt 1800 0 time step for sampling nodes Output files set cbc_final final_bcs cbc final boundary condition time series set final_error_casco final_error uxv final error from casco set final_eerror_casco final_error uxe final error from casco set final_error_saco final _error_saco uxv final error from saco set final_eerror_saco final _eerror_saco uxe final error from saco set samp_for_root samp_saco root file name saco output set samp_adj_root samp_moody root file name moody output Casco parameters set cuto
7. 41 Summary of Finite Element Matrices 53 The include files which define the data families 64 Functions which return useful numbers 73 vi Abstract Hannah C G D R Lynch and K Smith 2003 The CASCO4b User s Guide Can Tech Rep Hydrogr Ocean Sci 226 vii F80 p CASCO4b is an inverse model whose purpose is to determine the optimal time varying elevation boundary conditions from velocity and elevation observations in the interior of the model domain This is done by minimizing the misfit between the modelled and observed velocity and elevation The model is based on a finite element implemention of the three dimensional linear barotropic eguations of motion and the adjoint method CASCO4b contains two important improvements over its predecessor CASCO3e it al lows for the inversion of both elevation and velocity whereas CASCO3e only allowed for velocity observations and 2 the coastal boundary condition no flow normal to the coast has been implemented This report is a self contained description of CASCO4b and includes a complete description of the input and output structure three examples with the complete set of input files and plots of the output a description of the theory and details of the implementation in FORTRAN code R sum Hannah C G D R Lynch and K Smith 2003 The CASCO4b User s Guide Can Tech Rep Hydrogr Ocean Sci 226 vii 80 p CASCO4b es
8. The Momentum Equation matrices are from LW91 equations 28 32 Ni 26 ho di N dd eh f 8 di ia Zo i NB 664 lean ee 0 27 0 ho di _ ow gt Hd n 28 z Wi dti Yi 2 Qj l h a Wi Note that Z performs vertical averages and Y samples the velocity basis at the bottom Y is binary in our case Table 6 Summary of Finite Element Matrices 53 6 2 Optimal Fit of Model to Data We seek the least sguares fit Minimize AU W p E r Ws EW Ja subject to the forward model constraints To enforce these we introducing the Lagrange Multipliers 4 7 A and y and minimize the functional Q Woo 6 Walde Y GI k 1 k 1 N a CLP Ve Os Vis di 01 D Tk amp Us G1 amp D Il Mit TT a M IM Mz DR Z Via qui 2 Fe UY Via Uh E Il k Ax M Ce A Gai B amp 2 C q11 D vbr 1 Ex p Ra Il Mz E Il Mz I UNT Va Via F Ve 2 6 66 Gi Sk E Il The gradient of Q provides the first order conditions for a minimum The forward model is recovered by setting the gradient with respect to the Lagrange multipliers to zero The other gradient terms are 09 Ta Mga Wes 31 00 K YM aw ia LN s E Pisa Z Bia Y Hik Pel a ral 32 an n x dg Aa IUT Vet 33 00 n
9. Y V k 0 N 1 17 qk Z V k 0 N 1 18 Errors n dy P Vi Qu Vea k 1 N 19 Ek Th Ce Ur L k 1N 20 X indicates the transpose of X and Y indicates data and definitions are as follows C V k Vectors of Elevations Velocities at time k State variables in forward model qx transport vector vertical integral of Vy vb bottom velocity sampled from V V IC s assumed known with certainty p Elevation Velocity observations in time window k 1 k 6 x Elevation Velocity errors mismatch between observation and model T U measurement functional interpolates model to space time observation p the control variables Type I BC s Ex projects p onto the simulation as elevation BC s at time k R S Known forcing BC s wind BPG and the matrices are obtained with standard finite element assembly and are summarized in Table 6 52 Let W x y be the basis for and 4 z y z be the basis for Vp and Vy Let lt gt indicate integration over x y and lt gt indicate integration over the 3 d finite element From Lynch and Werner 1991 equation 27 we have the Wave Equation matrices At bAt 10 gh VW VW 21 W W 1 5 2 W Wi 1 0 At ghVW VW 22 At bAt 10 gh VW VW 23 W W 1 gt AP r WEE SW mo Wi Wirt Oh AP KW YE aw Se 25
10. iter gt 0 gt after each time step iter 1 gt end of simulation casco_iter the casco iteration count lt 0 gt casco finished CORR I RI AR IRA RA RI AIC IC a II I I 3 2 36 ok 2 2k a kk kk kk k 2 kkk kk k kkk kk kkk kk k kkk kkk Figure 8 Definition of outputm2 22 5 Examples There are three examples based on the simple examples from Lynch and Hannah 2001 The first one reproduces the Mode experiment using velocity observations for the inversion The second one uses the same forward solution but uses elevation observations for the inversion The third example uses the Mode B experiment from Lynch and Hannah 2001 to illustrate the use of velocity and elevation observations in the same experiment and the use of the land boundary condition 5 1 Mode A with velocity forcing Consider the idealized mid latitude shelf regime illustrated in Figure 9 The length scale is of order 10 km Bathymetry is 1 D of order 100 m with uniform cross shelf slope 1 8 x 1073 Horizontal resolution is uniform at 10 km There are 10 vertical elements everywhere with sinusoidal Az tapering to 1 meter at top and bottom Temporal resolu tion is of order 3 minutes Open water boundaries surround the computational domain The vertical mixing parameters reflect linearization within a tidal regime with mean speeds of order 10 cm s at the bottom 25 cm s aloft Six moorings are deployed on a regular 60 x 40 km grid as shown each is
11. The second operation represents temporal smoothing The end point conditions pio pi and p N 1 Pin are implied Similarly the space time operation W p by symmetry the same as p W is achieved by No W k yp W s 1i 3 pU k Woe 2 1 p k 1 2p i k pli k 1 for all i k with the same endpoint conditions 7 Software Design and Implementation Notes This section describes the design of the CASCO software The purpose is provide insight into the way the software is structured so that future developers can find their way about in the code Conceptually CASCO consists of three parts a linear forward model SACO the adjoint of SACO called MOODY and the support code which uses SACO and MOODY in an iterative procedure to compute the elevation boundary conditions which minimize the unexplained errors velocities and elevations This is the picture illustrated in Fig 1 The purpose of 5ACO is take the current estimate of the elevation boundary condi tions compute the velocities elevations and the errors The purpose of MOODY is to take the errors and compute the gradient of the cost function 9 with respect to the boundary elevations the control variables The support code has four functions 1 use the gradient information 09 0p to estimate new boundary conditions for SACO 2 coordinate the iteration procedure using SACO and MOODY 3 determine whether the iteration has converged and
12. are symmetric and sparse with at most three nonzero entries per row e Nix At for interior k 2 N _1 For the endpoints k 1 and k Ny Nik At 2 e Each row of represents the negative centered second time difference ln y _ Pika 2Pik Pik 1 0 At At with homogeneous Neumann initial and terminal conditions 2e 0 Pia Pi 2 Lp i0 Pi lpia M pio Pix Pi N 1 PiN Loin lt P Pin PiN At With these simplifications it is only necessary to store two small arrays to represent Wi e W contains the time invariant spatial smoothing operator 1 1 W s Tdi fds wo MT wi K At with fdt At N 1 The dimension of W is Ni Ns since it represents time invariant spatial smoothing N is the number of boundary nodes Because it is sparse its storage requires only 3M entries 2N if symmetry is exploited Wx contains the essential time smoothing information Because of the simplicity of uniform At this matrix is also time invariant and diagonal keyed to the spatial boundary node number 1 1 We Tdt ds 60 The general space time operation p W p is achieved by Ni No We e elt k W i j pU k k 1 i j 1 4 Walii DO plik elik 1 20lik pli k 1 The first operation represents the spatial smoothing The endpoint factor cx implements trapezoidal rule integration in time e 0 5 for k 1 or N otherwise ep 1
13. as stationaryc4 elevations4 verticals5 in the way they were called before in QUODDY4 The details of the calculation of the velocities and elevations have changed from QUODDY4 The most important change is that the wave eguation and velocity calcu lations are now explicitly matrix operations This was done to simplify the construction of the adjoint model and has led to the creation of several new stationary arrays As well the wave equation has been linearized and elevation is now a prognostic variable rather than total depth and the friction terms bottom stress coefficient and vertical eddy viscosity are now constant in time As a result of these changes the user subroutines which define the vertical grid vertgrids2 the bottom friction parameters lin_stress and the vertical eddy vis cosity lin_enz are only called during the initialization phase Notice that an independent SACO can be constructed using a subset of the subrou tines shown in Fig 32 Figure 33 shows the structure of SACO which is very close the source code for the version of SACO that was used to calculate the TRUTH solutions in Section 5 There are lots of new subroutines and most of them are stored one routine per file The genealogy who calls whom is given in Fig 34 As the subroutines stabilize they can be packaged into a smaller number of files 7 2 Data Structures An attempt was made to group the data into natural families data struc
14. 1 The velocities are partial vertical averages so the portion of the water column which was sampled must be part of the data stream 2 zminobs zmaxobs and zbotobs are defined in a coordinate system which is 0 at the surface of the ocean and positive downwards The values are typically positive This is the way in which observational data is typically reported This is the opposite to the vertical coordinate used in the hydrodynamics modules SACO and MOODY 3 Observed water depths and model water depths are generally not the same There fore CASCO uses the observed water depth to scale the model depths so that CASCO samples the same fraction of the water column as was observed 4 The time is with respect to a particular year where the year is given in the header Year boundaries are not a problem Let days in the previous year be negative and days in the following year be greater than 365 or 366 Year is a leap year 5 CASCO will reject data which falls outside the mesh and the time window 6 The file format is identical to the TRUXTON input and the routines which read the file are the same An exception is that CASCO4b stops if the mesh name in the uxv file is not identical to that in the ins file 7 The data does not have to be sorted with respect to time as CASCO sorts the data internally However the output will be sorted 2 3 Elevation Errors For the pressure variable the decision was made to deal with elevation the
15. 2 2k a kk kk kk k 2 kkk kk k kkk kk kkk kk k kkk kkk Figure 7 Definition of outputs2 21 Cost 4 2 kc 4 24 I TTI TTI CC IC a II I I 3 2 36 ok 2 2k a kk kk kk k 2 kkk kk k tol ak kkk 0000 adnanrgdaaanaaanaaaanaaaanaaaeananeaaeadanaana subroutine outputm2 iter kd seckd delt nn nnv ne x y in hdown amp zeta ubar vbar Z u v zetaold DomegaDxi DomegaDrho lrvs nr amp nsteps ntbc casco_iter Output of the adjoint variables index range i i nn j 1 nnv k 1 ne l 1 nev nn number of horizontal nodes ne number of horizontal elements Kd seckd gregorian calendar day and seconds since the beginning of that calendar day iter iteration number for the current time step delt time step size in seconds x i y i horizontal cartesian nodal coordinates hdown i horizintal nodal bathymetry in i k horizontal triangular element incidence list zeta i free surface elevation ubar i vbar i vertically averaged velocity z i j nodal coordinate locations in 3 d u i j v i j nodal values of the x and y components of velocity zetaold i free surface elevation DomegaDxi gradient of cost function in full time space DomegaDrho gradient of cost function in reduced time space lrvs map from boundary node space to global node space nr number of elevation boundary nodes nsteps number of hydrodynamic model time steps ntbc number of boundary condition time steps iter 0 gt initialization purposes
16. 4 deal with the input and output SACO is a linearization of QUODDY4 The model is 3 d barotropic with user specified vertical eddy viscosity and bottom friction The only forcing terms are the elevation boundary conditions The wind and baroclinic forcing have been removed 61 7 1 Control Structure The structure of CASCO4b is shown in Fig 32 The actual source code casco4b f looks much the same The basic flow is as follows All the input information is assembled and the covariance matrices W W Ws calculated Then there is an iteration loop which uses MOODY to calculate the gradient of the cost function computes trial boundary conditions uses those boundary conditions to run SACO estimates the optimal step size adjusts the solution and decides whether convergence has occurred Finally there is a final run of SACO and the answers are written out The code is based on the QUODDY4 source code The main QUODDY module was split into 3 modules startup station_calc saco_run and the elevation and velocity cal culation modules were stripped down the linear barotropic core elevations4 verticals5 The key design decision was to write the program so that the main module of CASCO never sees the data arrays It simply calls a series of modules which do the work call these the upper modules The data is passed between the upper modules via common blocks discussed below The upper modules call the next level of working modules such
17. CASCO4b genealogy Some routines that do not call other routines are not listed 66 7 3 Important Data Items Here we discuss the important data items where they are initialized and where they are updated Observations The data structure for the observations is defined in observes h The observations are read in by getdata_moody casco4b Errors The data structure for the errors is defined in errors h The errors are initialized in initial_conditions The trial error fields UERROR VERROR EERROR are updated during the SACO run outputs2 saco_uverror saco_eerror The best estimate of the error fields uerr_best verr_best eerr_best are calculated in adjust_soln The final error is written to a file by write_error Measurement functionals The data structure for measurement functionals are de fined in weights h The velocity calculations are done by saco_basis and ver tweights The elevation calculations are done by saco_basis_elev The subroutine measure func acts as a wrapper for saco_basis and saco_basis_elev and hides the details from casco4b Boundary Conditions The trial elevation boundary conditions are passed to SACO in an array DBC see deltabcs h DBC is the responsibility of trial_bcs The best estimate of the elevation boundary conditions DBC_BEST is calculated in adjust_soln after the SACO run The final boundary conditions DBC_BEST are written to a file in write_chc casco4b Related parameter
18. I I ee 4 User Defined Subroutines 4 1 Friction I IG I I I I I a 4 2 Vertical Grid GI I I II I I I I I I I I I I I Yg 4 3 Output I II I I GI I I I I I I I I I yg 5 Examples 5 1 Mode A with velocity forcing 5 2 Mode A with elevation forcing 53 Mode B 6 Model Theory and Implementation 6 1 Discrete Equations for the Forward Model 6 2 Optimal Fit of Model to Data 6 3 Gradient Descent Algorithm 6 4 Adjoint Equations 6 5 Covariance Matrices 6 5 1 Elevation lil vi vil E 01500 N E 12 12 13 15 15 16 16 16 19 23 23 39 41 6 5 2 Velocity 6 5 3 Boundary Conditions 6 5 4 Construction of W 7 Software Design and Implementation Notes 7 1 Control Structure 7 2 Data Structures 7 3 Important Data Items TA Legacy Features a 7 5 Timing FF GI I I I I I I I I I I I I I I Y 7 5 1 Model Time 7 5 2 Boundary Condition Time 7 5 3 Measurement Time 7 6 Other I I I I II I I G I I I I I I I I
19. In the same fashion as QUODDY SACO calls a user defined output subroutine outputs2 at every time step and MOODY calls outputm2 Samples of these routines are included with the distribution of CASCO4b Brief definitions are contained in Figures 7 and 8 When constructing the output routines the user must bear in mind that the routines are called at every time step of every CASCO iteration This can generate a lot of data As well on each iteration the subroutines must either overwrite the data files from the previous iteration or create unique file names on the fly As well the SACO output 19 routine outputs2 must be able to distinguish between being inside and outside the iteration loop The variable casco_iter provides this information The output routines provided with the CASCO4b distribution provide the option of either writing a standard set of files with unigue names every CASCO iteration or writing nothing outputs2 has the ability to distinguish between calls during the CASCO iteration procedure and those during the final SACO run A useful auxillary routine used by the sample output routines is tsl_write which writes time series of the basic variables 7 7 at specified nodes 20 Cost 4 2 kc 4 24 I TTI TTI CC IC a II I I 3 2 36 ok 2 2k a kk kk kk k 2 kkk kk k tol ak kkk AMAA QAQAQA AAA AA AAAANAAAAAARAAAAQAARQAAKIGCAdAQA 4X SUBROUTINE OUTPUTS2 ITER KD SECKD DELT MESH FILINQ FILICQ amp NN NNV NE X Y IN H
20. deviation from mean sea level rather than the pressure at the bottom of the water column The convention is again that the errors are unexplained observations defined as Observation Model At initialization the Model is zero so the errors are the observations read from an external file The file format used for the elevation input and output is the uxe file UneXplained Elevation 2 3 1 The UneXplained Elevation uxe file The uxe file is a specialized version of the generic m3d format The format is 0 an arbitrary number of comment lines 1 XXXX the first reguired line 2 mesh name e g test6 3 comment 4 Year e g 1995 4 the number of columns o 6 An arbitrary number of data records Each data record a line consists of 4 items e time decimal days in UTC0 format days from beginning of the year This is yrday0 utc in the USGLOBEC Georges Bank Data Thesaurus e xobs yobs horizontal positions in same coordinate system as the FEM mesh Units meters e eobs elevation observation Units m s Notes 1 The elevation is assumed to be the deviation about mean sea level 2 The time is with respect to a particular year where the year is given in the header Year boundaries are not a problem Let days in the previous year be negative and days in the following year be greater than 365 or 366 Year is a leap year 3 CASCO4 will reject data which falls outside
21. the velocity and elevation contributions are driven toward zero The regularization results in abandoning the true BC s in favor of smoothing decreasing the BC cost below the truth at the expense of a small nonzero velocity and elevation costs 90 6 Model Theory and Implementation We solve the linearized 3 D shallow water equations aC z Y MV 0 10 OV O OV with surface and bottom boundary conditions OV NE 0 2 0 12 RV 2 h 13 The vertical viscosity and linear bottom slip coefficient N x y z y are specified exogenously These equations are rearranged as in Lynch and Werner 1991 replacing equation 10 with aC OG tro V AVE f x q sV mg 0 14 q is the vertical integral of V Initial conditions are known Open water BC s are posed in terms of an unknown distribution This is adjusted to minimize velocity errors 6 1 Discrete Equations for the Forward Model These equations are discretized on linear finite elements triangular in x y terrain following bilinear in z Semi implicit time stepping is used leading to a 3 level in time version of equation 14 for surface elevation plus a conventional 2 level momentum equation In discrete form we have Initial Condition 6 Vai Va C Vo 6 V PTS Dynamic 51 M Cea A Ce B amp 1 CT gx D vbr Ek P HR k 0 N 1 15 N Viti Va F Vk 2 G ei Ge R k 0 N 1 16 vb
22. 6 to the end of the file One line per datum Each line contains the following time step boundary node number corresponding global node number elevation Data are grouped by time step i e the first set of lines gives all boundary conditions at time 1 the next set of lines is time 2 etc For an array F i j with i boundary node number i 1 nr j time index j 1 nt ig i global node number corresponding to boundary node i there will be nt nr lines written by do j 1 nt do i 1 nr write 3 i ig i F i j enddo enddo Note that the duration of the timeseries is nt 1 dt 3 2 Cost Function Information Let Oc 9 Qe be the costs associated with the associated with the boundary conditions the velocity errors and the elevation errors respectively Then the total cost Q O O Qe and from 1 Oe P Wplp 5 O X IEA 6 and L Qe aW a 7 E Il The factor of 1 2 in the formal definition of the cost function has been ignored in CASCO4b output as it was in CASCO3e The rms velocity error Veror 18 defined from N Vi or y 070 8 k 1 and the weighted rms velocity error is o The rms elevation error Ferro 18 defined from N Error y Ex 9 k 1 13 bank150 perturbation bcs dbc_best nr ntbc 10 04 1990 3600 0 25 3600 000 116 ER 0237034E 03 8881918E 03 2865048E 03 6548486E 02 2535568E 02 5906080E 02 0271669E 03 5983905E 02 2576023E 03 251539
23. 7E 03 5073619E 02 0 NODO Qn GW KM KA KA o 1 1 1 1 1 1 1 1 1 1 Figure 3 The first few lines of a sample cbc file The mesh is bank150 with 116 boundary nodes The timeseries comprises 25 1 hour points which therefore lasts 24 hours starting at 0100 hours on April 10 1990 and the weighted rms velocity elevation is Q A standard file has not been created to hold the information about the cost function and convergence criteria and the end of each iteration To retain flexibility and make it easy to change our mind about the convergence criteria the information about the cost function and the convergence criteria is written to standard output If the standard output is saved to a file it is easy to extract the reguired information using the unix utility grep Assume that the standard output has been saved to the file log Then a list of the cost function at the end of each iteration can be retrieved by grep COST log The output has 7 columns the string COST the iteration count O QI Varor QUA Eror The first line is an information line which described the columns The second line gives the initial values Information about the costs associated with the boundary conditions the velocity mis fit and the elevation misfit can be retrieved by grep COST bcs COST vel COST elev The three numerical columns are Qe O and Qe respectively Information about progress towards the convergence crit
24. B truth Velocity timeseries at the 6 sampling locations shown in shown in Fig 17 The sampling is perfect 1 e no observational error 45 200 50 4 hours 100 Time 50 0 1 wd uoeA8 z Inversion elevation after idual unexplained Mode B res Figure 26 s uio NDOJSA 200 50 4 hours 100 Time 50 Inversion velocity after i ined idual unexpla Mode B res 27 igure F 46 Mode B Velocity and elevation driven inverse solution elevation cm uvbar 20 cm s 0 L sl 0 L T T T 1 1 t 1 I 1 1 t Ll 4 lo Cb yh fa bee a CE RE SE PAIR i bo glio i ttt 50 F 7 50 F bod 4 RF y 4 4 N t U U 4 4 t t I 4 y U t t 1 1 i 1 t t 100 T J amp 100L PF 2 h gt das EVE St ttt A fui i i c ttt 150F J 150F p de fe oh o E i 1 t 1 t 1 200 J 200 0 50 100 0 50 100 x km x Km Figure 28 Mode B inverse solution A T Mode B Velocity and elevation driven inverse error elevation error mm velocity error 10 cm s E faye NR J PENI p 50 F 4 50 A Yo A a J E E ri Z 100 100 J gt gt 150 e 150 e J liste a es E EA 200 EN ES A 4 200 tt y J 0 50 100 0 50 100 x km x km Figure 29 Mode B inverse solution
25. Canadian Technical Report of Hydrography and Ocean Sciences 226 2003 The CASCO4b User s Guide by Charles G Hannah Daniel R Lynch Keston Smith Ocean Sciences Division Maritimes Region Fisheries and Oceans Canada Bedford Institute of Oceanography P O Box 1006 Dartmouth N S Canada B2Y 4A2 Dartmouth College Hanover N H Her Majesty the Queen in Right of Canada 2003 Cat No Fs 97 18 226E ISSN 0711 6764 Correct citation for this publication Hannah C G D R Lynch and K Smith 2003 The CASCO4b User s Guide Can Tech Rep Hydrogr Ocean Sci 226 vii F80 p il Table of Contents List of Figures List of Tables Abstract R sum 1 Introduction 2 Care and Feeding of CASCO4b 2 1 Hydrodynamics Input ins 2 2 Velocity Errors 2 2 1 The UneXplained Velocities uxv file lt lt 2 3 Elevation Errors 2 3 1 The UneXplained Elevation uxe file 2 4 Inversion Parameters 2 5 Friction II IG II I I I I I I I GI II I I I I I I I Yg 2 6 Vertical Grid 2 2 IG I GI I I I I I I I I I I I I I I I ug 3 CASCO Output 3 1 Perturbation Boundary Conditions 3 1 1 The CASCO Boundary Condition file chc 3 2 Cost Function Information 3 3 Final Error ee 34 Other I I I I I I I
26. DOWN AK ENZM amp ZETA UBAR VBAR Z U V ZETAOLD UBAROLD VBAROLD casco_iter This user subroutine acts as an interface for the user to output the results at a particular time in the simulation Index range I 1 NN J 1 NNV K 1 NE L 1 NEV NN number of horizontal nodes NE number of horizontal elements NNV number of vertical nodes NEV number of vertical elements NEV NNV 1 KD SECKD Gregorian calendar day and seconds since the beginning of that calendar day ITER Iteration number for the current time step DELT Time step size in seconds X I Y I horizontal cartesian nodal coordinates HDOWN I horizintal nodal bathymetry IN I K horizontal triangular element incidence list ZETA I free surface elevation UBAR I VBAR I vertically averaged velocity Z 1 J nodal coordinate locations in 3 D U I J V I J nodal values of the X and Y components of velocity ENZM I L elemental values of the vertical eddy viscosities ZETAOLD I free surface elevation UBAROLD I VBAROLD I vertically averaged velocity MESH Name of the mesh FILINQ Name if the inq file being used to control the simulation FILICQ Name of the icq4 file used to initialize the simulation ITER 0 gt initialization purposes ITER gt 0 gt after each time step ITER 1 gt end of simulation casco_iter the casco iteration count lt 0 gt casco finished CORR I RI AR IRA RA RI AIC IC a II I I 3 2 36 ok
27. G 4 Z 0 5 0 5 0 1 0 5 10 15 20 0 5 10 15 20 30 0 5 ge _ 0 4 o o o 2 o O O lt gt 03 1 5 3 gt S 0 2 TA 2 0 mec cc 0 0 5 10 15 20 0 5 10 15 20 Iteration Number Iteration Number Figure 22 Cost function as a function of iteration The horizontal bars represent the cost contributions of the prior estimate of the velocity and elevation noise and the true BC s The velocity cost is identically zero since there are no velocity observations Since there is no noise in the sampling and a perfect inverse exists the elevation contribution is driven toward zero The regularization results in abandoning the true BC s in favor of smoothing decreasing the BC cost below the truth at the expense of a small nonzero elevation costs 40 The solution is presented without discussion We simply note that while the solution inside the convex hull of the observation locations is guite good the ability of the model to extrapolate elevation field outside the observation locations is worse than when the velocity observations were used 5 3 Mode B The final test case is Mode B from Lynch and Hannah 2001 This case demonstrates the use of both velocity and elevation observations and the land boundary condition The parameters for the hydrodynamic model are provided in Table 5 The only difference from Mode A test case is the boundary element bel file The inversion parameters are identical to those in Table 4 Hydrodynamic Parameters
28. I I I I I I I 8 Final Comments References A Source Code and Test Cases 61 62 62 67 68 69 69 69 71 13 74 75 76 List of Figures DID cwm 49 D D D LI LI D A A AM D D D N D D D LA LAD WF KR RRR RR re HE O OUR WHMrF TWO IO L A ND SO Lo XX I D cn D D A O amp O Schematic showing the major modules and the data flow in CASCO4 2 A sample ins file 5 The first few lines of a sample chc file 14 Definition of linstress 17 Definition of linenz 18 Definition of vertgrids2 19 Definition of outputs2 21 Definition of outputm2 22 Test case geometry a 25 Mode A truth the steady state 27 Truth velocity time series for Mode A 28 Residual unexplained velocity time series after inversion 28 Inverse solution for Mode A with velocity forcing 29 Inverse solution error for Mode A with velocity forcing 30 Inverse solution boundary conditions 31 Mode A Cost function as a function of iteration a 32 Elevation observations truth for Mode A 36 Residual unexplained elevations after inversion 36 Mode A inverse solution with elevation forcing lt lt 37 Mode A inverse solutio
29. Lab Report NML 99 6 November 1999 http www nml dartmouth edu Publications internal_reports NML 99 6 html Lynch D and Hannah C 2001 Inverse model for limited area hindcasts on the continental shelf J Atmos Ocean Tech 18 962 981 Lynch D Ip J Naimie C and Werner F 1996 Comprehensive coastal circulation model with application to the Gulf of Maine Contin Shelf Res 16 815 906 Lynch D and Naimie C 2002 Hindcasting the Georges Bank circulation Part II wind band inversion Contin Shelf Res 22 2191 2224 Lynch D Naimie C and Hannah C 1998 Hindcasting the Georges Bank Circula tion Part I Detiding Contin Shelf Res 18 607 639 Lynch D and Werner F 1991 Three dimensional hydrodynamics on finite elements Part II Nonlinear time stepping model nt J Num Methods Fluids 12 507 533 Press W H Flannery B Teukolsky S and Vetterling W 1986 Numerical Recipes The Art of Scientific Computing Cambridge University Press 15 A Source Code and Test Cases All of the software and data files used to make the runs described in this document are contained in the distribution of the Casco4b code This distribution can be found under Casco in the Software section of the web page of the Numerical Methods Lab at Dartmouth College try http www nml dartmouth edu circmods gom html or http www nml dartmouth edu Software The matlab scripts which make the plots are also in
30. SACO the linearized forward model MOODY the exact adjoint of SACO and the code which uses these two models in an iterative fashion to construct the optimal time varying elevation boundary conditions This is illustrated in Figure 1 The forward portion of the model SACO is a linearization of a fully nonlinear 3D circulation model QUODDY Lynch et al 1996 Its inversion is achieved by standard descent algorithms using an exact algebraic adjoint MOODY and strong dynamical constraints The cost function is a weighted least sguares blend of velocity mismatch elevation mismatch and boundary elevation size slope and tendency The control pa rameters are the open water elevation boundary conditions Solution is achieved in the time domain as a complement to the frequency domain tidal inversion TRUXTON Lynch et al 1998 The intended use of CASCO is following the removal of a best prior circulation estimate which accounts for tide local baroclinicity and the direct response to local wind forcing The remaining subtidal velocity signal is then inverted to provide far field elevation forcing The version of CASCO described herein is CASCO4b released in November 2002 CASCO4b has two improvements over its predecessor CASCO3e 1 it allows for the inversion of both elevation and velocity whereas CASCO3e only allowed for velocity observations and 2 the coastal boundary condition no flow normal to the coast has been implemented The dev
31. ard model e time element 1 ends at k 1 e time element 2 ends at k 2 71 e time element N ends at k N Assume that the model either forward or adjoint has just completed a time step e In the forward model the time element just completed is te k e In the adjoint model the time element just completed is te k 1 The mapping between k and the observation time is shown in Fig 37 1 N 1 N 2 k k X X Time weights for observation which occurred in time element k X BASISOBS 6 i fraction time step completed when the observation occurred k BASISUBS 5 i time element in which observation occurred Figure 37 The relationship between the universal model time k and the time of the observations BASISOBS is an array defined for the velocity observations observes h A similar array called EBASISOBS exists for the elevation observations The interpolation of the forward model results SACO to the observation time is done by linear interpolation For observation 7 the model velocity is Vii 1 X Va X where X BASISOBS 6 i and k BASISOBS 5 i see Fig 37 For the inverse model MOODY the contributions to the adjoint forcing at time k come from errors which occur in in time elements k and k 1 Let F be an error that occurs in time element k then dx E x Xi x be the contribution from all errors which occur in time element k and dL Ei Xi r41 be the contribution fro
32. are comment lines that describe the data lines They are defined by the QUODDY4 standard The CASCO subroutines which collect the basic information for the hydrodynamics models are derived directly from QUODDY4 1 As such the basic input file ins is based on the inq file from QUODDY4 1 The ins file consists of lines 1 17 19 20 23 24 26 28 from the QUODDY4 1 inq file The interpretation of line 21 has changed 1 Comment 2 testH forward model truth through flow 3 day rampup in time 3 Case name L 4 test6 5 Boundary element incidence list 6 test6 4sides bel 7 Initial condition file 8 COLD START 06 04 1995 0 0 9 Echo file 10 testH1 echo 11 Simulation parameters 10 20 04 1995 0 0 12 SI UNITS units 13 1 00 1 00 x y and z scaling factor 14 43 500 degree latitude 15 10 0 minimum depth 16 15 04 1995 0 0 end date d m y and time sec of simulation 17 200 time step seconds 18 0 75 THETA 19 0 0002 TAUO 20 11 number of vertical nodes 21 0 nonlinear bottom stress factor NLBS 22 0 0100 minimum vertical viscosity diffusivity 23 0 0005 minimum bottom stress coefficient 24 0 50 time weighting factor for vertical stress diffusion Figure 2 A sample ins file The column of line numbers was added for presentation and is not part of the file format 2 2 Velocity Errors The CASCO convention is that errors are unexplained observations defi
33. artures At the northeast corner flow is incorrectly entering the domain and at the northwest corner too much flow is entering the domain More serious is the downstream boundary where a nearly impermeable wall has been artificially created the inverse solution is totally bogus there Figure 14 shows the discrepancy with truth Accuracy is high over much of the domain including the entrance region and well below the anticipated noise level in the interior The zones of obvious error are well above the noise level Generally interpolation accu racy is high as is extrapolation to portions of the boundaries But the balance of the boundary is grossly inaccurate These are patterns of mean error or bias Lynch and Hannah 2001 Examination of the true and inverse elevation boundary conditions is instructive Figure 15 displays both at four points in time Compared with truth the inverse solution has a tendency to smear the corners and flatten the outflow area The corner smearing lessens the mean sguare slope while still producing the same aggregate effects at the moorings The downstream artificial wall has little dynamical interaction with the flow at the moorings so the preferred solution is flat there to satisfy the regularization The cost function is plotted in Figure 16 It is clear that this solution has reduced the BC cost below that of the true BC s via the smoothing perturbations described above at little cost to the veloci
34. ase See Lynch and Hannah 2001 for a discus sion of accuracy and precision as it relates to inverse models The residual unexplained velocity data is reduced to a fraction of a cm s Figure 12 Note that the inertial oscilla tions are left in the residual a consequence of the heavy temporal smoothing At 72 23 Hydrodynamic Parameters ins mesh name test6 bel file test6 4sides bel start time dd mm yyyy s 06 04 1995 0 0 echo file testHl echo Latitude 43 5 N stop time dd mm yyyy s 20 04 1995 0 0 time step 200 s THETA 0 75 TO 2 0 x 1074 s7 NNV 11 NLBS 0 kmin 0 0005 m s Nin 0 0100 m s EPSN 0 5 Hidden Hydrodynamic Parameters Aza lm CLOSURE CONSTANT Inversion Parameters 0 001 2 0 01 Tmax 100 Wo 0 0 Wi 1012 Wa 0 0 Vicise 0 03 m s E noise 0 01 m At 72 hours Table 4 Parameters for the Mode A test case 24 of m ot d J AOS i SSSI 50 100 150 T DS NR 50 PES 50 RSS SSL NL O 100 RS 100 e e AUS SN 2150F DSS NN 150 J x N KO o o 200 L 200 L aj 0 50 50 100 Figure 9 Test case geometry Left 231 node triangular mesh Right bathymetry dash contours and mooring locations dots Corners are indicated for later reference 25 hr Figure 13 shows the complete solution at the end of the simulation Over much of the domain the solution is accurate there are notable dep
35. b module Fig 32 These functions are listed in Table 8 and are collected together in the file numbers f The details of the gradient descent algorithm was changed for CASCO4b In CASCO3b the magnitude of the trial boundary conditions was scaled with the amplitude of the cost function smaller steps with smaller cost In CASCO4b the amplitude of the trial bound ary conditions is constant The scaling is as follows The trial boundary conditions are scaled so that the largest element is 1 Then the trial boundary conditions are multiplied by the parameter scale DoDr In the distribution scale DoDr 0 1 which means that the largest element of the trial boundary conditions is 10 cm This parameter is set in the file casco params h We think that what is really desired is to have the trial boundary conditions generate elevations and velocities in the SACO run that are of the same order as the misfits This will reduce the likelihood of roundoff errors dominating the descent algorithm A determined attempt was made to use explicit typing of all the variables This allows one to find typing errors by compiling with the u option which finds all the untyped 13 variables 8 Final Comments This document provides the user with complete description of the input and output structure of CASCO4b and three examples with the complete set of input files output files and plots of the output These should help new users to get started and help to ver
36. bc call write_error call dbc_best_to_dbc call saco_run 1 call error_to_best call write_error write the perturbation bcs write the final errors version 1 IMove the bcs around IDo the final saco run IMove errors to arrays to write from Write the final errors version 2 Figure 32 The structure of CASCO4b 63 Table 7 The include files which define the data families CASCO DIM casco_params h constants h covariance h deltabcs h errors h friction h gradcost h mesh h misc h numbers h observes h saco h sacotime h startstop h statarrays h statistics h weights h wrho h h Table 7 The name of the common block is usually the same as the name of the All the array size specifications The casco specific input parameters All the physical constants and time stepping parameters The variance and covariance of the elevation and velocity errors The perturbation boundary conditions calculated by MOODY The velocity and elevation errors The bottom friction and vertical eddy viscosity information The gradient of the cost function All the input related to mesh geometry and boundary conditions Miscellaneous stuff Definition of functions which return useful numbers not used yet The time location and values of the unexplained observations which are the forcing for casco The prognostic variables for 5ACO Timing information for the forward run and the mapping from the model time points to the
37. both script to run ModeB with elevation and velocity forcing casco_ModeB_basecase_justelev script to run ModeB with elevation forcing casco_ModeB_basecase_justvel script to run ModeB with velocity forcing cleanup script to delete output files do_the_manual m runs the test cases and plots the output Readme basecase information saco_ModeA_basecase script to run modeA TRUTH saco_ModeB_basecase script to run modeB TRUTH sample_output the output from the runs shown in the manual test6 3sides bel bel file for modeB test6 4sides bel bel file for modeA test6 bat mesh files required for model runs test6 bnd test6 ele test6 lev test6 nod run sample_output Sample output from three test cases presented in the manual run sample_output modeA_truth Output from modeA truth run run sample_output modeB_truth Output from modeB truth run run sample_output modeA_justelev Output from modeA elevation inversion run sample_output modeA_justvel Output from modeA velocity inversion run sample_output modeB_both Output from modeB elevation and velocity inversion src source code for casco4b and saco3 adjust_soln f bcs_saco_moody f 78 BuildWrho3 f cascoi_fixsubs f casco2_sprspak f casco4a f CASCO DIM casco info casco_info f casco_input f casco_params h constants h convergence2 f covariance h dbc_best_to_dbc f deltabcs h domega_write f elevationm4 f elevations4 f EMATpreMULT f errors h error_to_best
38. boundary condition time points The start and stop times for the model runs The stationary arrays calculated by stationaryq4 The current values of the basic statistical measures The interpolation functions for the observations which define the measurement functionals The weight function for the elevation boundary conditions include file but not always The include files contain the variable declarations common blocks and documen tation for the data structures The documentation includes assigning responsibility for filling the data structures All the data is hidden from the main part of CASCO The common blocks are used to pass data between the upper modules defined as those that appear in Figure 32 There is one family per include file so that modules only see data that they need The upper modules communicate with the next level by a combination of common blocks and subroutine arguments 64 SACO Main Module startup IRead the mesh files and such station calc ICalculate all the constant and common arrays 1 get friction iter IRead the friction values set_times IPrecalculate all the timing information getdata moody Get the observations measure_func ICompute the measurement functionals 1 saco_run iter IRun the forward model Figure 33 The structure of 5ACO based on CASCO modules To sucessfully run SACO using this structure the user will need to modify the subroutines that provide the elevation bou
39. c 0 4 al Velocity cm s 0 50 100 150 200 Time hours Figure 12 Residual unexplained velocity time series after inversion of Mode A using velocity observations 28 uvbar y driven inverse solution o o Vv gt lt a BS Ss w gt o wy A RUN CE O A A eke o NA TE ANT to hay th AU Ob ND Pots tu hse V L FIN E F a IHR A o o o LO o LO o l Kt HIS ST wy A nes 4 9 T L 2 G o o gt o LO o LO o i H n 9 100 50 x km 100 50 x km Figure 13 Inverse solution for Mode A with velocity forcing 20 Mode A Velocity driven inverse error elevation error mm velocity error 10 cm s E S ET an Co ll NA NS a Gae WD mn te A La 50 F 4 _50 L o deo e 4 4 Ho owe La NS se oan oV E 100 d amp zool gt gt N Ni NY 150F 7 150F RE 4 NV NE A NA Sd s ad A A MS VV NA _200 L 4 200 Ei T Td 0 50 100 0 50 100 x km x km Figure 14 Inverse solution error for Mode A with velocity forcing The difference between the Truth solution Fig 10 and the inverse solution Fig 13 is plotted Notice that the elevation error is plotted in mm whereas the solution is in cm 30 Boundary Elevations Boundary Elevations
40. capable of sampling and vertically averaging the entire water column with error The observations at these locations will be inverted to obtain the pressure signal on the boundary A forward model run of SACO forced by boundary conditions constitutes truth The most basic mode is along shelf throughflow mode A of Lynch and Hannah 2001 The 5ACO truth is shown in Figure 10 This solution was started from rest with BC s ramped up linearly over 72 hours and held constant thereafter The slow spinup is effectively quasi steady as illustrated in Figure 11 The effects of topography and friction are visible as horizontal gradients but they are not dominant and small damped inertial oscillations are similarly present Effectively we have a nearly 2D x z mode The parameters for 5ACO are given in the upper part of Table 4 For the inversion the velocity noise is presumed totally uncorrelated with RMS size Vioise 3 cm s The expected boundary condition slope is 1079 compatible with the observed 10 cm s velocity The time scale for boundary conditions At 72 hours and the overall length of simulation is 216 hours Penalties for boundary condition size and tendency are zero the latter consistent with the prior smoothing by the 72 hour time step These and all other case parameters are summarized in Table 4 For this example no noise is added to the data The resulting inversion is a measure of the accuracy of the inversion for this c
41. cluded Some of the scripts make use of the OPNML matlab tool box written by Brian Blanton which can be found at http www opnml unc edu The version used are the Matlab5 versions Below is the text of the file readme contents included with the November 2002 distribution of Casco4b Hadad k kk k kk k kk k k TI TALI TI ada LI TI kk RO ROF k LI Distribution of Casco4b November 2002 Source Code casco An inverse model for estimating the optimal time varying barotropic boundary conditions from velocity errors saco An independent version of the linear forward model embedded in casco The only forcing is the pressure boundary conditions This code was used to create the TRUTH for the test cases modeA and modeB Directories bcs_forward boundary conditions for the test problems inputs forcing data and input mesh mesh files for the test problems mfiles matlab file for plotting the casco output run the home of the test cases src casco4b the casco4b source code src casco4b saco3 source code for saco3 uses code from casco Useful Information Files bcs_forward readme bcs describes the elevation bc forcing files inputs Readme inputs describes the contents of the inputs directory readme contents this file run Readme basecase describes the scripts for running the test cases FO E k k kk k k kk k kk k k kk k kk GK LI CI LS Contents of the directories bcs_forward boundary con
42. ction between the two methods is done at run time see Table 2 Convergence of the CASCO4b iteration is achieved when any of the following condi tions is reached e The objective function is small relative to its initial value Q lt a Og 10 e The most recent change in the objective function is small relative to its current value AQ lt 690 e The iteration count exceeds a maximum Imaz More sophisticated convergence rules are reguired One possibility is to consider the convergence of the cost due to the boundary conditions Notice however that the cost due to the boundary conditions is not reguired to decrease monitonically only the total cost 9 is reguired to do that 2 5 Friction There are two friction parameters in CASCO the linearized bottom friction parame ter x and the vertical eddy viscosity N The bottom friction parameter can vary horizontally z y and the vertical eddy viscosity can vary in all three dimensions N v y z In the present implementation they are both independent of time In general z y and N z y z are taken from a run of the full nonlinear model and are brought into CASCO via the user subroutines lin_stress and lin_enz However they can also be set to the minimum values specified in the ins file This is controlled as follows If the parameter NLBS 0 line 21 in ins then k Kmin Where Kmin is defined on line 23 of the ins file If NLBS 1 then lin_stress is call
43. ditions for the test problems 76 readme bcs test6 testH atw ak0 0005 resbcs 3sides s2r test6 testH atw ak0 0005 resbcs s2r test6 testl atw ak0 0005 resbcs 3sides s2r test6 testl atw ak0 0005 resbcs s2r inputs observational files for test problems from the casco4b manual inversel adcp_out n0 0 uxe inversel adcp_out n0 0 uxv inversel adcp_out n0 uxe inversel adcp_out n0 uxv inversel modeB land uxe inversel modeB land uxv Readme inputs sample_node list testHi ins testIl ins mesh Mesh geometry and bc type specification files test6 3sides bel test6 4sides bel test6 bat test6 bnd test6 ele test6 lev test6 nod testgeom m testgeom ps mfiles Matlab files for plotting casco output The plotting routines should be called from within the run directory edges2chains m plot_cbc_s2r m plot_contourmat m plot_cost_function2 m plot_cost_function m plot_final_elevation_error m plot_final_velocity_error m plot_input_elevation_data m plot_input_velocity_data m plot_inverse_error m TT plot_inverse_solution m plot_truth m read_cbc m read_m3d m read_uxe m scattered_contour m suptitle m run The working directory from which the test cases can be run casco_ModeA_basecase_justelev script to run ModeA with elevation forcing casco_ModeA_basecase_justvel script to run ModeA with velocity forcing casco_ModeA_basecase_both script to run ModeA with elevation and velocity forcing casco_ModeB_basecase_
44. e inversion process such as the expected errors the regularization terms and the convergence criteria Section 2 4 5 The bottom friction and vertical eddy viscosity Section 2 5 6 The specification of the vertical grid Section 2 6 CASCO4b can be run with elevation data and velocity data together or either data type separately Table 1 shows the input stream required by CASCO4b ignoring inputs required by the output routines which come after The symbols are defined in the rest of this section The examples in Section 5 will illustrate the complete data stream including data required by the sample output routines 2 1 Hydrodynamics Input ins The basic input to the hydrodynamics models SACO and MOODY is the ins file file ins in Table 1 a sample of which is shown in Figure 2 In short e Line 2 is a comment line e Line 4 defines the mesh name The required files are test6 nod test6 ele test6 bat test6 lev where test6 is the mesh name e Line 6 defines the name of boundary element file e Line 8 defines the start time of the simulation The start time is read using a for matted read statement that allows for a single space after COLD START and then the date in DD MM YYYY seconds format The FORTRAN read statement has the following format i2 1x i2 1x i4 1x f21 12 Note that all CASCO runs require a COLD START e Line 10 defines the echo file where basic diagnostic information about the hydro dynamics
45. e one of two methods Method of Steepest Descent or Conjugate Gradient Method Press et al 1986 Golub and VanLoan 1983 In the steepest descent method p is simply the negative of the gradient of the cost function p 009 0p In the CGM we identify a sequence of gradients g 00 0p and directions hy with indicting iteration number The direction is computed as a blend of the current gradient and the previous direction hia gi yhi 47 _ Im 9 gi 48 gog p is then computed as an increment of arbitrary length in this direction Pisa Ria giga yh 49 The method reduces to the steepest descent method if y is arbitrarily set to zero Initially ho go is sufficient to get things started The gradient descent algorithm is summarized as follows e Given p 6 and the gradient Di compute the direction of descent p A single forward model calculation with 0 1 is sufficient to calculate 6 and e e Compute the optimal value of 0 e By superposition update p and c This completes the new improved forward solution It is followed by an adjoint solution forced by the new values of 6 and e resulting in a new gradient oH and the cycle continued to convergence suitable convergence rule is the achievement of vanishingly small values of the gradient thereby satisfying the final first order condition for an optimal solution Note that only a single forward and adjoint solution is regui
46. ed to fill z y The calling subroutine get_friction applies the filter max k Kmin after calling lin_stress Values of NLBS other than 0 and 1 cause CASCO4 to stop The vertical eddy viscosity is handled in a similar fashion The parameter CLOSURE is specified as a parameter in the file friction h If CLOSURE CONSTANT then N 2 y 2 Nain Where Nmin is determined by line 22 in the ins file If CLOSURE SPATIAL then N x y 2 is specified by the user subroutine lin enz The calling subroutine get_friction applies the filter N max N Ninin after calling lin enz Values of CLOSURE other than CONSTANT and SPATIAL cause CASCO4 to stop These two user subroutines are discussed further in Section 4 The input data stream for lin_stress and lin_enz have not been dealt with It is assumed that the users will modify the CASCO input stream to meet the need 2 6 Vertical Grid The vertical grid is defined by the user subroutine vertgrids2 Section 4 The number of vertical nodes NNV is defined in the ins file and is passed to vertgrids2 If CASCO is being used in concert with a nonlinear circulation model to construct a nonlinear inverse then perhaps the vertical grid should be taken from the nonlinear circulation model In development it has been simpler to use the standard subroutine sinegrid with ll its controlling parameter dzbl the distance between the lower or upper two nodes specified in vert
47. efinitions of the the array sizes and needs to be modified by the user The user subroutines are defined further in the rest of this section Some auxillary output subroutines which are part of the distribution are also listed 4 1 Friction As described above the user needs to specify two friction parameters the linearized bottom friction coefficient x and the vertical eddy viscosity N These can be spatial constants or defined by the user subroutines lin_stress and lin_enz If NLBS 1 then the subroutine lin_stress is called to specify k x y If CLOSURE SPATIAL then the subroutine lin_enz is called to specify N a y z Figures 4 and 5 contain the definitions of lin_stress and lin_enz After calling lin_stress and lin enz CASCO3 applies the minimums given in the ins file 4 2 Vertical Grid The distribution of the vertical nodes is determined by the user subroutine vertgrids2 The subroutine is called once during the initialization phase and the number of vertical nodes NNV as defined in the ins file and is passed vertgrids2 is defined in Figure 6 If CASCO is being used in concert with a nonlinear circulation model to construct 16 COAG nn k kk kk k EE 2k kk kk 2k kkk ok 2k C Read linearized bottom friction coefficients from data files SUBROUTINE lin_stress kd seckd ncall ak akmin nn INCLUDE CASCO DIM integer ncall integer kd real seckd real ak nndim real akmin integer nn Index range I 1 NN
48. elopment and applications of CASCO3e are described by Hannah and Lynch 1999 Lynch and Hannah 2001 and Lynch and Naimie 2002 This document is intended to be a self contained description of CASCO4b Sections 2 5 provide a complete description of the input and output structure of CA5CO4b and three examples with the complete set of input files and plots of the output These should help new users to get started and help to verify that the code is working in a new environment Sections 6 and 7 contain descriptions of the theory and implementation that should help the advanced user to modify the code Appendix A contains a list of the files that constitute the distribution of CASCO4b CASCO4 is written in FORTRAN77 The code attempts to be portable So far it has been compiled and run successfully on the following operating systems SGI Linux on an Intel processor and Sun Solaris The development of CASCO and the distribution discussed herein assumes the existence of a unix like operating system for features such as long file names shell scripts and tools such as awk grep and make As well this document assumes that the user has some experience with finite element circulation Turbulence Unexplained Estimate Data Output V C Error Boundary Condition Supplement Forcing Cov Errors BC s BC s Basis Convergence Criterion Figure 1 Schematic showing the major modules a
49. eria can be obtained by grep 14 CONV log There are 5 lines for each iteration The first line labelled CONVER1 shows progress toward the e criterion Table 2 the second column is the iteration count the third 2 00 and the fourth a The second line labelled CONVER2 shows progress towards the ez criterion the third column is AQ Q and the fourth is The next three lines are information about the relative change in Qs Qu and 9 which are not yet used The sixth line reports whether convergence has been achieved Information about whether the iteration count has exceeded the maximum value is not reported explicitly However the iteration count is reported 3 3 Final Error The final error that portion of the observations which cannot be explained by the per turbation boundary conditions is written to files specified by the user This action is independent of the user subroutines The file formats are the uxe and uxv formats 3 4 Other The standard output from CASCO should be saved to a file the log file while CASCO is running The output contains the convergence and cost function information see above and lots of other diagnostic output Both the log file and the echo file Section 2 1 should be checked when CASCO crashes The user subroutines distributed with CASCO4b can write output time series of elevation and depth integrated velocity for both the forward and inverse model runs This can generate a lot of outp
50. ero velocity cost Since there are no elevation observations the elevation cost is identically zero 32 The complete shell script for running the inversion above is given here The script file name is casco ModeA_basecase_justvel November 2002 distribution bin csh script to run casco to compute the optimal elevation bcs and then run saco to compute the final estimate of the linear solution Input files set ins_file inputs testHi ins saco input file set inp_vel_error inputs inversel adcp_out n0 uxv set inp_elev_error none set samp_nodes_list inputs sample_node list set samp_dt 1800 0 time step for sampling nodes Output files set cbc_final final_bcs cbc final boundary condition time series set final_error_casco final_error uxv final error from casco set final_eerror_casco final_error uxe final error from casco set final_error_saco final_error_saco uxv final error from saco set final_eerror_saco final_eerror_saco uxe final error from saco set samp_for_root samp_saco root file name saco time series set samp_adj_root samp_moody root file name moody time series CASCO parameters set cutoff1 0 001 cutoff on cost function set cutoff2 0 01 cutoff on change in cost function 30 maximum number of iterations set itermax set w0 0 0 weight for elevations set wi 1 0E12 weight for elevation slopes set w2 0 0 time weight set bcs_tscale 72 0 smoothing t
51. error The difference between the Truth and the inverse solution is plotted Notice that the elevation error is plotted mm compared with cm for the solution in Fig 28 48 Elevation cm Elevation cm 4 Boundary Elevations t 0hr 21 41 61 Boundary Node Number Elevation cm Elevation cm 2 4 Boundary Elevations t 72 hr oo 21 41 Boundary Node Number 61 Figure 30 Inverse solution boundary conditions dots at four points in time Also shown are the truth BC s line Initial conditions are zero for both truth and inverse solution The a b c d represent the location of the corners of the mesh as shown in Fig 9 Elevation values are not shown for the part of the boundary that is land a to b 49 3 2 5 2 5 2 r D a 815 o gt Toi 1 5 o o 1 E 1 o gt p de 0 0 S 6 6 606 5 0 5 10 15 0 5 10 15 1 0 25 _ 0 8 _ 0 2 0 6 0 15 O bs 04 W ui A 0 24 0 05 0 0 0 5 10 15 0 5 10 15 Iteration Number Iteration Number Figure 31 Cost function as a function of iteration The horizontal bars represent the cost contributions of the the prior estimate of the velocity and elevation noise and the true BC s The expected elevation cost is 1 and is hidden in the frame of the plot Since there is no noise in the sampling and a perfect inverse exists
52. f FEMPAKA f friction h get_bcstime f getdata_moody4 f getelevobs f get_errvar f get_friction1 f get_nsteps f getobs_4c f getobs_8c f get_sacotime f getvelobs f get_Wrho f gradcost h indexx f initial_conditions f inverse_parameters h lin frict f main tex makefile measure_func f mesh h misc f misc h moody_run f 9 newnorms f numbers f numbers h observes h outputm2 f outputs2 f reorderinbe f rms f saco_basis_elev f saco_basis f saco h saco_run3 f sacotime h set_times f startstop h start_subs f startup f statarrays h stationaryc3 f station_calc2 f statistics f statistics h stepsizeb f sup_norm f time_basis f time_manip f trial_bcs_cgm f trial_bcs f vertgrids2 f verticalm4 f verticals5 f weights h wrho h write_cbc f write_eerror f write_error f write_uverror f src saco3 source code for saco3 uses code from casco bcs2 f makefile saco3 f 80
53. ffi 0 001 cutoff on cost function set cutoff2 0 01 cutoff on change in cost function set itermax 100 maximum number of iterations set itermax 30 maximum number of iterations set w0 0 0 weight for elevations set wi 1 0E12 weight for elevation slopes set w2 0 0 time weight set bcs_tscale 72 0 smoothing time scale for bcs hours set exp_rms_vel_err 0 03 expected rms velocity error m s set set set exp_rms_elev_err 0 01 expected rms elevation error m method_descent 2 casco_exe src casco4b casco_exe lt lt EOD cutoff1 cutoff2 itermax wO w1 w2 exp_rms_vel_err exp_rms_elev_err bcs_tscale method_ descent ins_file inp_vel_error inp_elev_error samp_ samp_ samp_ samp_ samp_ nodes_list dt adj_root nodes_list dt 4 samp_for_root cbc_ final final_error_casco final_eerror_casco samp_nodes_list samp_dt samp_for_root final_error_saco final_eerror_saco EOD The truth solution Fig 23 and the standard figures are provided without discussion See Lynch and Hannah 2001 for a discussion of the velocity only inversion of mode B The meeting of the land with the downstream open boundary is not as smooth as one would wish Downstream is defined relative to shelf wave propagation It is hard to see but in the lower left corner of Fig 28 and Fig 29 there are larger than desired velocity vectors parallel to the coast We believe that the followin
54. g explanation holds The downstream boundary does not affect the solution at the observation locations therefore it is entirely determined by the regularization In the solutions with four wet boundaries the regularization provides some smoothness from one boundary segment to the next However the introduction of the land boundary on the left weakens the smoothness constraints in the lower left corner This results in larger elevation gradients and large velocities in that corner Clearly a better downstream boundary condition is reguired 43 y Km Mode B Truth elevation cm uvbar 20 cm s Or 1 or Tp de Litta fe cle a rt t tof ppi rt to tf I v t ttt 1 I i 1 t t t 50 7 50 bid 0 needs K TERE Et DI sei ga St Oh oh da C 1t t ft 1 1 1 i t T E tt 100 F 7 lt 100 o HA pe Y Wed mU gt na CRE RE i i t t t t t t Ry i t t t 150 F 4 150 dan A RAR dort 1 t t t t t 1 t t i n ti t 1 200F 200 x AA 0 50 100 0 50 100 x km x km Figure 23 Mode B truth 44 Elevation cm O DST 0 50 100 150 200 Time hours Figure 24 Mode B truth Elevation timeseries at the 6 sampling locations shown in shown in Fig 17 The sampling is perfect i e no observational error o T gx O lt c al ym Velocity cm s al 0 50 100 150 200 Time hours Figure 25 Mode
55. grids2 Thus the vertical grid information does not appear in the data stream 3 CASCO Output The outputs from CASCO fall into four primary categories 1 The answer the perturbation boundary conditions which minimize the errors Sec tion 3 1 2 Information about the cost function and the convergence of the solution Sec tion 3 2 3 The final errors at the observation locations Section 3 3 4 Information about the forward and adjoint solutions Section 3 4 3 1 Perturbation Boundary Conditions When CASCO4b is finished the optimal boundary conditions are written to a file spec ified by the user The format of the file is part of the CASCO specification and the file is written out without regard to the user subroutines 3 1 1 The CASCO Boundary Condition file cbc The file for CASCO boundary condition output cbc contains time series for the ele vation at the boundary nodes of the mesh along with information on the start time the time step for the timeseries and the total number of time points in the file The unit of elevation is metres A sample cbc file is shown in Fig 3 line I CASENAME line 2 comment line line 3 Start time for the timeseries in Human Units day month year seconds 3 integers one real all separated by one blank Seconds is the elapsed time since the day began line 4 number of time points boundary condition time step At line 5 number of boundary nodes 12 line
56. her string Just beware of spaces The string cutoffl says put the value of the variable cutoff1 here So the line cutoff1 cutoff2 is replaced by 0 001 0 01 before being given to CASCO4b The first eight lines are the CASCO4b input stream with the order as given in Table 1 The next six lines supply the time series sampling information for the output subroutines Clearly this will change as the output subroutines develop The time series output can be suppressed by setting the root file name samp_for_root and or samp_adj_root to the string none This reguires that the next two lines in the input stream be deleted See the output subroutines for more details The next two lines are the file names for the perturbation boundary conditions and the final error The next three lines feed the time series sampling information to the output subrou tines for the final SACO run after the iteration process is complete The next to last line is the file name for the final error from the final SACO run The final line EOD closes the input stream 5 2 Mode A with elevation forcing The second test case is the same as the first except that elevation rather than pressure is used for the observations The shell script for running this test case is practically identical to the previous one The script file name is casco ModeA_basecase_justelev The only differences are in the input file section where the elevation input file has a name and the ve
57. ify that the code is working in a new environment This document also provided a description of the theory and implementation that should help the advanced user to modify the code CASCO4b is a straight forward implementation of the equations presented in Sec tion 6 We believe that the coding of CASCO4b is correct However there may be bugs hiding in corners that we have not explored Very few extras have been provided to make the code easy to use The directions for future enhancements await some practical experience with the model in real applications There are many areas where CA5CO could be improved These include e Convergence rules Section 2 4 e Details of the estimate of the trial boundary conditions Section 7 6 Improved covariance models Section 2 4 e Improved boundary conditions in regions that do not influence the observations Section 5 Allowing for nondiagonal covariance structure next paragraph There is one area for future development that will reguire both theoretical develop ments and practical experience In their one experiment with a short boundary condition time step 3 hr Lynch and Hannah 2001 found that inertial oscillations developed in both the forward and inverse models They were difficult to control and were only re moved after several hundred iterations The CASCO3e User Guide Hannah and Lynch 1999 also provides an example of this problem in the example of the use of the weight funct
58. ime scale for bcs hours set exp_rms_vel_err 0 03 expected rms velocity error m s set exp_rms_elev_err 0 01 expected rms elevation error m set method_descent 2 1 gt steepest descent 2 gt conjugate gradient The executable program set casco_exe src casco4b casco_exe lt lt EOD cutoff1 cutoff2 itermax wO w1 w2 exp_rms_vel_err exp_rms_elev_err 33 bcs_tscale method_ descent ins_file inp_vel_error inp_elev_error samp_nodes_list samp_dt samp_adj_root samp_nodes_list samp_dt samp_for_root cbc_final final_error_casco final_eerror_casco samp_nodes_list samp_dt samp_for_root final_error_saco final_eerror_saco EOD The first half of the script defines the variables which are used to feed CASCO4b The second half runs CASCO4b and feeds it the required input The first line specifies that the script is a c shell script this defines the syntax used to define the variables Note that defines a comment The three lines after the comment Input files define the input files the ins file the uxv file the uxe file and a sample node list In this case there are no elevation obser vations and the dummy file name none is used to indicate this CASCO4b understands none and NONE The sample node list is required by the user subroutines outputs2 and outputm2 it defines the nodes where time series of elevation and depth averaged velocity will be reported The output is writ
59. ions to provide temporal smoothing This example Section 5 2 of Hannah and Lynch 1999 was not provided in this document Lynch and Naimie 2002 found that a practical way around this problem for their application was to write the boundary condition time series as a convolution of the wind stress Another possibility is to run CASCO with a coarse boundary condition time step as a first cut and then refine the time step 74 We believe that the heart of the problem is that the misfits or errors are presented to the inverse model as a series of impulses at the observation time and these impulses excite inertial oscillations One way forward would be smear the misfits in time The most general way to implement this would be to allow for non diagonal covariance structure for the misfits which would reguire revisions to the code of the inverse model There may be other useful methods Acknowledgements We thank Chris Naimie for his efforts on the software engineering of OUODDY4 His work made the development of CASCO much easier than could have been We also thank Chris for the measurement functional code This work was supported in part by the U S National Science Foundation NSF the U S National Oceanic and Atmospheric Administration NOAA References Golub G and VanLoan C 1983 Matrix Computations The Johns Hopkins University Press Hannah C and Lynch D 1999 Casco3 User Guide Dartmouth Col lege Numerical Methods
60. is written e Lines 13 defines some scale factors which are best set to 1 0 e Lines 14 to 17 define the latitude the minimum depth the simulation end time and the time step for the hydrodynamics mode The simulation end time is read with an unformatted read statement therefore Line 16 reguires 3 integers day month and year and a floating point number seconds The time in the example uses the same format as reguired by Line 8 e Line 18 defines the parameter THETA which controls the timing stepping algorithm in the wave equation THETA 0 75 is a good value e Line 19 defines another wave equation parameter 70 and the reader is referred to the literature 7 0 0002 is standard e Line 20 defines the number of vertical nodes e Line 21 defines the parameter NLBS which controls how bottom friction is han dled If NLBS 0 then the linear bottom friction parameter is constant and taken as the value specified in Line 23 Otherwise the user subroutine lin_stress is called See Sections 2 5 and 4 e Linc 22 defines the minimum vertical eddy viscosity e Line 23 defines the minimum value of the linear bottom friction coefficient e Line 24 defines the time weighting factor for the solution of the vertical structure of the velocity 0 5 is a good value In CASCO4b this variable is not used The time weighting factor is hardcoded to be 0 5 See Section 7 4 for more details The lines with text inside the such as Comment
61. locity input file is now none The input file section of the script is shown here Input files set ins_file inputs testH1 ins saco input file set inp_vel_error none set inp_elev_error inputs inversel adcp_out n0 uxe set samp_nodes_list inputs sample_node list set samp_dt 1800 0 time step for sampling nodes 39 Elevation cm NO sia NANO LS o 50 100 150 200 Time hours Figure 17 Elevation observations truth for Mode A at the 6 sampling locations Fig 9 The sampling is perfect i e no observational error 0 4 EN minn Tmn ia Elevation cm Q D T I o LS o 50 100 150 200 Time hours Figure 18 Residual unexplained elevations after inversion of Mode A with elevation forcing 36 y Km Mode A Elevation driven inverse solution elevation cm uvbar T T T T T T 20 ol ol E DUE ea bh E DA Bt Uu U 1 Kr 4 ok E We Be A he GE ae oa 50 L J 50 F O e A Y Ek K sth RO 4 ch a e i BP LE ae dI E U U 4 1 4 4 y 100F J lt 100 L KLT ha E L E LDE 2 EES e heed p D POOR CEL UE OA BO Mob Sie lh ae of be Abs Ne yb 150 7 150 Sal op Ye E ka x 9 U L 4 L y X NU Am Gh WM ye BRE a A DS sai AS 200F J 200F _ L L L L L L 0 50 100 0 50 100 x km x km Figure 19 Mode A inverse solution with elevation forcing 37
62. ltre time step for the boundary conditions hours Convergence Convergence when Q lt 100 9 Convergence when AQ lt 60 Imaz maximum number of iterations Table 2 Input parameters related to the inversion process in CASCO4b There is a second control on the time domain obtained by specifying a time interval for the boundary conditions Atgc which is greater than the time step in the hydrodynamic model This control mechanism has two effects it limits the time variability and it reduces the size of the solution space The construction of W requires the parameters W0 W1 W2 which penalize am plitude slope and tendency respectively The scaling for these parameters are given in Table 2 where p is the boundary elevation 0p 0n is the slope along the boundary Op Ot is the time derivative and lt gt is the expected value W1 can be related to the expected size of boundary inflows via the geostrophic relationship fV 90p 0s where V is the expected size of the boundary inflows and 0p 0s is the slope along the boundary Thus WI Jm 10 9 V 4 It may be useful to scale V to the observed velocity data variance adjusted for the noise level Vy V snal Vise In CASCO4b there are two algorithms available for finding the minimum of the cost function steepest descent and the conjugate gradient method The conjugate gradient method generally finds the minimum in fewer iterations and is usually preferred The sele
63. m all errors which occur in time element k 1 where the dx appear in the original formalation Sections 6 1 and 6 2 Then the total contribution to the adjoint forcing the right hand side of the equations or RHS at time k is 12 Table 8 Functions which return useful numbers nelev_bcs number of elevation boundary nodes NR nbc_time_steps number of time steps in the boundary condition time series ntime steps number of time steps in a single run of saco or moody nnodes number of horizontal nodes in the mesh NN nelems number of horizontal elements NE nobservations number of velocity observations nelevobservations number of elevation observations max iterations maximum number of iterations descent_code The code to determine the descent algorithm From the point of view of the observations an error E which occurs in time element k contributes to the adjoint forcing at two different times k and k 1 The contributions to the forcing are RHS k E X RHS k 1 Ej 1 X 7 6 Other We have tried to implement a policy of only providing subroutines with the data they actually need data hiding Often a subroutine needs access to a number such as the number of boundary nodes nr but not any of the other data in the data family that nr is associated with To deal with this we have implemented functions which return many of the important numbers in the code Two examples are max_iteration and descent code in the main casco4
64. n error with elevation forcing 38 Inverse solution boundary conditions 39 Cost function as a function of iteration 40 Mode B truth 44 Mode B truth Elevation timeseries at the 6 sampling locations 45 Mode B truth Velocity timeseries at the 6 sampling locations 45 Mode B residual unexplained elevation after inversion 46 Mode B residual unexplained velocity after inversion 46 Mode B inverse solution 47 Mode B inverse solution error 48 Inverse solution boundary conditions 49 Cost function as a function of iteration 50 The structure of CASCO4b 63 The structure of SACO based on CASCO modules 65 CASCO4b genealogy 66 Definition of the model time in CASCO 70 Definition of the boundary condition time in CASCO 70 The relationship between model time and observation time 72 List of Tables D I OO O MI The input data stream for CASCO4b 3 Input parameters related to the inversion process in CASCO4b 10 List of user controlled files and subroutines 2 16 Parameters for the Mode A test case a 24 Parameters for the Mode B test case
65. n h This makes it much easier to find and change The variable scale DoDr see Section 7 6 has been moved out of the subroutine casco_input f and into the include file casco_params h This should make it much easier to find and change The subroutine statistics computes the cost function statistics In the current version the cost function components and the weighted mean sguare errors are computed separately even though they are closely related Section 3 2 This is simply a residual from earlier versions and should be cleaned up in the future 7 5 Timing 7 5 1 Model Time Time nodes in the Forward and Adjoint Models are represented with the same universal time index k used in the equations presented herein k increases in the natural forward time direction The unknown Forward and Adjoint variables are at k 1 N The Forward Model Initial Conditions ICs are at k 1 0 The Adjoint Model Terminal Conditions TCs are at k N 1 N 2 The model timestep DT is uniform Time steps are numbered in the natural manner The first time step goes from k 0 to k 1 etc Generally time step k ends at time point k time tk t0 k DT This is summarized in Fig 35 The crude notation such as t0 and DT is used in the text to allow for easy reference to the ASCII graphics in Figs 35 37 The forward model stops when the time at the end of the current time step is greater than or equal to the user specified stop time The time check is af
66. nd the data flow in CASCO4 The ovals are the primary data objects and the rectangles are processes which act on the data ABC are the perturbation boundary conditions and e represents the errors or misfits between the observations and model results models in particular the QUODDY family QUODDY4 FUNDY5 TRUXTON6 and the standard file types associated with those models The model names such as QUODDY FUNDY and CASCO are taken from geograph ical features around the Gulf of Maine they are not acronyms Nevertheless we have capitalized the names because capitalizing model names seems to be the accepted con vention 2 Care and Feeding of CASCO4b The inputs for CA5CO4b fall into five primary categories 1 The standard inputs for the hydrodynamics models including mesh files and time information Section 2 1 2 The velocity errors data which need to be minimized Section 2 2 1 2 Imaz WO W1 W2 Vioise E noise At method file ins file uxv file uxe Table 1 Convergence criteria Maximum number of iterations Regularization terms Expected rms noise level for velocity m s Expected rms noise level for elevation m Time step for the boundary conditions hours Code for the descent method Hydrodynamics input file Velocity input file Elevation input file The input data stream for CASCO4b 3 The elevation errors data which need to be minimized Section 2 3 4 The parameters which control th
67. ndary conditions bcs2 and write the output outputs2 65 CASCO Who calls Whom startup echosets3 initializes3 gday up date getfriction lin stress lin_enz station_calc stationaryc2 sprsbldi ematbuild bansoltr sprsrowzap cthomas set_times get_sacotime gday_utcO get bcstime time_basis getdata_moody getvelobs getelevobs nobservations nelevobservations getvel_obs getobs_8c basis2d utcO_gday indexx sortDRL getelev_obs getobs 4c basis2d utcO_gday indexx sortDRL measure _func saco basis basis2d poival3d vertweights poival3d getWrho BuildWrho3 sprsbldwi statistics nobservations nelevobservations nelev_bcs rms_2vec rwms_2vec rms_ivec dWDELTAd xWRHOx xWRHOy moody run up date2 EMATpreMULT vertavg xWRHO elevationm3 sprspremltadd bansoltr outputm2 tsimoody write Domega2 write Donegal write verticalm4 sprspremlt cthomas trialbcs sup norm2d trial bcs_cgm sup norm2d directionCGM saco_run bcs2 dmy vertavg up date elevations4 sprsmltadd bcs2 EMATMULT bansoltr outputs2 tsi saco write saco_uverror poival3d verticalsb sprsmlt cthomas adjust_soln sup_norm2d sup_normid nobservations nelevobservations nelev_bcs stepsize5 diWDELTAd2 dWDELTAd xWrhox xWrhoy error_to_best nobservations nelevobservations dbc_best_to_dbc nelev_bcs write_cbc dmy cstring nelev_bcs write_error write_uverror write_eerror write_uverror gday_utc0 write_eerror gday utcO Figure 34
68. ned as Observa tion Model At initialization the Model is zero so the errors are the observations read from an external file The file format used for the velocity input and output is the uxv UneXplained Velocities The file format assumes that the velocity observations have been integrated over some vertical interval The upper and lower limits of the vertical interval and the observed water depth is reguired for each data record The format was constructed for ADCP data but current meter data can be added by assuming integration over some small interval e g 1 m 2 2 1 The UneXplained Velocities uxv file The uvx file is a specialized version of the generic m3d format The format is 0 an arbitrary number of comment lines E XXXX the first required line 2 mesh name e g test6 gt WwW Year e g 1995 o comment 8 the number of columns 6 An arbitrary number of data records Each data record a line consists of 8 items e time decimal days in UTCO format days from beginning of the year This is yrday0 utc in the USGLOBEC Georges Bank Data Thesaurus e xobs yobs horizontal positions in same coordinate system as the FEM mesh Units meters e uobs vobs horizontal velocity in E N coordinates Units m s e zminobs upper limit of the observed water column Units m e zmaxobs lower limit of the observed water column e zbotobs observed water depth Notes
69. oundary condition time j is computed in set_times as follows For each k we need 1 the next j in the list i e min j s t time k lt time j 2 theta defined as follows e let the time from j 1 to j be DTBCS e let the time from j 1 to k be DT e let theta DT DTBCS e thus if time k time j then THETA 1 We have insisted on e time k 0 time j 1 e time k N time j Jmax e For time k lt time j 1 we set j 1 THETA 1 e For time k lt time j Jmax we set j NTBC THETA 1 The mapping from the model time k to the boundary condition time j is written to fortran logical units 80 and 81 by subroutine set_times The boundary condition time and the model time are only guaranteed to coincide at the endpoints Linear interpolation is used in all cases to map the boundary elevations between the two time bases 7 5 3 Measurement Time The mapping of the time of the observations to the universal model time k is the function of the subroutine measure_func The actual work is done by saco basis velocity and saco_basis_elev elevation The information is stored in the variables BASISOBS and EBASISOBS which are defined in weights h The description below uses the velocity variable names but the description is applies to both types of observations In general the observation time is not at the k s but fall in the intervals between These are the time elements which are numbered in the natural sense of the forw
70. red per iteration 56 6 4 Adjoint Equations There is a natural correspondence among the prime forward and dual adjoint vari ables L 50 V y 51 g lt v 52 Ub y 53 It is useful to recast the adjoint equations by a change of nomenclature Wr Via N Via F d ZT A Y 54 Sk Ws Pe Sa Ws Qr k 1 N Sk M Skri A Ga2 B 55 VE Va TIO ek W Tk ekpa We Ur k LN dra Cru C k O N 1 56 vbi Gi D k 0 N 1 57 with Terminal Conditions L Ek Ck Vi 0 k gt N 58 and Gradient Do p W y Cr Ex 59 P k 1 In this form the strong correspondence among forward and adjoint equation sets is more visible Essentially all the algorithms for solving the forward model can be recycled for the adjoint with due regard for the detailed differences e g reversal of first derivatives and the C V coupling 57 6 5 Covariance Matrices 6 5 1 Elevation Elevation measurements are assumed to be sampled at a point in y t Modeled elevation is sampled in the same way Elevation errors are assumed to be homogeneous and independent W is configured to compute the weighted mean sguared elevation error 2 2 y Et oise where 2 is the expected value of the squared error the elevation noise variance and Noss is the number of elevation observations Accordingly W is diagonal W e
71. res This section discusses several odd features which may cause some confusion These features are a legacy from the development pathway taken with CASCO While annoying these are not bugs The following features are a result of the fact that the CASCO input routines are based on the QUODDY4 1 input routines and the conversion has not been pushed to completion e The lev file is required even though CASCO does not use baroclinic pressure gradients e The bottom friction parameters NLBS and AKMIN are in the ins file However their only use is in get_frictionl where NLBS determines whether to use AKMIN for the bottom friction parameter or call the user subroutine lin_stress The friction related parameters should probably be relegated entirely to the user subroutines However the present structure parallels the QUODDY structure and has proven very useful for development purposes e The variable EPSN is the time weighting factor for the solution of the vertical structure of the velocity It is specified in line 24 of the ins file Fig 2 In CASCO3e and CASCO4b this variable is not used In subroutine stationaryc2 ESPN 0 5 is specified This input parameter is required for the software to run but it is ignored 68 These types of features should be removed as the software develops One item that was on this list for CASCO3b Hannah and Lynch 1999 was the variable CLOSURE This variable has been moved into the include file frictio
72. s are DELTTBCS the time interval for the boundary conditions in DBC calculated in set_times based on the averaging time scale specified by the user NTBC the number of time levels for the bcs that MOODY passes to SACO This is calculated in set_times Error Covariance The error covariance structure is assumed and diagonal and in the present implementation is assumed constant Section 2 4 The expected rms ve locity error and elevation error are specified by the user and the error variances are calculated in get_errvar CASCO4b The data structures are defined in covari ance h 67 BC Covariance The boundary condition covariance is defined by the three numbers WO WI W2 and the calculation in BuildWrho3 get_Wrho W0 provides ampli tude control W1 provides slope control and W2 provides tendency control The data structures are defined in wrho h Friction The friction parameters are user specified parameters defined by calls to user subroutines bottom friction lin_ak vertical eddy viscosity lin_enz By assumption the friction is constant in time The horizontal viscosity has been eliminated Vertical Grid The vertical node positions ZMID are computed using using vertgrids2 initializes3 The vertical grid is fixed in time I think all references to ZOLD and ZNEW have been removed One conseguence of fixing the vertical grid is that ZMID I NNV cannot be used to determine the elevation 7 4 Legacy Featu
73. t un modele inverse dont le but est de d terminer les conditions de fronti re optimales et d l vation variables dans le temps partir des observations de vitesse et d l vation l int rieur du domaine du mod le Ceci est fait en minimisant la diff rence entre la vitesse et l l vation model es et observ es Le mod le est bas sur une impl mention l ments finis des quations tridimensionnelles lin aires et barotropiques du mouvement et de la m thode d adjoint CASCO4b contient deux am liorations importantes par rapport son pr d cesseur CASCO3e il permet l inversion de l l vation et de la vitesse tandis que CASCO3e prenait en compte seulement les observations de vitesse et 2 la condition de fronti re c ti re coulement nul normal la c te a t introduite Ce rapport est une description exhaustive de CASCO4b et inclut une description compl te de la structure d entr e et de sortie trois exemples avec l ensemble complet des fichiers d entr e et des sorties graphiques une description de la th orie et des d tails de l impl mentation en code de Fortran vil 1 Introduction CASCO is an inverse model whose purpose is to determine the optimal time varying elevation boundary conditions from velocity and elevation observations in the interior of the domain This is done by minimizing the misfit between the modelled and observed velocity and elevation CASCO has 3 major components
74. ten to a tsl file The variable samp_dt defines the time step for the output time series The five lines after Output files define the basic output files The first three file names are required the cbc file the final error at the observation locations from the final iteration and the error from the final SACO run The final errors should be identical and this provides a consistency check The last two lines are part of the user subroutines samp_for_root is the root file name for the time series output for SACO the forward model during the iteration process In this case the time series output during iterations 4 and 51 for example will be written to the files samp_saco 04 ts1 and samp_saco 51 ts1 samp_for_root is the root file name for the time series output for MOODY the adjoint model In this case the time series output during iterations 4 and 51 for example will be written to the files samp moody 04 ts1and samp_moody 51 ts1 34 The nine lines after CASCO Parameters define the inversion parameters reported in Table 4 The line after The executable program defines the location and name of the pro gram which will be run The lines after Run CASCO actually run the code The first line executes the program The lt EOD is a c shell feature which says read the following lines until you find the string EOD and treat the lines as if they were typed in from the keyboard The pair of EOD s can be replaced by any ot
75. ter the time step has been completed Thus the Forward Model will generally run beyond the requested stop time by a fraction of a time step This is turn determines the location of the Terminal Conditions for the Adjoint Model since the time nodes coincide exactly The real world time is kept track of by two numbers kd seckd kd is the day counter and seckd counts the elapsed seconds in the current day 7 5 2 Boundary Condition Time The MOODY boundary conditions which are sent to SACO are defined at a larger timestep than the model time step identical in both SACO and MOODY Let the boundary condition time counter be j Figure 36 shows the relationship between j and the universal time counter k 69 t tk t0 k DT 0 at the beginning of the forward model run 1 at the end of the first forward model time step N at the end of the last forward model time step N 1 at the beginning of the adjoint backwards model run N at the end of the first adjoint model time step 1 at the end of the last adjoint model time step Figure 35 Definition of the model time in CASCO k is the universal time index ICs are the initial conditions and TCs are the terminal conditions k 1 N 1 N N 1 N 2 j j 1 1 coincides with k Jmax coincides with k Figure 36 Definition of the boundary model time in CASCO k is the universal time index Fig 35 and j is the boundary condition time index 70 The mapping between the model time k to the b
76. the mesh and the time window 4 CASCO4 stops if the mesh name in the uxe file is not identical to that in the ins file 5 The data does not have to be sorted with respect to time as CASCO4 sorts the data internally However the output will be sorted 2 4 Inversion Parameters The inversion process reguires parameters distinct from those reguired by the hydrody namics These include information on e The expected value of the errors associated with the input velocity time series e The expected value of the errors associated with the input elevation time series e The information for construction of the regularization terms reguired for an un derdetermined system e The convergence criteria e The descent algorithm The 10 parameters required by CASCO4a Table 1 are shown in Table 2 with more details Further details are given below The cost function Q which is minimized in CASCOA is n ln RR Y cw 1 where p are the boundary elevations W contains the regularization terms are the velocity errors Ws is the covariance matrix for the velocity errors ex are the elevation errors W is the covariance matrix for the elevation errors The construction of the covariance arrays for the velocity errors and the elevation errors assumes that the errors are uncorrelated and egual The resulting covariance matrix W is diagonal and requires a single input the velocity noise level V oise 1 1 wre 2 noise
77. tical eddy viscosity ekmin minimum value for the vertical eddy viscosity Kd seckd Gregorian day and elapsed seconds in the day ncall O at first call ncall lt gt 0 otherwise C Called once at startup COO ORO OG He I I EEE 2 2k 4 24 21 2k 24 3 21 2k 1 3 2k 24 2K ak 2k 24 2 ak 2k 3 2k Ek kak 2k kkk af af ak Figure 5 Definition of lin_enz 18 Ok ke ke ak ak ak ak 3k 3k 3k 3k 3k HE le le 3k EEE DE ak ak DE HE EC ccc EEE 2K K ak K ak 2K 2k 2K HE ff 2k 2K 2K 2K SUBROUTINE VERTGRIDS2 NN NNV HDOWN ZETA Z C Include file for parameter statements INCLUDE CASCO DIM C Fixed global variables INTEGER NN NNV REAL HDOWN NNDIM ZETA NNDIM REAL Z NNDIM NNVDIM C Pron nn nn nn nn nn e C This subroutine must be modified by the user to assign the vertical C grid under each horizontal node There must be assignments to C specify the location of the nodes in the 3 D array Z 1 J C C Index range I 1 NN J 1 NNV C NN number of horizontal nodes C NNV number of vertical nodes C C Z and J are positive upward C Bottom Z HDOWN J 1 C Surface Z ZETA J NNV C C HDOWN I is the water depth positive down C ZETA I is the elevation and should be identically zero for CASCO3 C Called once at startup Figure 6 Definition of vertgrids2 a nonlinear inverse then perhaps the vertical grid should be taken from the nonlinear circulation model However this has not been an issue during development 4 3 Output
78. tures and define one or two common blocks per family Each family is defined in an include file 62 CASCO4b Main Module iter 0 call casco_input itermax max_iterations method descent_code call startup call station_calc call get_friction iter call set_times call getdata_moody call measure_func call get_errvar call get_Wrho call initial_conditions call statistics iter converged false Get the casco user inputs IEstablish the maximum of iterations Get descent method Read the mesh files and such ICalculate all the constant and common arrays lRead the friction values IPrecalculate all the timing information Get the forcing data Compute the measurement functionals Compute the error variances IBuild weight function for the elevation bcs ISet the initial conditions ICompute initial value of the cost function do while iter lt itermax and not converged iter iter 1 call moody_run iter if method 1 then call trial_bcs iter else if method call trial_bcs_cgm iter else STOP due to error endif call saco run iter call adjust_soln iter call statistics iter ICompute the gradient of the cost function Compute the trial bcs for the saco run IGradient descent IConjugate gradient method Compute velocity errors based on trial bcs Compute optimum step size bcs amp errors converged convergence iter Test for convergence end while call write_c
79. ty error at the moorings Essentially where the BC s exert little influence over the observations the regularization terms dominate These BC s in turn create the most serious extrapolation errors The cost of the elevation errors is identically zero since there were no elevation observations The details of the solution in the regions with large errors are noticeably different in CASCO4b compared with CASCO3e These changes result from the changes in the details of the descent algorithm in CASCO4b Section 7 Since the velocity misfit is well below the specified error level the program is clearly rummaging around trying to minimize the regularization terms 26 Mode A Truth uvbar elevation cm green Ele ese maine PE n el LEI o fut Sor De io sabe On oe Aa ie SA EA La e o Cd o e co sn lu SN de SD o Cal ds ea A o gt o LO o LO o Kr H S wy p 3 2 p O 9 50 L 100 150 200 100 50 x km 100 50 x km Figure 10 Mode A truth the steady state 27 Velocity cm s I 5 D ON o O1 o 100 150 200 Time hours Figure 11 Truth velocity time series for Mode A at the 6 sampling stations shown in Fig 9 Small damped inertial oscillations are superposed upon the guasi steady solution The sampling is perfect i e no observational error 0 6 T T T T T lt
80. ut since files are generated during each iteration CASCO4b writes some diagnostic information to fortran logical unit numbers without defining the file names This results is the system generating default file names which vary between compilers The timing information related to the perturbation boundary conditions is written to logical unit 81 in Gregorian data format Information about how to map the hydrody namic model time to the boundary condition time is written to logical unit 80 For example the information about the measurement functionals for the velocity ob servations and the velocity covariance matrix Ws are written to logical unit 22 The information on the measurement functionals for the elevation observations can be write to a file by uncommented the relevant lines in saco_basis_elev 15 Array Sizes CASCO DIM Include file that defines the array sizes User Subroutines lin_stress bottom friction lin_enz vertical eddy viscosity vertgrids2 vertical nodes outputs2 output for SACO runs outputm2 output from MOODY runs Auxillary Output Subroutines tsl_saco_write time series of elevation and velocity tsl_moody_write time series of adjoint elevation and velocity domega2_write writes 00 0p domegal write writes 00 0 Table 3 List of user controlled files and subroutines 4 User Defined Subroutines The list of files and subroutines which the user needs to control is given in Table 3 The file CASCO DIM contains the d
81. where NN number of horizontal nodes ak I linear bottom friction parameter akmin minimum value for the bottom friction parameter kd seckd the Gregorian day and elapsed seconds in the day ncall 0 at first call ncall lt gt 0 otherwise Called once at startup acao aacooaoaoaoaaoo Cost 4 2 kc 4 24 I TTI TTI CC IC a II I I 3 2 36 ok 2 2k a kk kk kk k 2 kkk kk k tol ak kkk Figure 4 Definition of lin_stress 17 Cost 4 2 kc 4 24 I TTI TTI CC IC a II I I 3 2 36 ok 2 2k a kk kk kk k 2 kkk kk k tol ak kkk C Read vertical eddy viscosity coefficients from data files SUBROUTINE lin_enz kd seckd ncall enz enzmin nn nnv ARAQAQAQAAARA AA QAR a cacraooaooo C Problem G4 and thus SACO use the vertical eddy viscosity as a vertical element based guantity F5 and thus TRUXTON use a vertical node based vertical eddy viscosity Thus a different en s3r is reguired for the different applications This code does not distinguish since in Q4 ENZM is dimensioned ENZM NNDIM NNVDIM and the same reading code can serve both masters IF the input is formatted correctly That is if in 04 applications the element values are written to the first NNV 1 positions and a zero is inserted at NNV INCLUDE CASCO DIM integer ncall integer kd real seckd real enz nndim nnvdim real enzmin integer nn nnv Index range I 1 NN where NN number of horizontal nodes J 1 NNV where NNV number of vertical nodes enz I J ver
82. x Dvb A LD Hk 1 34 00 x x D r Ws Jk 35 00 x Ze da ek We Aj 36 54 a po YA 37 Setting the first six of these to zero gives the Adjoint Model Yk Ymi N Vi FT Dex ZT Uin Y 38 Og Ws Pr d Ws Qr k 1 N Ak M Akpa A Ak42 B 39 Ve F Vka Z G ee Wel T char We Ora k LL Pest il C k O N 1 40 DL 41 D k 0 N 1 41 with Terminal Conditions Oks Ek k Yk 0 k gt N 42 Note that and have been eliminated The last gradient vector 00 0p equation 37 gives the direction of steepest de scent It is set to zero by iteration During iteration it directs the improvement in the boundary conditions p 6 3 Gradient Descent Algorithm We search for the minimum Q in the parameter space p From any given point p there are 2 decisions which direction to move in p and how far to go in that direction 0 before changing direction p p p 43 Any selection of p will produce new errors 6 06 44 e de 45 99 with 6 e the current errors and cT the change in the errors achieved with 9 1 In turn the objective function will become a simple function of 9 The optimal is obtained by simple minimization p PU Wolo 0 Wild e Wde 46 AA w e This minimizes Q along the search direction p The direction of descent p is selected according to th

The CASCO4b User's Guide

Contents

Download Pdf Manuals

Related Search

Related Contents