Subscript error regarding BED array in stratified cohesive sediment application (with COAWST)

[c.drinkorn]

Hi all,

I am experiencing a strange error regarding the number of bed layers in a stratified cohesive sediment application:

Subscript #3 of the array BED has value 31 which is greater than the upper bound of 30

The error is traced back to:

coawstG 0000000001B1CFCF sed_bed_cohesive_ 581 sed_bed_cohesive.f90
coawstG 0000000001AC5196 sed_bed_cohesive_ 66 sed_bed_cohesive.f90

where the first is the line where tau_cr is being updated:

Code: Select all

                  bed(i,j,Ksed,ibtcr) = bed(i,j,Ksed+1,ibtcr)-                 &
     &                 cff2*(bed(i,j,Ksed+1,ibtcr)-bed(i,j,Ksed,ibtcr))

and the latter is the call to sed_bed_cohesive_tile, respectively.

The setup has 30 bed layers so 31 are indeed not correct. The error only occurs with cohesive sediment.
This is my initial setup in ana_sediment:

Code: Select all

!  Set bed layer properties.
!
          DO k=1,Nbed
            bed(i,j,k,iaged)=time(ng)
! Each layer is 1cm thick and has porosity of 0.66
          bed(i,j,k,ithck)=0.001_r8
        bed(i,j,k,iporo)=0.66_r8
            DO ised=1,NST
              bed_frac(i,j,k,ised)=1.0_r8/REAL(NST,r8)
            END DO
          END DO

and here is my cohesive sediment section in sediment.in:

Code: Select all

!------------------------------------------------------------------------------
! Cohesive / Mixed bed critcal shear
! Bed critical shear stress values if cohesive or mixed Bed is active
! These values are applied to the entire bed. [1:Ngrids] values expected
!------------------------------------------------------------------------------

! Minimum shear for erosion

MUD_TAUCR_MIN = 0.01

! Maximum shear for erosion

MUD_TAUCR_MAX = 1.0

! Tau_crit profile slope

MUD_TAUCR_SLOPE = 4.0

! Tau_crit profile offset

!MUD_TAUCR_OFFSET = 1.0
MUD_TAUCR_OFFSET = 12.0

! Tau_crit consolidation rate Tc

MUD_TAUCR_TIME = 8800 !16 hours

I hope there is someone who recognizes this error or has an idea where it comes from. I have tried for a while to solve this myself but I did not succeed.
Thank you!
Cate

[aretxabaleta]

Looks like an interesting problem.
In the cohesive bed routine, computations happen in both the Ksed layer and in the Ksed+1, but that should not exceed the number of bed layers.
“! Ksed is used for grain props, and Ksed+1 for bed_tau_crit at bottom of layer Ksed”
The error message seems to be happening when the layer info gets updated after erosion has happened.
Two thoughts:
- Maybe an issue happens when the entire bed is eroded and the Ksed reaches the maximum bed layer. I think the code has a catch for this already, but I wonder if it is failing. Could you please check to see if you are eroding the entire bed?
- It seems that it happens during debug mode, I wonder if it is something that gets skipped during normal compilation and that is why we have not seen it.
Could you also check what your minimum layer thickness (newlayer_thick in sediment.in) is?

a

[c.drinkorn]

HI,

yes, the problem only occurs in Debug mode. In normal Mode the model runs into stability problems after some time.
Unfortunately, I can not look into what happens with the bed layers because the error occurs already at the first time step so that both rst and dia are empty. In the first his file everything still looks normal, of course.
newlayer_thick is now set to zero but it was set to all kinds of values from 0.001 to 0.003 before, with no effect on the outcome.
I will try to run the model with much thicker layers that would prevent any complete erosion and see if this makes a difference...

Thank you for your help!

Cate

[c.drinkorn]

UPDATE:

Ok, so I ran into the same error with layer thicknesses of 30cm, i.e. 30 layers à 30cm which makes 9m of sediments. I doubt that it was possible to erode this within one time step... It seems as it is rather a numerical problem instead that may have been overlooked because it is ignored / compensated for in normal simulation mode...

[aretxabaleta]

