I am working with CROCO to study Hamilton Harbour, a smaller scale embayment in Ontario with length scales of 1-10 km. I’ve successfully been able to run the model by supplying a .nc grid file derived from a bathymetry file (40 meter resolution) I have, and using an analytical temperature initialization that I modified in ana_initial.F
# elif defined HAMILTON_HARBOR
# define Z0 (-10.)
# define THKNSS 1.
# define SLOPE_DEG .04 /* thermocline slope angle in degrees */
# define TEMPSURF 20.
# define TEMPBOT 10.
do k=1,N
do j=JR_RANGE
do i=IR_RANGE
# ifdef TEMPERATURE
!
! Tilted thermocline: thermocline gets shallower moving east
! z_r is the cell-center depth (negative, 0 at surface)
! Z0 is the thermocline center depth at x = xr_min
! A slope angle of 0.02 roughly translates to a 3 meter increase in Tcline depth
!
slope_val = tan(SLOPE_DEG * pi / 180.0)
t(i,j,k,1,itemp) = TEMPBOT + 0.5*(TEMPSURF - TEMPBOT)
& *(1. + tanh((z_r(i,j,k) - Z0
& - slope_val * (xr(i,j) - xr(1, 1)))
& / THKNSS))
t(i,j,k,2,itemp) = t(i,j,k,1,itemp)
! t(i, j, k, 1, itemp) = 20
! t(i, j, k, 2, itemp) = t(i, j, k, 1, itemp)
# endif /* TEMPERATURE */
# undef Z0
# undef THKNSS
# undef SLOPE_DEG
# undef TEMPSURF
# undef TEMPBOT
# ifdef SALINITY
t(i,j,k,1,isalt) = 0
t(i,j,k,2,isalt) = t(i,j,k,1,isalt)
# endif
enddo
enddo
enddo
The initialization is that of a tilted thermocline. The model runs without any compile issues or errors. However, I find that at time 0, there is an odd effect at the boundaries which set the temperature at all depths to be constant. At the next time step, all of that constant temperature changes to another constant value. The interior evolves normally, but the ends stay fixed.
Is this happening due to the way ghost points are done? I plot a section below at two time steps to show the problem.
Here is my cppdefs options:
/*
! Hamilton Harbor tilted thermocline
! ======== ======= ====== ===========
*/
/* Parallelization */
# undef OPENMP
# undef MPI
/* Grid */
# undef SPHERICAL
# define CURVGRID
# define MASKING
# define NEW_S_COORD
/*Dynamics*/
# define UV_ADV
# undef UV_COR
# define UV_VIS2
# define SOLVE3D
# define TS_DIF2
# define BODYFORCE
/* Initial and Boundary Conditions */
# define ANA_INITIAL
# define ANA_SMFLUX
# define ANA_STFLUX
# define ANA_BTFLUX
# define NO_FRCFILE
And these are the gridpoints information I provide in param.h. The bathymetry file is 345 x 146 points.
#elif defined HAMILTON_HARBOR
parameter (LLm0=343, MMm0=144, N=20) ! HAMILTON_HARBOR
These are the vertical coordinates I am working with but I’ve made Hc smaller (2m, 5m) with no resolution of this issue. I’ve also tried different values for theta_s and theta_b.
S-coord: THETA_S, THETA_B, Hc (m)
4.0d0 2.0d0 10.d0
There was unrelated issue documented here: BLOW UP error on custom bathymetry dataset - Hamilton Harbour in case someone may want more context.

