Emulating barotropic model

Dear All,

I have some queries on the workings of NEMO. Can I seek the help of the community to correct me if I’m wrong in the understanding below:

  1. I was trying to create a two-layer setup. It created the domain cfg fine, however, the bathy_meter variable are showing double the depth from the original bathy_meter file. Please refer to the parameters in the namelist_cfg below.

  2. If we try to edit the initial and boundary conditions with a constant T & S, and turn off the model parameters in the namelist that will control the generation/dissipation, theoretically we could emulate a barotropic mode if only the tides are used as forcing. Agree? In this sense it does not matter how many layers is used.

&namcfg        !   parameters of the configuration
!-----------------------------------------------------------------------
   !
   ln_e3_dep   = .true.    ! =T : e3=dk[depth] in discret sens.
   !                       !      ===>>> will become the only possibility in v4.0
   !                       ! =F : e3 analytical derivative of depth function
   !                       !      only there for backward compatibility test with v3.6
   !                       !
   cp_cfg      =  "oifs"   !  name of the configuration
   jp_cfg      =     029   !  resolution of the configuration
   jpidta      =     1942  !  1st lateral dimension ( >= jpi )
   jpjdta      =     1222  !  2nd    "         "    ( >= jpj )
   jpkdta      =      2   !  number of levels      ( >= jpk )
   jpiglo      =     1942  !  1st dimension of global domain --> i =jpidta
   jpjglo      =     1222  !  2nd    -                  -    --> j  =jpjdta
   jpizoom     =       1   !  left bottom (i,j) indices of the zoom
   jpjzoom     =       1   !  in data domain indices
   jperio      =       0   !  lateral cond. type (between 0 and 6)
/

&namzgr_sco    !   s-coordinate or hybrid z-s-coordinate                (default F)
!-----------------------------------------------------------------------
   ln_s_sh94   = .true.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)|
   ln_s_sf12   = .false.   !  Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied
   ln_sigcrit  = .false.   !  use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
                           !  stretching coefficients for all functions
   rn_sbot_min =   10.0    !  minimum depth of s-bottom surface (>0) (m)
   rn_sbot_max = 10200.0    !  maximum depth of s-bottom surface (= ocean depth) (>0) (m)
   rn_hc       =  0.0    !  critical depth for transition to stretched coordinates
                        !!!!!!!  Envelop bathymetry
   rn_rmax     =    1.0    !  maximum cut-off r-value allowed (0<r_max<1)
                        !!!!!!!  SH94 stretching coefficients  (ln_s_sh94 = .true.)
   rn_theta    =    0.0    !  surface control parameter (0<=theta<=20)
   rn_bb       =    0.8    !  stretching with SH94 s-sigma
                        !!!!!!!  SF12 stretching coefficient  (ln_s_sf12 = .true.)
   rn_alpha    =    4.4    !  stretching with SF12 s-sigma
   rn_efold    =    0.0    !  efold length scale for transition to stretched coord
   rn_zs       =    1.0    !  depth of surface grid box
                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
   rn_zb_a     =    0.024  !  bathymetry scaling factor for calculating Zb
   rn_zb_b     =   -0.2    !  offset for calculating Zb
                        !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above]
   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1)
/
1 Like

Hello,

Since you print the namzgr_sco block, I suspect you want to create a mesh having s-coordinates and 2 vertical levels (you certainly have ln_sco=.true. in your namelist). Nothing wrong about that, I even think this is the best way to go for a barotropic model rather than assuming partial steps.
Having jpkglo=jpkdta=2 is something a bit unusual, and trying to create this kind of grid, I indeed faced two problems:

  • Although one expects to have bottom_level=1 everywhere in the domain_cfg.nc file, it is equal to 2. Using jpkglo=jpkdta=3 works however fine giving the expected bathymetry and bottom_level=2 in output. There’s something wrong going on for that special case jpkglo=2. I need some time to find out and will let you know.

  • The second problem I found is linked to an online modification of the enveloping topography (which should be your “real bathymetry” with rmax=1) but only if your domain crosses the equator. This should not be done if rmax=1 (and not done in any case unless the user choose it I think). You can check in your ocean.output if you get this message:

    ' s-coordinates are tapered in vicinity of the Equator'
    

