SSW_LOGINT: z1 goes negative

[nganju]

I have been having problems with shallow cells, especially when using SED_MORPH to accelerate bed change. The domain is a relatively shallow subembayment of SF Bay, with tides and freshwater input. Larger SED_MORPH causes an earlier crash, in the same cell.

After some troubleshooting, I determined that I was crashing in ssw_bbl, at SSW_LOGINT because z1 became negative at a particular cell. z_r is the elevation of the rho location, and z_w is the elevation of the vertical velocity location, correct? So how can z_r(i,j,1)<z_w(i,j,0)?

I put in a catch for now, but does this mean that the rho point is under the bed? Would a lot of deposition (relative to the original cell depth) cause this? The original depth is 1.2 m, and deposition has reduced that to about 0.75 m. No cell has gone dry at the time of the crash (zeta is well over 1 m in all areas during the crash, "bath" is no less that 0.5 m). Any ideas are appreciated...

Code: Select all

#ifdef SSW_LOGINT
!
!  If current height is less than z1ur, interpolate logarithmically
!  to z1ur. (This has not been updated wrt ssw...uses outdated zo defs)
!
IF (Zr(i,j).lt.sg_z1min) THEN
            DO k=2,N(ng)
             z1=z_r(i,j,k-1)-z_w(i,j,0)
              z2=z_r(i,j,k  )-z_w(i,j,0)
              IF ((z1.lt.sg_z1min).and.(sg_z1min.lt.z2)) THEN
                fac=1.0_r8/LOG(z2/z1)
                fac1=fac*LOG(z2/sg_z1min)
                fac2=fac*LOG(sg_z1min/z1)
                Ur_sg(i,j)=fac1*u(i,j,k-1,nrhs)+fac2*u(i,j,k,nrhs)
                Vr_sg(i,j)=fac1*v(i,j,k-1,nrhs)+fac2*v(i,j,k,nrhs)
                Zr(i,j)=sg_z1min
              END IF
            END DO
          END IF

[jcwarner]

Neil-

it sounds like wet_dry is not the cause of these problems, because you state that "no cell has gone dry at the time of the crash." The morphology causes a feedback to the momentum equations by allowing h to evolve. The update to h occurs in set_depth, at the same time that the new water level is used to compute the new cell thicknesses. Both z_r and z_w are computed with the same information. So I do not see how z_r(i,j,1)-z_w(i,j,0) can be negative.
Are you sure that the instability is created at that cell, or was it caused by something else?

[nganju]

Wetting and drying is off right now...

I put in a catch and the restarted model went for another few days before a different problem emerged (rapid bed deposition in the same spot, with large negative sediment concentrations?).

Back to the first crash: bustrc and buvstrc were Nan's in the time step before the crash, and I traced that to Ur_sg, which was getting a bad "fac" and "fac2" because of z1.

I had ssw_bbl write out the values of z_r and z_w at the problem cell, and they are very close to each other, and z_r exceeds z_w after a few time steps, then the model crashes (z_r(i,j,1)=-0.723, z_w(i,j,0)=-0.7229).

Perhaps the first and second crash are related? I will keep looking...

[nganju]

Another resolved problem to report for this one...in set_scoord.F

Code: Select all

# if defined WET_DRY
      hc(ng)=MIN(MAX(hmin(ng),0.0_r8),Tcline(ng))
# else
      hc(ng)=MIN(hmin(ng),Tcline(ng))
# endif

If you are using SED_MORPH, hmin should technically be variable, but it looks like it is pulled once, at the beginning of the run. If Tcline and hmin are greater than zero, then as sediment accumulates and updated depth reduces to less than hc, you run into problems with the calculation of z1, z_r, z_w, etc in set_depth.F and ssw_bbl.h. The easiest solution is to set Tcline to zero or hardwire hc=0. This may account for some persistent blowups when using SED_MORPH, especially in rapidly depositing areas.