WET_DRY and Dcrit

[kate]

The only code for keeping the water depth to at least Dcrit is here in step2d*.h:

Code: Select all

# ifdef WET_DRY
!  Modify new free-surface to Ensure that depth is > Dcrit for masked
!  cells.
# endif

# if defined WET_DRY && defined MASKING
          zeta(i,j,knew)=zeta(i,j,knew)+                                &    
     &                   (Dcrit(ng)-h(i,j))*(1.0_r8-rmask(i,j))
# endif

except that in practice, this code does nothing where rmask is one, including all the wet and dry cells.

[jcwarner]

good. that is what it is supposed to do. it is only to catch places where the landmask=0. if a cell has a landmask = 0 (ie it is land), then we can not let the water level go to 0. So that line of code makes the water level be Dcrit higher than the land. This is required because we can not let zeta go below h. When that happens, then the sigma levels become negative and a lot of bad things happen.

[kate]

Fine, but what's keeping the water column depth from going negative where rmask is one? Bad things happen there too.

[jcwarner]

well, where rmask = 1 is more complicated. but that is the whole wet dry schema. when h+zeta < dcrit, then the ubar/vbar are not allowed to flow out of that cell. this happens in step2d, with a call to wetdry.F. It sounds like i should take a look at one of your setups to dig in and see what is happening. is there an easy way to do this?
-john

[kate]

This domain has a lot of files and I'm still hand-tuning the bathymetry. It's just as easy for me to dig in and watch what's going on. I've mostly been concentrating on the 3-D fields, but it's time I learn what's going on in 2-D. I'm adding a print statement so I can find out where the water is leaking out.

[kate]

With all of John's bits of WET_DRY code I do indeed not get much shallower than Dcrit. It can go a little shallower due to the rzeta term.

Back to the question of tracers in the dry cells. John's code prevents the tracers from updating during the predictor part of the timestep, but they still get advection during the corrector part of the timestep from a nonzero Huon,Hvom even where umask_wet and vmask_wet are zero. To prevent that, you can try:

Code: Select all

diff --git a/ROMS/Nonlinear/step3d_uv.F b/ROMS/Nonlinear/step3d_uv.F
index 30c7058..a511e21 100644
--- a/ROMS/Nonlinear/step3d_uv.F
+++ b/ROMS/Nonlinear/step3d_uv.F
@@ -1251,6 +1251,9 @@
         DO k=1,N(ng)
           DO i=IstrP,IendT
             Huon(i,j,k)=Huon(i,j,k)-DC(i,k)*FC(i,0)
+#  ifdef WET_DRY
+            Huon(i,j,k)=Huon(i,j,k)*umask_wet(i,j)
+#  endif
           END DO
         END DO
 !
@@ -1470,6 +1473,9 @@
           DO k=1,N(ng)
             DO i=IstrT,IendT
               Hvom(i,j,k)=Hvom(i,j,k)-DC(i,k)*FC(i,0)
+#  ifdef WET_DRY
+              Hvom(i,j,k)=Hvom(i,j,k)*vmask_wet(i,j)
+#  endif
             END DO
           END DO
         END IF

The application of srflx in pre_step3d.F can also be masked out in the dry cells.

With these changes, the dry cell tracers remain completely static unless there is flow into the dry cell. In practice, there's some instability whereby always-dry cells can gain salt through the intermittent salt flux on their wet side. Perhaps it is the model's way of telling me I should have masked those cells.

If the cells never get wet, why should I care? Their wacky tracer values infect the neighboring wet cells through horizontal diffusion.

[nacholibre]

Hi Kate,
Did this change in step3d_uv.F work to prevent the tracers from updating in your case? I was wondering whether this worked or did not work and this change was rolled back, because I did not see it in the your repository.
Thank you
Zafer

+# ifdef WET_DRY
+ Hvom(i,j,k)=Hvom(i,j,k)*vmask_wet(i,j)
+# endif

+# ifdef WET_DRY
+ Huon(i,j,k)=Huon(i,j,k)*umask_wet(i,j)
+# endif

[kate]

It turns out I want the tracers to update in the dry cells. Who knew? Well, I let them partially update - I hacked in some code to prevent them from diffusing with their wet neighbors due to an instability with the rotated mixing tensors. I got the most stable run by turning off horizontal diffusion globally...

[arango]

I can see how the rotated tensor for tracer diffusion can go unstable near wetting and drying areas. We cannot numerically develop a wetting and drying algorithm from first principles in a model like ROMS, but we can use an engineering approach. I think that if you need horizontal diffusion in your application, you can avoid diffusion in areas where the wetting and drying is taking place by activating the tracer sponge switch LtracerSponge and adding the diff_factor to the application GRID NetCDF file. Recall that now in ROMS the diffusion coefficients are computed as:

Code: Select all

                  diff2(i,j,itrc) = diff_factor(i,j) * diff2(i,j,itrc)
                  diff4(i,j,itrc) = diff_factor(i,j) * diff4(i,j,itrc)

So you can manipulate the values of diff_factor regardless if you want an sponge or not. For example, you can set diff_factor to be zero in areas where the bathymetry is shallower than an specified value (say, h < 50m) and a value of one anywhere else. It is also possible to have a sponge next to the open boundary, if so desired.

See the following trac ticket fro more details.

[patrickm]

The implicit correction of Lemarié et al. (2012) fixes the problem for us (in ROMS_AGRIF). We implemented it for both interior diffusion (bi-laplacian) and sponges (laplacian) and WET_DRY runs without any problem in all configurations (e.g., with wetting-drying cells in sponge areas).

[wpgong]

Hi, Patrick:
could you tell where and how you did the modification in the source code to successfully run with wetting and drying?
thanks.

10/15

[nacholibre]

Thanks to everybody who contributes to the discussion. I have a simulation which blows up frequently and/or I see unreasonably high tracer concentrations in shallow water and that is why this topic interested me.

Even though we updated the code earlier to apply masking to dry cells with BULK_FLUXES the SOLAR_SOURCE was overlooked, since it is applied separately in pre_step3d.F. I think that is a good catch. Thank you, Kate. But the issue was still there without the SOLAR_SOURCE.

Using SPONGE (or activating LtracerSponge for the latest version) would be something to look into but in my case I don't have any horizontal diffusion of tracers defined in my domain (TNU2=0.0).

Lemarié et al. (2012) sounds to be a promising solution to this issue, but I assume it requires a senior ROMS developer to undertake that responsibility, though I would be happy to get involved and help if I can. From what I understand, this approach relaxes the constraint due to the gentle slope approximation and a provides a better overall solution to the tracer instability over the large slopes. I also posted about this at another topic (See:wetting and drying problem). I assume this will help with the cases with frequent wetting and drying in the domain as well, where this issue is more common. Additional treatment might be necessary for drying cells (such as Kate's approach to not updating them or else).

I would be happy to hear more on this topic.
Zafer

[jcwarner]

?

[kate]

Have you downloaded the latest Rutgers ROMS patch on this? It will probably be next week before I can get through all the merge conflicts to try it out.