Concerning your second point, what do you mean by “generation/dissipation” ? If you just have tides (no heat neither freshwater forcing, eventually wind and pressure), this should work as is (and with a linear equation of state). That’s not optimal though since tracers are being updated anyway.

1 Like

Seems I found the problem in domzgr.F90

  1. Comment lines 1423-1431

  2. Edit line 1554 from

       IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
    

    to

       IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) )   mbathy(ji,jj) = MAX( 1, jk )
    

I’ll create a ticket and modify the code accordingly as soon as it works on your side.

Jérôme

2 Likes

yes

I am not sure if its right if I set to the below settings, I am indeed getting the 2 levels as you have described. bottom_level=2.0 for everywhere except the land where it is 0.0

   jpkdta      =      2    !  number of levels      ( >= jpk )
   rn_theta    =    20.0    !  surface control parameter (0<=theta<=20)
   rn_bb       =    1    !  stretching with SH94 s-sigma
   rn_hc       =  0.0    !  critical depth for transition to stretched coordinates

This indeed was found in the output. Is it suggested to lower rmax to 0?

understood, as such the T and S has to be forced to be constant at the initial file, and only turning on tides at the LBCs.

the domzgr.F90 I used seems to be different.

Hi,

You get bottom_level=2 with jpkdta=jpkglo=2 which is the problem ! You should have 1 everywhere except on land (=0).

The message you obtain should not appear which is corrected by the first change proposed above in domzgr.F90.

Indeed, your domzgr.F90 comes from an older version. Still changes I propose should work.
These correspond to lines 2080-2088 and 2216 in your domzgr.F90.

Jérôme

1 Like

Thanks a lot for the clear steps, allow me some time to recompile and test this.

1 Like

Hi Jérôme,

The new domain seems correct now, the bottom_level is now 1. Also I no longer have the s-coordinates are tapered in vicinity of the Equator in the ocean.output.

Thank you.

Good. Hope your bathymetry is ok too.

You could maybe crop your maximum bathymetry to something lower than rn_sbot_max = 10200 m. This could save some barotropic iterations (in case of a split-explicit free surface) or enable you to run with a larger time step (with explicit free surface). Removing very deep trenches (>6000 m) will probably not degrade that much your result.

Related to your question about initial T/S profiles, a useful tip in case of a “barotropic” experiment:

  • Choose a linear equation of state, e.g. set ln_seos = .true.

  • Although this would be a bit long to explain here, I recommend that you choose constant temperature and salinity values that match reference values used in the equation of state (hence, by construction, model density equals the model reference density). These values are hardcoded and currently are:

Tref = 10._wp
Sref = 35._wp

Other things to know:
To save some computation time, you may disable any horizontal diffusion processes on tracers (but you must keep tracer advection, in particular with a non-linear free surface).

Don’t hesitate to post again here issues on that subject (barotropic modelling with NEMO). To my knowledge, this has not been documented.

1 Like

I would test ln_seos, also noted on the reference T and S.

I would also like to document some of the settings I have toggled from a previous baroclinic run. Shown below are for the barotropic model. Feel free to share your views or if I have missed any.

From the namelist_cfg:

ln_rstart   = .false.      !  start from rest (F) or from a restart file (T)

Intent to turn off bulk formulation, as such flux has to be on

ln_flx      = .true.   !  flux formulation 
ln_blk      = .false.    !  Bulk formulation
ln_traqsr   = .false.   !  Light penetration in the ocean
ln_apr_dyn  = .false.   !  Patm gradient added in ocean & ice Eqs.
ln_traadv_OFF = .true. !  No tracer advection
ln_traadv_fct = .false. !  FCT scheme
ln_traldf_OFF   = .true.   !  No explicit diffusion
ln_traldf_lap   = .false.   !    laplacian operator

