tide + met + WAD

[mpayopayo]

WAD makes the waterline move back and forth. It also allows to relax the hardwired condition for minimum depth (currently 10m).

1st create new branch: feature/metWAD

Steps:

1. switch on WAD. This is the namelist namwad ! Wetting and Drying (WaD) in the namelist_cfg file
2. test run with WAD switched on
3. regenerate bathymetry to relax the hardware condition of 10m minimum depth
4. run

[mpayopayo]

I tried running with the namwad namelist from the namelist_ref file

!-----------------------------------------------------------------------
&namwad        !   Wetting and drying                                   (default F)
!-----------------------------------------------------------------------
   ln_wd       = .true.   !  T/F activation of wetting and drying
   rn_wdmin1   =  0.1      !  Minimum wet depth on dried cells
   rn_wdmin2   =  0.01     !  Tolerance of min wet depth on dried cells
   rn_wdld     =  20.0     !  Land elevation below which wetting/drying is allowed
   nn_wdit     =  10       !  Max iterations for W/D limiter
/

but I get the following error:

  ===>>> : E R R O R

          ===========

 misspelled variable in namelist namwad in configuration namelist iostat =  4324

So I think is its the wrong namwad

The namwad variables in the manual (below) are different from those in the namelist_ref

!-----------------------------------------------------------------------
&namwad ! Wetting and Drying (WaD) (default: OFF)
!-----------------------------------------------------------------------
   ln_wd_il     =  .false.   ! T/F activation of iterative limiter
   ln_wd_dl     =  .true.   ! T/F activation of directional limiter
   ln_wd_dl_bc  =  .false.   ! T/F Directional limiteer Baroclinic option
   ln_wd_dl_rmp =  .true.   ! T/F Turn on directional limiter ramp
   rn_wdmin0    =  0.30      ! depth at which WaD starts
   rn_wdmin1    =  0.2       ! Minimum wet depth on dried cells
   rn_wdmin2    =  0.0001    ! Tolerance of min wet depth on dried cells
   rn_wdld      =  2.5       ! Land elevation below which WaD is allowed
   nn_wdit      =  20        ! Max iterations for WaD limiter
   rn_wd_sbcdep =  5.0       ! Depth at which to taper sbc fluxes
   rn_wd_sbcfra =  0.999     ! Fraction of SBC fluxes at taper depth (Must be <1)
/

With this namwad from the manual I get the following error:

  ===>>> : E R R O R

          ===========

 iom_varid, file: domain_cfg.nc, var: rn_wd_ref_depth not found

So domain_cfg.nc needs to be modified

NOTE rn_wdld may need to play around

[jpolton]

@endaodea Can you give @mpayopayo some pointers about the domain_cfg.nc requirements to get wetting and drying going?

This Severn estuary domain is a sort of prototype for getting an all singing and dancing REP configuration going.

[endaodea]

For W&D the reference depth is set as a single value in the domain cfg file (rn_wd_ref_depth). E.G. For AMM15 I have a clearance of 10 m which should mean any point that is above MSL by less than 10 m will be OK for the run.

double rn_wd_ref_depth ;

data:
rn_wd_ref_depth = 10 ;

If in your case all bathy points are <0 (no negative bathy) you could just add the variable directly the domain cfg file with the value of 0.

If on the other hand you do have negative depths then the cfg file has to be set up wit the reference height in mind, so that all values are less than the reference height. I'd imagine for the Severn there are many such points.

Basically if say a point has a bathy of -5 m and you used a reference height say of 10 m then the cfg value is +5m, and the model thicknesses are say for sigma coordinates 5 m/N levels where N is the number of levels.

At time step 0 if the model is set to rest, SSH=0, at these such points the SSH = Min Wet depth, say a few cm and the initial wet thicknesses are 5 cm/N Levels.

[mpayopayo]

Thank you @endaodea,

