Deriving scale factors if ln_isfcav = true

[jamesharle]

Context

When employing ln_isfcav = true and generating a domain_cfg.nc file using the DOMAINcfg tool the vertical scale factors are altered in the shallowest wet cell as follows:

in domzgr.f90 (ln1828):

e3w_0  (ji,jj,ik  ) = e3t_0  (ji,jj,ik)

this is to accommodate the use of partial cells under ice sheets. The consequence of this is that an incorrect gdept_0 is generated using e3_to_depth (in zgr_read) in the active (wet) cells wherever the ice sheet is present. 

Dynamically this is not an issue as far as I can tell as gdept_n is corrected for this fact in dom_vvl_init, whilst still perpetuating the e3w_n based on e3t_0(ji,jj,ik) - (also in iscplrst.F90).

The consequence of this approach is that the thickness factors e3w_0 and e3w_n, and depth coordinate gdept_0 are incorrect in domain_cfg.nc, meshmask.nc and any output files. The only way the user can back out the true depth at t-points is to update diawri.F90 to output gdept_n or use the ice draft depth with 0.5*e3w_0 at the first wet cell. I thought it would be cleaner and clearer if the treatment of the scale thickness of the first wet cell under the ice sheet was treated locally within the dynamics code rather than globally (corrupting the vertical coordinates).

As the alteration to e3w at the first wet cell under the ice sheet was not immediately apparent, this proved interesting when deriving domain_cfg.nc files for both sh94 and hybrid z-s coordinates. Below are my changes to domzgr.f90 (in tools/DOMAINcfg) to remove this adaption and to the other code in src/OCE that has to be updated accordingly.

I have tested both the old code and the changes below on an ISOMIP-like set up. Both yielding identical run.stat files.

I’ve raised this ticket on the suggestion of Pierre, who says he’ll be revisiting the ice cavity coordinates in his 2019WP.

Proposal

Firstly to retain the full e3w thicknesses in tools/DOMAINcfg/src

domzgr.f90

@@ -1825 +1825 @@
                   e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 
                ENDIF 
             !       ... on ik / ik-1 
-               e3w_0  (ji,jj,ik  ) = e3t_0  (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 
                e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
-! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 
-               gdept_0(ji,jj,ik-1) = gdept_1d(ik-1)
+               e3w_0  (ji,jj,ik) = e3w_0  (ji,jj,ik) + ( gdept_0(ji,jj,ik) - gdept_1d(ik) )
             ENDIF 
          END DO 
       END DO 

In domvvl we therefore have to retain the conventional way of defining gdept_n without an adjustment for the uppermost wet cell under the ice sheet (this will be handled locally in dyn_hpg_isf).

domvvl.F90

@@ -116 +116 @@
       !!----------------------------------------------------------------------
       INTEGER ::   ji, jj, jk
       INTEGER ::   ii0, ii1, ij0, ij1
-      REAL(wp)::   zcoef
       !!----------------------------------------------------------------------
       !
       IF(lwp) WRITE(numout,*)
@@ -161 +160 @@
       DO jk = 2, jpk                               ! vertical sum
          DO jj = 1,jpj
             DO ji = 1,jpi
-               !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
-               !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
-               !                             ! 0.5 where jk = mikt     
-!!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ??
-               zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) )
                gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
-               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
-                  &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk)) 
+               gdept_n(ji,jj,jk) = gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk) 
                gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
                gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
-               gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  &
-                  &                + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk)) 
+               gdept_b(ji,jj,jk) = gdept_b(ji,jj,jk-1) + e3w_b(ji,jj,jk) 
             END DO
          END DO
       END DO
@@ -589 +581 @@
       INTEGER, INTENT( in ) ::   kt   ! time step
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zcoef        ! local scalar
       !!----------------------------------------------------------------------
       !
       IF( ln_linssh )   RETURN      ! No calculation in linear free surface
@@ -645 +636 @@
       DO jk = 2, jpk
          DO jj = 1,jpj
             DO ji = 1,jpi
-              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
-                                                                 ! 1 for jk = mikt
-               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
                gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
