For questions and comments concerning HAMSOM please send an e-mail to huebner(at)ifm.uni-hamburg.de.

20.06.2002

Johan van der Molen wrote:

Dear Jan Backhaus,

a few years ago, Meinte Blaas obtained a version of your HAMSOM ocean model for research purposes at our institute. We are presently developing an addition to the model to compute sediment transport and sea-bed change (morphodynamics). The first results are very promising: we have obtained a stable code that is capable of generating realistically looking seabed features in an idealised coastal sea geometry. Unfortunately, there is a small programming or discretisation error in HAMSOM that influences the results. Could you please let us know if this is a known bug, and if there is a remedy?

The problem is that HAMSOM generates slight deviations in the flow velocities (both u and v) near closed (dry) boundaries. A morphodynamical code is very sensitive to this kind of errors. The model is run in the 2D mode (1 layer), with only tidal forcing. Density and wind forcing are excluded. We have carried out a number of tests by subsequently switching off terms in the equations of motion as allowed by switches in hamsom.input. The problem remains when horizontal diffusion and advective terms are switched off. The velocity deviations seem to disappear when Coriolis forcing is switched off (we use the Wais rotation matrix option, switching this off makes the model numerically unstable). We therefore suspect that the error is in the discretisation of the Wais matrix near dry cells, but we cannot yet identify an error in the source code.

I hope that you recognise this problem, and have found a remedy. If not, could you please advise us on how to find its cause?

With best regards,
Johan van der Molen

Our idea about the computation of the pressure gradients was right, the problem is solved. The following piece of code gives an example of the repairing IF statement. Such a statement should be inserted twice in MOTONE, and twice in MOTTWO.

For all other geometrical situations, the original code should be correct (if no coriolis is computed for the case of a coast with a 1 cell wide entrance to a channel, i.e.: only three wet cells out of 6).

We have not tested it for the case of density anomalies (which are taken into account in MOTONE, but not in MOTTWO !???), or 3D flow.

With best regards,
Johan

IF (iwais)
THEN! external press. gradient (according to Wais)
pyv = fdvy(i)*sinfvcc(i)
pyu = cosfvcc(i)*cfiv1(i)/(dpx(k)*dlvu(i))
dpv = pyv*(tmp(i,k)-tmp(ip1,k))
facpe = 1./REAL( & jc(1,ip1,k)+jc(1,ip1,kp1)+jc(1,i,k)+jc(1,i,kp1))
facpw = 1./REAL( & jc(1,ip1,k)+jc(1,ip1,km1)+jc(1,i,k)+jc(1,i,km1))
pxeast = facpe* & dble(tmp(ip1,k)+tmp(ip1,kp1)+tmp(i,k)+tmp(i,kp1))
pxwest = facpw* & dble(tmp(ip1,k)+tmp(ip1,km1)+tmp(i,k)+tmp(i,km1))
dpu = pyu*(pxeast-pxwest)
IF (abs(facpe-0.5)<1E-4 .or. abs(facpw-0.5)<1E-4)
THEN !bugfix JM/JG 20-6-02: correct gradient along coasts
dpu=dpu*2.
ENDIF
! end bugfix
dpv = dpv-dpu
ELSE ! external pressure without rotation matrix
dpv = fdvy(i)*(tmp(i,k)-tmp(ip1,k))
ENDIF