internal wave-induced mixing

ln_zdfiwm   = .false.      ! internal wave-induced mixing

Choose not to use external data and use tidal harmonic forcing

nn_dyn2d_dta   =  2 !(previously was 3)

After changing these, I realise I also had to edit parts of the the namelist_ref as it reads in amm12bdy files in the reference. Please let me know if this step is incorrect or if there is a better way to do this.

From the namelist_ref:
I had to change the variables in bdydta and namsbc_flx to 'NOT USED'

With the above changes, I have now reached the step where it crashes without an error:

 dom_vvl_sf_swp : - time filter and swap of scale factors
interpolate scale factors and compute depths for next time step

Do you have any explicit message of the problem in your job summary ? If not, you could try to compile with debug options, e.g. with array bound checking, overflow detection, etc… (examples in arch directory)
What is your version of NEMO ?

1 Like

Sorry its a bit long but it starts with the following

forrtl: severe (174): SIGSEGV, segmentation fault occurred
Image              PC                Routine            Line        Source             
nemo               00000000018AB765  Unknown               Unknown  Unknown
nemo               00000000018A9527  Unknown               Unknown  Unknown
nemo               0000000001848BA4  Unknown               Unknown  Unknown
nemo               00000000018489B6  Unknown               Unknown  Unknown
nemo               00000000017D6606  Unknown               Unknown  Unknown
nemo               00000000017D9C20  Unknown               Unknown  Unknown
libpthread.so.0    00002AAAB045D7E0  Unknown               Unknown  Unknown
nemo               0000000001010A04  Unknown               Unknown  Unknown
nemo               00000000017447EE  Unknown               Unknown  Unknown
nemo               00000000016BF482  Unknown               Unknown  Unknown
nemo               00000000013B5FDF  Unknown               Unknown  Unknown
nemo               00000000015E6A67  Unknown               Unknown  Unknown
nemo               00000000015E6D5C  Unknown               Unknown  Unknown
nemo               0000000001613BF7  Unknown               Unknown  Unknown
nemo               00000000010B36EE  Unknown               Unknown  Unknown
nemo               00000000012D0839  Unknown               Unknown  Unknown
nemo               0000000000CFCFDE  Unknown               Unknown  Unknown
nemo               000000000074C850  domvvl_mp_dom_vvl         682  domvvl.f90
nemo               0000000000481AF3  step_mp_stp_              235  step.f90
nemo               0000000000432931  nemogcm_mp_nemo_g         149  nemogcm.f90
nemo               0000000000432840  MAIN__                     18  nemo.f90
nemo               000000000043280E  Unknown               Unknown  Unknown
libc.so.6          00002AAAB090ED20  Unknown               Unknown  Unknown
nemo               00000000004326E9  Unknown               Unknown  Unknown
forrtl: severe (174): SIGSEGV, segmentation fault occurred
Image              PC                Routine            Line        Source             
nemo               00000000018AB765  Unknown               Unknown  Unknown
nemo               00000000018A9527  Unknown               Unknown  Unknown
nemo               0000000001848BA4  Unknown               Unknown  Unknown
nemo               00000000018489B6  Unknown               Unknown  Unknown
nemo               00000000017D6606  Unknown               Unknown  Unknown
nemo               00000000017D9C20  Unknown               Unknown  Unknown
libpthread.so.0    00002AAAB045D7E0  Unknown               Unknown  Unknown
nemo               0000000001010A04  Unknown               Unknown  Unknown
nemo               00000000017447EE  Unknown               Unknown  Unknown
nemo               00000000016BF482  Unknown               Unknown  Unknown
nemo               00000000013B5FDF  Unknown               Unknown  Unknown
nemo               00000000015E6A67  Unknown               Unknown  Unknown
nemo               00000000015E6D5C  Unknown               Unknown  Unknown
nemo               0000000001613BF7  Unknown               Unknown  Unknown
nemo               00000000010B36EE  Unknown               Unknown  Unknown
nemo               00000000012D0839  Unknown               Unknown  Unknown
nemo               0000000000CFCFDE  Unknown               Unknown  Unknown
nemo               000000000074C850  domvvl_mp_dom_vvl         682  domvvl.f90
nemo               0000000000481AF3  step_mp_stp_              235  step.f90
nemo               0000000000432931  nemogcm_mp_nemo_g         149  nemogcm.f90
nemo               0000000000432840  MAIN__                     18  nemo.f90
nemo               000000000043280E  Unknown               Unknown  Unknown
libc.so.6          00002AAAB090ED20  Unknown               Unknown  Unknown
nemo               00000000004326E9  Unknown               Unknown  Unknown
forrtl: severe (174): SIGSEGV, segmentation fault occurred
Image              PC                Routine            Line        Source             
nemo               00000000018AB765  Unknown               Unknown  Unknown
nemo               00000000018A9527  Unknown               Unknown  Unknown
nemo               0000000001848BA4  Unknown               Unknown  Unknown
nemo               00000000018489B6  Unknown               Unknown  Unknown
nemo               00000000017D6606  Unknown               Unknown  Unknown
nemo               00000000017D9C20  Unknown               Unknown  Unknown
libpthread.so.0    00002AAAB045D7E0  Unknown               Unknown  Unknown
nemo               0000000001010A04  Unknown               Unknown  Unknown
nemo               00000000017447EE  Unknown               Unknown  Unknown
nemo               00000000016BF482  Unknown               Unknown  Unknown
nemo               00000000013B5FDF  Unknown               Unknown  Unknown
nemo               00000000015E6A67  Unknown               Unknown  Unknown
nemo               00000000015E6D5C  Unknown               Unknown  Unknown
nemo               0000000001613BF7  Unknown               Unknown  Unknown
nemo               00000000010B36EE  Unknown               Unknown  Unknown
nemo               00000000012D0839  Unknown               Unknown  Unknown
nemo               0000000000CFCFDE  Unknown               Unknown  Unknown
nemo               000000000074C850  domvvl_mp_dom_vvl         682  domvvl.f90
nemo               0000000000481AF3  step_mp_stp_              235  step.f90
nemo               0000000000432931  nemogcm_mp_nemo_g         149  nemogcm.f90
nemo               0000000000432840  MAIN__                     18  nemo.f90
nemo               000000000043280E  Unknown               Unknown  Unknown
libc.so.6          00002AAAB090ED20  Unknown               Unknown  Unknown
nemo               00000000004326E9  Unknown               Unknown  Unknown