Could you please modify a line in the sed_bed_cohesive.F code?
Could you change line 705 from:

Code: Select all

DO k=2,Nbed

to

Code: Select all

DO k=2,Nbed-1

That specific loop should now look like:

Code: Select all

      
      DO k=2,Nbed
                IF (thck_avail.lt.thck_to_add) THEN
                  thck_avail=thck_avail+bed(i,j,k,ithck)
                  Ksed=k
                END IF
      END DO

That should deal with the debug issue and then we can figure out why you encounter an issue at runtime

[c.drinkorn]

That was a really good idea! The result is stunnning:
Of all 30 layers à 30cm thickness, the lower 29 have a thickness of zero at the first time step already! The top layer maintains a thickness of several meters (EDIT: I just realize it now after I have read it another time - several meters!? OK, so it seems there is something wrong with the dynamics of pushing layers up maybe?) except for some scattered small regions. The model crashes after 4 time steps of 10 minutes with the error:

[m11000:6671 :0] Caught signal 8 (Floating point exception)
...
==== backtrace ====
2 0x00000000000575cc mxm_handle_error() /scrap/jenkins/workspace/hpc-power-pack/label/r-vmb-rhel6-u7-x86-64-MOFED-CHECKER/hpcx_root/src/hpcx-v1.8.0-gcc-OFED-3.18-redhat6.7-x86_64/mxm-v3.6/src/mxm/util/debug/debug.c:641
3 0x000000000005773c mxm_error_signal_handler() /scrap/jenkins/workspace/hpc-power-pack/label/r-vmb-rhel6-u7-x86-64-MOFED-CHECKER/hpcx_root/src/hpcx-v1.8.0-gcc-OFED-3.18-redhat6.7-x86_64/mxm-v3.6/src/mxm/util/debug/debug.c:616
4 0x0000000000032570 killpg() ??:0
5 0x0000000001aebdb3 sed_bed_cohesive_tile() /mnt/lustre01/scratch/b/b380636/Arctic20km/obcatmforceriversinisedbulk100_clim_bering5/Build/sed_bed_cohesive.f90:344
6 0x0000000001ac5196 sed_bed_cohesive() /mnt/lustre01/scratch/b/b380636/Arctic20km/obcatmforceriversinisedbulk100_clim_bering5/Build/sed_bed_cohesive.f90:66
7 0x0000000001a75c72 sediment() /mnt/lustre01/scratch/b/b380636/Arctic20km/obcatmforceriversinisedbulk100_clim_bering5/Build/sediment.f90:100
8 0x0000000000429cc3 main3d() /mnt/lustre01/scratch/b/b380636/Arctic20km/obcatmforceriversinisedbulk100_clim_bering5/Build/main3d.f90:383
9 0x000000000040e982 roms_run() /mnt/lustre01/scratch/b/b380636/Arctic20km/obcatmforceriversinisedbulk100_clim_bering5/Build/ocean_control.f90:181
10 0x000000000040d2dc ocean() /mnt/lustre01/scratch/b/b380636/Arctic20km/obcatmforceriversinisedbulk100_clim_bering5/Build/master.f90:108
11 0x000000000040ca9e main() ??:0
12 0x000000000001ed20 __libc_start_main() ??:0
13 0x000000000040c9a9 _start() ??:0

I looked up line 344 in sed_bed_cohesive.f90 and it is the following line:

Code: Select all

                cff1 = bed(i,j,1,ibtcr) -                               &
     &             bottom(i,j,idnet) *                                  &
     &             (bed(i,j,2,ibtcr)-bed(i,j,1,ibtcr)) /                &
     &             (bmz(1)-bottom(i,j,idnet))

So I guess bmz(1)-the net deposition becomes zero and causes the floating point exception error? So if cff1 has no solution, in the following lines where

Code: Select all

                cff2 = MAX( rho0*tau_w(i,j) , tcr_min(ng) )
                bed(i,j,1,ibtcr) = MIN(bed(i,j,1,ibtcr),                &
     &                                 MAX( 0.5_r8*(cff1+cff2), cff2 ))

ibtcr will always be the minimum value which indeed is also the case in my output.