Is rn_wd_ref_depth defined on the configuration namelist? And how does rn_wd_ref_depth relate to rn_wdld?

That minimum wet depth you refer in your last sentence, is it rn_wdmin1?

[endaodea]

Hi @mpayopayo , No rn_wd_ref_depth (its a hang over from when it was a namelist variable) is now only assigned as a single float in the netcdf configuration file. WE thought we should put it there as the data there wouldn't make any sense without it. We could put it in the namelist but that opens the door to the data being not used correctly. Perhaps we should rename it so its not seen as a namelist entry anymore.

unfortunately there is an older version of the WAD scheme( ln_wd_il) that will be removed in future versions. The value of rn_wdld should have no effect when using ln_wd_dl

The minimum wet depth or critical depth is the value or rn_wdmin1 usually a few cm .

For the dl case what you usually need is:

!-----------------------------------------------------------------------
&namwad ! Wetting and Drying (WaD) (default: OFF)
!-----------------------------------------------------------------------
   ln_wd_il     =  .false.   ! T/F activation of iterative limiter
   ln_wd_dl     =  .true.   ! T/F activation of directional limiter
   ln_wd_dl_bc  =  .true.   ! T/F Directional limiteer Baroclinic option
   ln_wd_dl_rmp =  .true.   ! T/F Turn on directional limiter ramp
   rn_wdmin1    =  0.2       ! Minimum wet depth on dried cells
   rn_wd_sbcdep =  5.0       ! Depth at which to taper sbc fluxes
   rn_wd_sbcfra =  0.999     ! Fraction of SBC fluxes at taper depth (Must be <1)
/

omitted settings are for the older il scheme. (set false here)

[mpayopayo]

I get the bathy_meter variable in the domain_cfg.nc capped at zero values (no negative values) despite having negative values in the bathy_meter.nc file. I compared the f90 files I have with those of Jason who ran nemo3 for LPL Bay and while Jason's have WAD variables the ones I'm using don't have so I'm pretty sure this is the problem. I've asked Enda for the version of the code I should need.

[mpayopayo]

@jpolton @micdom The domain generated with the tools available only has positive bathymetry (>0) despite the bathymetry file having negative values, and there are no WAD variables in the domain file.

It seems that the problem is on the domain tools since there are no WAD variables in domain.f90, dom_oce.f90 ordomzgr.f90. I compared the tools we have with those that Jason for the Liverpool Bay test. The largest differences in the domain generation tool files are in domzgr.f90.

However, it seems that we are not using at the moment the version of the code that includes the latest modifications on wad @endaodea (thank you!) pointed me towards the branch that has some changes needed for wad https://forge.ipsl.jussieu.fr/nemo/browser/NEMO/branches/NERC/NEMO_4.0.4_CO9_package_tides Differences related to wad are in files dom_oce.f90, domain.f90, dtatsd.f90, and dynspg_ts.f90. Other differences are the updated version of the zdfgls.f90.

For reference the differences in the tool files used by Jason

$diff domain.f90 SRCfromJholt/domain.f90
47c47
<    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
---
>    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
49c49
<    !! Software governed by the CeCILL licence        (./LICENSE)
---
>    !! Software governed by the CeCILL licence        (NEMOGCM/NEMO_CeCILL.txt)
162c162
<          &             ppa2, ppkth2, ppacr2
---
>          &             ppa2, ppkth2, ppacr2, rn_wd_ref_depth
273a274
>          WRITE(numout,*) '      WAD Reference depth)              rn_wd_ref_depth     = ', rn_wd_ref_depth
432a434,444
> !CEOD Just force it to be this for now
>       z2d(:,:) = hbatt(:,:) ! add back on reference height to get appox dep
>                                  !this is later corrected for with specified min depth bg user for above greoid
>                                  ! WAD points
>       !where (z2d   (:,:).lte.1e-5)  z2d(:,:) = -10.0
>       where (tmask  (:,:,1).eq.0)  z2d(:,:) = 0.0
>       !CEOD CALL iom_rstput( 0, 0, inum, 'rn_wd_ref_depth'   , rn_wd_ref_depth, ktype = jp_i4  ) ! replace this later with variable
>       CALL iom_rstput( 0, 0, inum, 'rn_wd_ref_depth'   , rn_wd_ref_depth  ) ! replace this later with variable
>       CALL iom_rstput( 0, 0, inum, 'ht_wd', z2d )        !    ht_wd
> !CEOD
>
483,489d494
<       DO jj = 1,jpj
<          DO ji = 1,jpi
<             z2d (ji,jj) = SUM ( e3t_0(ji,jj, 1:mbkt(ji,jj) ) ) * ssmask(ji,jj)
<          END DO
<       END DO
<       CALL iom_rstput( 0, 0, inum, 'bathy_meter'   , z2d , ktype = jp_r4 )