Let me figure out the debugging options, my version is 4.0.6.

Something goes wrong at line 682 of domvvl.f90 (located in YOURCONFIG/BLD/ppsrc/nemo/). Have a look at that line.
Meanwhile, you could test with ln_linssh = .TRUE. and see if the model proceeds further.

Hi, using

ln_linssh = .TRUE.

it proceeded briefly and aborted at time step 2 with

STOP
   stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN en
 counter in the tests
 
 kt 2 |ssh| max  1.7041E+12 at i j    86   6   MPI rank  1
 kt 2 |U|   max  3.6444E+11 at i j k  85   4 1 MPI rank  1
 kt 2 Sal   min   35.00     at i j k  43   2 1 MPI rank  1
 kt 2 Sal   max   35.02     at i j k  85  52 1 MPI rank  1
 
        ===> output of last computed fields in output.abort* files
 
 
 huge E-R-R-O-R : immediate stop

Hi, during the troubleshooting, I figured the crash happens due to the 2 layers setup. I tested by using

nn_dyn2d_dta   =  0

and as such the |U| still blows up but not the |ssh|. This led me to believe that one of the layers was causing a very high velocity.

I have checked by using the original domain with 50+ layers and it ran fine.

I am keen to further figure out if this can be fixed by tuning the make DOMAINcfg settings such as rn_theta, rn_bb or rn_hc. For now I have tried both

   rn_theta    =    20.0    !  surface control parameter (0<=theta<=20)
   rn_bb       =    1    !  stretching with SH94 s-sigma
   rn_hc       =  0.0    !  critical depth for transition to stretched coordinates