However, I still don't understand why all the sediment layers become zero.

[aretxabaleta]

Well, I think it might be because you have newlayer_thick set to zero, so it will create new layers of zero thickness. Could you change that?

[c.drinkorn]

Yes, of course.
But now that newlayer_thick is larger than zero the layers are all exactly that value.

[aretxabaleta]

The idea of the new layer thickness is that is a way to preserve the stratigraphy of the bed after deposition has occurred and new layers need to be created. You might need to fine tune the choice of new_layerthick so that it makes sense to your application.

[c.drinkorn]

So, I could trace back the reason for this error and find a rudimentary workaround for it:

The problem in sed_bed_cohesive is that the catch for the case when the available sediment is not enough for the amount needed to add is not taking care of the variable Ksed. If all the sediments are not enough to accumulate thick_avail in the necessary amount, Ksed will reach the value of Nbed. Then, further down in the code when tau_cr is being updated, an index of i,j,Ksed+1 is being addressed. This is, of courses, beyond the limits of the bed layers defined in the model set up.
I added the line Ksed=Ksed-1 to that catch. However, this is just a workaround and might violate several conservation restrictions. So, I would be very grateful for any feedback on that!

Why the available sediments are being reduced so fast is another question but it could just be due to some spin up artefacts which would take up too much time to be looked into...

On another note: I also added a restriction to a bed(ibtcr) value of zero in sed_fluxes since this would cause division by zero. So I changed the line

Code: Select all

cff = rho0/bed(i,j,1,ibtcr)

to

Code: Select all

cff = rho0/MAX(bed(i,j,1,ibtcr),eps)

(I added the variable eps in the header.)
Since bed(ibtcr) is being initialized as 0.0 this seems to be reasonable.
Again, any feedback on that is welcome.

[jcwarner]

i think this path your are on will help, and you have identified several issues that need to be addressed.
However we are in hope that your statement:
"Why the available sediments are being reduced so fast is another question but it could just be due to some spin up artefacts which would take up too much time to be looked into..."
We hope you actually do spend time to figure this out! I think this is a key issue, and hopefully there is something simple like winds with 9999 fill value (or something crazy like that) causing large stresses. In the end, the issues you bring up need to be addressed, but if you could get the model to take a few steps, save every step to a his file (like 10 steps), and then look to see why there is such a big change - would be helpful.
thanks
j

[c.drinkorn]

Hi Alfredo and John

the initial excessive removal of sediments in my model is physical: it is due to a deep sea water fall inducing strong erosion at its bottom (they correlate with high vertical local velocities). I could stabilize this erosion a bit by shortening consolidation times and increasing the model time steps, however, sooner or later the simulation runs into the same problem. Now that I could workaround the Ksed+1 issue I can also see where this really comes from. Debug mode reveals a floating point exception error in this line:

Code: Select all

                cff1 = bed(i,j,1,ibtcr) -                               &
     &             bottom(i,j,idnet) *                                  &
     &             (bed(i,j,2,ibtcr)-bed(i,j,1,ibtcr)) /                &
     &             (bmz(1)-bottom(i,j,idnet))

So it seems that if there is just as much erosion as sediment in the upper layer the numerator will become zero. Since this is being limited so that erosion cannot exceed this amount of material further up in the code, heavy erosion can cause this error easily. What about restricting the numerator here:

Code: Select all

                cff1 = bed(i,j,1,ibtcr) -                               &
     &             bottom(i,j,idnet) *                                  &
     &             (bed(i,j,2,ibtcr)-bed(i,j,1,ibtcr)) /                &
     &             MAX((bmz(1)-bottom(i,j,idnet)),eps)

In fact, I am just trying this and will report if it helped within the next hours.

Generally, if I have a stratified sediment application with a high vertical resolution I would expect the model to allow erosion beyond layers. The stratification shall only provide the parameterization of compaction / swelling. Or maybe I am misunderstanding it?

[jcwarner]

ok. glad you are digging in.
What is a deep sea water fall? can you show a pcolor of the bottom stresses? how high do they get?
i am sure we can come up with something, but let me see what the domain looks like.
thanks,
j