The model has been used in some studies, see [8,5,4,7,15,16,14,2]. In the development of the code especially the ability of the model to maintain the properties of different water masses has been in focus and in  the model capabilities in this respect are documented using observations from the SKAGEX experiment in 1990 . Model results from the present code are also compared to model results from the Blumberg and Mellor model .
In version 1.1 of the code the statement ETA = ETAP in Main.f90 has been corrected to ETAP = ETA after Yngve Heggelund found this serious bug. In version 1.2 the effects of wind and bottom stress are included in the 2-D part of modesplit.f90. This had negligible effects on the results for the Osterøy case, but when simulating the storm surge described in  it proved to be important to include these effects also in the 2-D steps.
In version 2.0 the Coriolis term averaging described in Espelid et al.  is introduced to avoid instabilities that may or will occur when applying 1/4-weights for these terms. The superbee-limiter method is also used for advection of momentum. One may select the POM-form on the viscosity terms by setting a logical. One may choose between the UNESCO-formulation as given in Gill  or the simplified equation of state found in Wang . A number of routines have been made more efficient on shared memory parallel computers by introducing OPEN MP-directives and rewriting loops. Some routines are made more efficient by introducing arrays that enable looping over wet points only. The work on parallelization and improving the efficiency has been done by Helge Avlesen, Parallab, University of Bergen.
In version 3.0, the time stepping methods are changed. In the 2-D steps, the forward-backward method is replaced with the implicit -method. For = 0.5, the method is 2nd order accurate. For = 1.0, the method is 1th order and fully implicit. In the 2-D steps, the leapfrog method is used as a predictor of values at the new time step. For less than 0.5, this method is unstable. In the 3-D steps, one may still run the model with an explicit technique (except for the vertical mixing). There is also the option of using the -method in the 3-D step. In this case, the Euler forward method is used as a predictor. The motivation for introducing more implicitness into the model is to avoid noise and growth of numerical errors and thereby allow simulations with less viscosity and more physical free waves. The cost is more computer time. If is larger than 0 in both the 2D steps and the 3D steps, the computing time will almost be doubled since the method is based on computing first a predictor and then a corrector. In version 3.0 there is also an option for using the 4th order advection scheme presented in . This method is more accurate than the default TVD-scheme with a superbee limiter. However, over- and under-shooting near fronts may occur. The method is not monotonic.
In version 4.0, the time stepping of the 3-D steps is changed. The implicit -method is replaced by a predictor-corrector method with the leapfrog method as the predictor and the Adams-Moulton 2-level method as the corrector. This method is used in ROMS [37,36] and it is also analyzed in . The method has a stability region that includes parts of the imaginary axis and the right half complex plane, and it is third order accurate. It requires the storage of some variables at three time levels, but since it is more accurate and stable than previous methods the 3D time step may be substantially increased in many applications. The phase errors are reduced, allowing experiments with lower values of viscosity. When the 3D time steps are increased, the number of 2D time steps per 3D step may be increased to keep the 2D Courant number at an acceptable level. Also in version 4.0, it is introduced as an option to use the fourth order method described in  for estimation of internal pressure gradients. In some test cases, the internal pressure gradient errors are substantially reduced when using higher order methods.