and

   rn_theta    =    0.0    !  surface control parameter (0<=theta<=20)
   rn_bb       =    0.8    !  stretching with SH94 s-sigma
   rn_hc       =  0.0    !  critical depth for transition to stretched coordinates

and both crashes.

Have you tried without bdy (ln_bdy = .false.) ? As far as I understand, you just have tidal forcing, so unless your tidal potential is on (ln_tide_pot=.true.), your model should stay at rest. That should be a good start. If, so, then we have to dig into your bdy settings.

Hi, I have been using it with ln_bdy = .false. without any difference.

Just an update that I was looking up an alternative approach in modifying the code based on the following work:

As there is an example config it will be easier to troubleshoot if we get this working. The included example domain also has 2 levels which is what we needed.

I am trying to port over the changes over to r4.0.6. For a start, I have used a reference version close to the r8814 that was used to first detect the the changes in code for the AMM7_surge project. Highlighting the difference in the files showed changes in getting rid of the vertical physics to help stabilize the single layer run and also remove the tracers.

Addition of ln_2d option, various modifications

Summary

diawri.F90
dom_oce.F90
domain.F90
domzgr.F90
dtatsd.F90
dynhpg.F90
dynkeg.F90
dynspg.F90
dynvor.F90
dynzad.F90
istate.F90
nemogcm.F90
sbc_oce.F90
step.F90

Adding tide constituents, optional

Summary

tide.F90
tide_mod.F90
tideini.F90

I seem to have implemented everything into r4.0.6, however there were quite a large change in nemogcm.F90 where by

CALL zdf_init        ! namelist read

is now

CALL zdf_phy_init    ! Vertical physics

I realised I could not just add IF (.NOT. ln_2d) in front of it as further looking into zdf_phy_init has certain functions that is not related only to vertical physics. I have gone through zdf_phy.F90 to seek out those lines that were previously in zdf_init. The code compiles but the example still crashes at the first time step, blowing up for both the ssh and U.

At the moment I believe I will need to go through the code again to seek if I can remove any other vertical physics that will cause the stability issues. If anyone has successfully done this before it will be great if you can share how to port the changes in AMM7_surge over to 4.0.6.

Hi, after incorporating the changes into the code I have reached the stage where it runs conditionally based on the number of domain levels. I have never been able to get 2 levels working as the ssh and velocity would blow up early on. 3 levels is fine for now. I was wondering if I could refine the namelist_cfg anymore to make 2 levels working. My current settings:

   jpkdta      =      3    !  number of levels      ( >= jpk )
   ln_s_sh94   = .true.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)|
   ln_s_sf12   = .false.   !  Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied
   ln_sigcrit  = .false.   !  use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
                           !  stretching coefficients for all functions
   rn_sbot_min =   10.0    !  minimum depth of s-bottom surface (>0) (m)
   rn_sbot_max = 10200.0    !  maximum depth of s-bottom surface (= ocean depth) (>0) (m)
   rn_hc       =  0.0    !  critical depth for transition to stretched coordinates
                        !!!!!!!  Envelop bathymetry
   rn_rmax     =    0.1    !  maximum cut-off r-value allowed (0<r_max<1)
                        !!!!!!!  SH94 stretching coefficients  (ln_s_sh94 = .true.)
   rn_theta    =    0.0    !  surface control parameter (0<=theta<=20) ! 3 layers min cannot use 8,20
   rn_bb       =    0.0    !  stretching with SH94 s-sigma ! 3 layers min cannot use 1
                        !!!!!!!  SF12 stretching coefficient  (ln_s_sf12 = .true.)
   rn_alpha    =    4.4    !  stretching with SF12 s-sigma
   rn_efold    =    0.0    !  efold length scale for transition to stretched coord
   rn_zs       =    1.0    !  depth of surface grid box
                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
   rn_zb_a     =    0.024  !  bathymetry scaling factor for calculating Zb
   rn_zb_b     =   -0.2    !  offset for calculating Zb
                        !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above]
   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1)