-               gdept_n(ji,jj,jk) =    zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk) )  &
-                   &             + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk) ) 
+               gdept_n(ji,jj,jk) = gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)
                gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
             END DO
          END DO

The correction term is also updated in iscplrst.F90:

iscplrst.F90

@@ -222 +222 @@
                END DO
                IF (mikt(ji,jj) > 1) THEN
                   jk = mikt(ji,jj)
-                  gdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * e3w_n(ji,jj,jk)
+                  gdept_n(ji,jj,jk) = gdept_0(ji,jj,jk-1) + e3w_n(ji,jj,jk)
                   gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
                   gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn   (ji,jj)
                END IF

In dyn_hpg_isf the definition of e3w in the first wet cell under the ice sheet is added:

dynhpg.F90

@@ -581 +581 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk ) ::  zhpi, zhpj
       REAL(wp), DIMENSION(jpi,jpj,jpts) ::  zts_top
       REAL(wp), DIMENSION(jpi,jpj)      ::  zrhdtop_oce
+      REAL(wp), DIMENSION(jpi,jpj)      ::  ze3w_top
       !!----------------------------------------------------------------------
       !
       zcoef0 = - grav * 0.5_wp   ! Local constant initialization
@@ -592 +593 @@
 
       ! compute rhd at the ice/oce interface (ocean side)
       ! usefull to reduce residual current in the test case ISOMIP with no melting
+      ! compute the half wet cell thickness just under the ice sheet. NB this is not
+      ! 0.5 * e3w_n(:,:,ikt) due to the partial step in the ice sheet
       DO ji = 1, jpi
         DO jj = 1, jpj
           ikt = mikt(ji,jj)
-          zts_top(ji,jj,1) = tsn(ji,jj,ikt,1)
-          zts_top(ji,jj,2) = tsn(ji,jj,ikt,2)
+          zts_top (ji,jj,1) = tsn(ji,jj,ikt,1)
+          zts_top (ji,jj,2) = tsn(ji,jj,ikt,2)
+          ze3w_top(ji,jj  ) = gdept_n(ji,jj,ikt) - gdepw_n(ji,jj,ikt) 
         END DO
       END DO
       CALL eos( zts_top, risfdep, zrhdtop_oce )
 
-!==================================================================================     
 !===== Compute surface value ===================================================== 
 !==================================================================================
       DO jj = 2, jpjm1
@@ -611 +614 @@
             iktp1j = mikt(ji,jj+1)
             ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure
             ! we assume ISF is in isostatic equilibrium
-            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w_n(ji+1,jj,iktp1i)                                    &
+            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( ze3w_top(ji+1,jj)                                                 &
                &                                    * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) )   &
-               &                                  - 0.5_wp * e3w_n(ji,jj,ikt)                                         &
+               &                                  - ze3w_top(ji,jj)                                                   &
                &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          &
                &                                  + ( riceload(ji+1,jj) - riceload(ji,jj))                            ) 
-            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w_n(ji,jj+1,iktp1j)                                    &
+            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( ze3w_top(ji,jj+1)                                                 &
                &                                    * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) )   &
-               &                                  - 0.5_wp * e3w_n(ji,jj,ikt)                                         & 
+               &                                  - ze3w_top(ji,jj)                                                   &
                &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          &
                &                                  + ( riceload(ji,jj+1) - riceload(ji,jj))                            ) 
             ! s-coordinate pressure gradient correction (=0 if z coordinate)

[mathiot]

For the time being, the computation of the deepest e3w_0 has been changed in the ice shelf in r10490. According to this change a modification in the computation of the isf load was required (r10491).
Now the e3_to_depth give correct answer and e3 and depth in the ice shelf point ('land point') in domain_cfg and mesh_mask are correct. ISOMIP run using user_def or domain_cfg give the same answer and after 30 years results are as expected.

With James, we will revisit this in the WP2019 as part of his action on sigma under the iceshelf and my action on ice shelf coupling, geometry and cleaning/simplification.

[mathiot]

Fixed in WP2019