$ diff dom_oce.f90 SRCfromJholt/dom_oce.f90
32a33
>    REAL(wp), PUBLIC ::   rn_wd_ref_depth !: Reference depth for Wet and Dry
291c292

$diff domzgr.f90 SRCfromJholt/domzgr.f90
55d54
<    LOGICAL  ::   ln_s_melange      ! use mix of ZPS and hybrid s-sig
62d60
<    INTEGER  ::   nn_sig_lev
1912c1910
<       INTEGER  ::   ios, mmin                      ! Local integer output status for namelist read
---
>       INTEGER  ::   ios                      ! Local integer output status for namelist read
1914,1918c1912,1913
<
<       REAL(wp) ::   rn_hc_bak   ! temporary scalars
<       REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
<
<       REAL(wp) ::   zrfact   ! temporary scalars
---
>       REAL(wp) ::   zrfact
>       !
1923c1918
<          &                rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b, nn_sig_lev, ln_s_melange
---
>          &                rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
1964,1982c1959,1978
<
<       IF(ln_s_melange) THEN
<       CALL zgr_zps          ! Partial step z-coordinate
<       IF(lwp)WRITE(numout,*) 'zgr_sco :   called zgr_zps JDHA'
<
<       ! Scale factors and depth at U-, V-, UW and VW-points
<       DO jk = 1, nn_sig_lev                        ! initialisation to z-scale factors above ln_s_sigma to remove any zps
<          gdept_0(:,:,jk) = gdept_1d(jk)
<          gdepw_0(:,:,jk) = gdepw_1d(jk)
<          e3t_0  (:,:,jk) = e3t_1d  (jk)
<          e3w_0  (:,:,jk) = e3w_1d  (jk)
<          e3f_0  (:,:,jk) = e3t_1d  (jk)
<          e3u_0(:,:,jk) = e3t_1d(jk)
<          e3u_0(:,:,jk) = e3t_1d(jk)
<          e3v_0(:,:,jk) = e3t_1d(jk)
<          e3uw_0(:,:,jk) = e3w_1d(jk)
<          e3vw_0(:,:,jk) = e3w_1d(jk)
<       END DO
<       ENDIF
---
>       ! CEOD
>       ! Add on some reference level to make nagative bathy positive to create
>       ! sensible e3t_0's
> !CEOD MAKE points on the perimeter assuming a bondary is here on shallower than
> !10m so that tide can be specified here.
> !      DO jj = 1, jpj
> !         DO ji = 1, jpi
> !            !CEODIF( mig(ji) <= 2 .or. mjg(jj) <= 2  .or. mig(ji) >= jpiglo-2 .or.  mjg(jj) >=jpjglo-2) then
> !            IF( mig(ji) <= 10   .or. mjg(jj) <=10  .or.  mjg(ji) >=jpiglo-10 .or. mjg(jj) >=jpjglo-10   ) then
> !               IF( bathy(ji,jj) > 0.0)   then
> !                  bathy(ji,jj) = max(5.0,bathy(ji,jj)) !
> !               ENDIF
> !               IF( bathy(ji,jj) <= 0.0)   then
> !                  bathy(ji,jj) =  -rn_wd_ref_depth*2._wp ! Make sure it is land
> !               ENDIF
> !            ENDIF
> !         END DO
> !      END DO
>       bathy(:,:) = bathy(:,:) + rn_wd_ref_depth !  Arbitrary reference level, This needs
>                                       !  to be accounted for in istate in the actual model and bdys etc.
2003,2009d1998
<       IF(ln_s_melange) THEN
<       DO jj = 1, jpj
<          DO ji = 1, jpi
<            zenv(ji,jj) = MIN( bathy(ji,jj), gdepw_1d(nn_sig_lev + 1) )
<          ENDDO
<       ENDDO
<       ENDIF
2110d2098
<       !hbatz(:,:) = zenv(:,:)
2120,2122d2107
<       zrmax = MAXVAL( hbatt(:,:) )
<       IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
<       IF(lwp) WRITE(numout,*) 'MAX HBATT ', zrmax
2247,2260d2231
<       IF(ln_s_melange) THEN
<       IF( lwp )  WRITE(numout,*) 'JDHA - ADJUSTING MBATHY'
<       IF( lwp )  WRITE(numout,*) nn_sig_lev, gdepw_1d(nn_sig_lev + 1)
<       DO jj = 1, jpj
<          DO ji = 1, jpi
<             IF( bathy(ji,jj) < gdepw_1d(nn_sig_lev + 1) ) THEN ! should this be hbatt or bathy
<                DO jk = 1, nn_sig_lev
<                   IF( scobot(ji,jj) >= gdept_n(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
<                END DO
<             ENDIF
<          END DO
<       END DO
<       ELSE
<       IF( lwp )  WRITE(numout,*) 'JDHA - ADJUSTING MBATHY SHOULD NOT BE HERE'
2268,2269d2238
<       ENDIF
<
2326,2327c2295
<                mmin = min(mbathy(ji,jj),nn_sig_lev)
<                DO jk = 1,mmin
---
>                DO jk = 1, mbathy(ji,jj)
2329c2297
<                  IF (e3w_0(ji,jj,jk) <= 0._wp .OR. e3t_0(ji,jj,jk) <= 0._wp ) THEN
---
>                  IF (e3w_n(ji,jj,jk) <= 0._wp .OR. e3t_n(ji,jj,jk) <= 0._wp ) THEN
2332,2333c2300,2301
<                     WRITE(numout,*) 'e3w',e3w_0(ji,jj,:)
<                     WRITE(numout,*) 'e3t',e3t_0(ji,jj,:)
---
>                     WRITE(numout,*) 'e3w',e3w_n(ji,jj,:)
>                     WRITE(numout,*) 'e3t',e3t_n(ji,jj,:)
2337c2305
<                  IF( gdepw_0(ji,jj,jk) < 0._wp .OR. gdept_0(ji,jj,jk) < 0._wp ) THEN
---
>                  IF( gdepw_n(ji,jj,jk) < 0._wp .OR. gdept_n(ji,jj,jk) < 0._wp ) THEN
2340,2341c2308,2309
<                     WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:)
<                     WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
---
>                     WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:)
>                     WRITE(numout,*) 'gdept',gdept_n(ji,jj,:)
2345c2313
<                  IF( gdepw_0(ji,jj,jk) > hbatt(ji,jj) ) THEN
---
>                  IF( gdepw_n(ji,jj,jk) > hbatt(ji,jj) ) THEN
2348c2316
<                     WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,jk-1:jk+1), bathy(ji,jj),  hbatt(ji,jj), jk, mbathy(ji,jj)
---
>                     WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:)
2353,2354c2321
<          !     DO jk = 1, mbathy(ji,jj)-1
<                DO jk = 1, mmin-1
---
>                DO jk = 1, mbathy(ji,jj)-1
2356c2323
<                 IF( gdept_0(ji,jj,jk) > hbatt(ji,jj) ) THEN
---
>                 IF( gdept_n(ji,jj,jk) > hbatt(ji,jj) ) THEN
2359c2326
<                     WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
---
>                     WRITE(numout,*) 'gdept',gdept_n(ji,jj,:)
2390d2356
<       REAL(wp) ::   rn_hc_bak   ! temporary scalars
2408,2416c2374,2375
<                   IF(ln_s_melange) THEN
<                   z_gsigw3(ji,jj,jk) = gdepw_1d(jk)/gdepw_1d(nn_sig_lev + 1)
<                   z_gsigt3(ji,jj,jk) = gdept_1d(jk)/gdepw_1d(nn_sig_lev + 1)
<                   ELSE
<                   !JDHA z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
<                   !JDHA z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
<                   z_gsigw3(ji,jj,jk) = -fssig1_jdha( REAL(jk,wp)-0.5_wp, rn_bb, hbatt(ji,jj) ) !JDHA
<                   z_gsigt3(ji,jj,jk) = -fssig1_jdha( REAL(jk,wp)       , rn_bb, hbatt(ji,jj) ) !JDHA
<                   ENDIF
---
>                   z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
>                   z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
2420,2423d2378
<                   IF(ln_s_melange) THEN
<                   z_gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(nn_sig_lev,wp)
<                   z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(nn_sig_lev,wp)
<                   ELSE
2426,2427c2381
<                   ENDIF
<                END DO
---
>                   END DO
2443,2474d2396
<
<             IF(ln_s_melange) THEN
<
<             DO jk = 1, nn_sig_lev+1
<                IF( bathy(ji,jj) < gdepw_1d(nn_sig_lev + 1) ) THEN ! should this be bathy or hbatt?
<
<                   zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(nn_sig_lev,wp)
<                   zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(nn_sig_lev,wp)
< !                 zcoeft = ( REAL(MIN(jk,nn_sig_lev),wp) - 0.5_wp ) / REAL(nn_sig_lev-1,wp)
< !                 zcoefw = ( REAL(MIN(jk,nn_sig_lev),wp) - 1.0_wp ) / REAL(nn_sig_lev-1,wp)
<
<                   rn_hc_bak = rn_hc
<                   rn_hc = MIN( MAX ( &
<                 &            (hbatt(ji,jj)-gdepw_1d(nn_sig_lev + 1)) / (1._wp - (gdepw_1d(nn_sig_lev + 1)/rn_hc)) &
<                 &                   ,0._wp) ,rn_hc)
<
<                   gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
<                   gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
<                   gde3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
<
<                   rn_hc = rn_hc_bak
<
<                   IF( gdepw_0(ji,jj,jk) < 0._wp ) THEN
<                      WRITE(*,*) 'zgr_sco :   gdepw  at point (i,j,k)= ', ji, jj, jk, (z_gsigw3(ji,jj,jk)*10000._wp-zcoefw*10000._wp)
<                   ENDIF
<
<                ENDIF
<
<             END DO
<
<             ELSE
<
2482,2483d2403
<
<             ENDIF !melange
2490,2563d2409
<
<             IF(ln_s_melange) THEN
<
<             ! extended for Wetting/Drying case
<             ztmpu  = hbatt(ji,jj)+hbatt(ji+1,jj)
<             ztmpv  = hbatt(ji,jj)+hbatt(ji,jj+1)
<             ztmpf  = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
<             ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
<             ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
<             ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
<                    & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
<
<             IF( bathy(ji,jj) < gdepw_1d(nn_sig_lev + 1) ) THEN ! should this be bathy or hbatt?
<             DO jk = 1, nn_sig_lev+1      ! scale factors should be the same in both zps and sco when H > Hcrit??
< !           DO jk = 1, jpk
<
<                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
<                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
<                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
<                   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
<                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     &
<                   &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   &
<                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
<                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
<                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
<                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
<                   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
<                !
<
<                rn_hc_bak = rn_hc
<                rn_hc = MIN( MAX( &
<               &            (hbatt(ji,jj)-gdepw_1d(nn_sig_lev + 1)) / (1._wp - (gdepw_1d(nn_sig_lev + 1)/rn_hc)) &
<               &                   ,0._wp) ,rn_hc)
< !              e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) )
< !              e3w_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) )
<                e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) )
<                e3w_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) )
<
<
<                rn_hc = MIN( MAX( &
<               &            (hbatu(ji,jj)-gdepw_1d(nn_sig_lev + 1)) / (1._wp - (gdepw_1d(nn_sig_lev + 1)/rn_hc_bak)) &
<               &                   ,0._wp) ,rn_hc_bak)
<                e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) )
<                e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) )
<                e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) )
<                e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) )
<
<
<                rn_hc = MIN( MAX( &
<               &            (hbatv(ji,jj)-gdepw_1d(nn_sig_lev + 1)) / (1._wp - (gdepw_1d(nn_sig_lev + 1)/rn_hc_bak)) &
<               &                   ,0._wp) ,rn_hc_bak)
< !              e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) )
< !              e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) )
<                e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) )
<                e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) )
<
<
<
<                rn_hc = MIN( MAX( &
<               &            (hbatf(ji,jj)-gdepw_1d(nn_sig_lev + 1)) / (1._wp - (gdepw_1d(nn_sig_lev + 1)/rn_hc_bak)) &
<               &                   ,0._wp), rn_hc_bak)
< !              e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) )
<                e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) )
<         !
<
<                rn_hc = rn_hc_bak
<
<
<             END DO
<
<             ENDIF
<
<             ELSE
<
2598,2600d2443
<
<             ENDIF ! melange
<
2871,2902d2713
<  FUNCTION fssig1_jdha( pk1, pbb, phbatt ) RESULT( pf1 )
<       !!----------------------------------------------------------------------
<       !!                 ***  ROUTINE fssig1 ***
<       !!
<       !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
<       !!
<       !! ** Method  :   the function provides the non-dimensional position of
<       !!                T and W (i.e. between 0 and 1)
<       !!                T-points at integer values (between 1 and jpk)
<       !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
<       !!----------------------------------------------------------------------
<       REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
<       REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
<       REAL(wp), INTENT(in) ::   phbatt!JDHA
<       REAL(wp)             ::   pf1   ! sigma value
<       REAL(wp)             ::   pbb_mod!JDHA
<       !!----------------------------------------------------------------------
<       !
<       pbb_mod = pbb * ( 1._wp + ( 1._wp - EXP( ( phbatt / rn_sbot_max ) ) ) / EXP(1._wp) ) !JDHA
<       IF ( rn_theta == 0 ) then      ! uniform sigma
<          pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
<       ELSE                        ! stretched sigma
<
<          !JDHApf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
<          !JDHA   &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
<          !JDHA   &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
<          pf1 =   ( 1._wp - pbb_mod ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
<             &  + pbb_mod * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
<             &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
<       ENDIF
<       !
<    END FUNCTION fssig1_jdha
2917d2727
<
2919d2728
<
2922d2730
<
2926d2733
<
2928,2930c2735,2736
<            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
<           &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
<
---
>             &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
>             &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )

[endaodea]

Just a comment to say most of the package branch changes concern a NEMO bug (in my opinion) with regards to how depths at U/V points are derived. In the original code it uses the initial hu_0 hv_0 to infer the now u/v depths using ssh_n. But that does not actually work.
The top panel is at the initial state of rest when all is well with ssh=0 The second panel is at some future state where the SSH has dropped.

The fix in the package branch just updates e3un and e3vn explicitly instead of relying on hu_0 and hv_0