I will further test suggested settings.

Thank you.

1 Like

Dear jhc,
I want to reproduce the AMM7 experiment in a NEMO version 4.0.6, like you, did you manage to configure said version in barotropic mode? Or do you have any recommendations? If you can help me I would appreciate it very much.

Elsy

Hi jhc and Jérôme,

I’m trying to developp a barotropic configuration for the Caribbean Sea and Gulf of Mexico with NEMO4.0.6 or more.

Thanks to your comments, I reach the point where I have a correct domain_cfg.nc file with 3 layers. Using AMM7_surge code (NEMO3.6), I have run a barotropic simulation without any problem with tides, atmospheric forcings (winds and pressure), the domain_cfg file, constant temperature (=10ºC) and salinity (=35psu).

I have tried to port the changes over the r4.0.6 version as you did, modifying the 14 routines + zdfphy.F90, adding usrdef_istate.F90 and usrdef_sbc.F90. The code compiles, but the model crashes after 50 time steps (ssh>20m) in the Bahamas where there are very steep slopes (tide only run). I can see in the output.abort.nc file that the vertical velocity (vovecrtz) reaches 0.06m/s in the problematic area, so I suppose I didn’t remove all the vertical physics ?

Do you have any idea of where the problem might be or where I could modify the code ? As the ssh doesn’t blow up very early, I’m also wondering if the problem is not coming from the complex area (but it worked with v3.6…). I have also tried a simulation at rest, it works.

Thank you for your help,

Alisée

here are the main modifications I made :

grep "ln_2d" *
domain.F90:      NAMELIST/namdom/ ln_linssh, rn_isfhmin, rn_rdt, rn_atfp, ln_crs, ln_meshmask, ln_2d
domain.F90:         WRITE(numout,*) '      2D mode                                 ln_2d      = ', ln_2d
domain.F90:         IF(ln_2d) WRITE(numout,*) '   2D mode active: All tracer processes and vertical diffusion turned off'
dom_oce.F90:   LOGICAL , PUBLIC ::   ln_2d          ! (=T) run in 2D barotropic mode (no tracer processes or vertical diffusion)
domzgr.F90:      IF (ln_2d .AND. .NOT.ln_sco)   CALL ctl_stop( ' 2D mode must be used with ln_sco' )
dtatsd.F90:         IF( ln_2d ) WRITE(numout,*) '   2D ocean - ocean will be started at rest and T&S = 0'
dynhpg.F90:      IF (ln_2d) THEN
dynhpg.F90:      END IF  ! if ln_2d
istate.F90:         !IF (ln_2d) THEN
nemogcm.F90:   USE dom_oce , ONLY : ln_2d
nemogcm.F90:      IF (.NOT. ln_2d)     CALL ldf_tra_init      ! Lateral ocean tracer physics
nemogcm.F90:      IF (.NOT. ln_2d) THEN  ! No tracers in 2D mode
step.F90:      IF ( .NOT. ln_2d ) THEN
step.F90:      IF ( .NOT. ln_2d )    CALL eos    ( tsn, rhd, rhop, gdept_n(:,:,:) )  ! now in situ density for hpg computation
step.F90:      IF (.NOT. ln_2d)      CALL dyn_zdf    ( kstp )  ! vertical diffusion
step.F90:      IF (.NOT. ln_2d) THEN  ! No tracers in 2D mode
step.F90:      ENDIF ! not ln_2d
step.F90:      IF (.NOT. ln_2d)   CALL tra_nxt       ( kstp )  ! finalize (bcs) tracer fields at next time step and swap

+ remove the calls to vertical physics in *zdfphy.F90*

Hi Alisée and others,

According to me, not having the 2-layers case working (I mean jpk=2, i.e. just one active layer) is really something to solve.
I assume that from the discussion above, the domain_cfg.nc file is now ok.
Alisée, would it be possible to track the changes you made on git or elsewhere ?

Jérôme

1 Like