Opened 18 months ago

Last modified 2 months ago

#2177 new Defect

Deriving scale factors if ln_isfcav = true

Reported by: jamesharle Owned by: systeam
Priority: low Milestone:
Component: tools Version: trunk
Severity: minor Keywords: DOMAINcfg, TOOLS, domzgr, ISF
Cc:

Description

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

     
    18251825                  e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
    18261826               ENDIF  
    18271827            !       ... on ik / ik-1  
    1828                e3w_0  (ji,jj,ik  ) = e3t_0  (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
    18291828               e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
    1830 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code  
    1831                gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 
     1829               e3w_0  (ji,jj,ik) = e3w_0  (ji,jj,ik) + ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) 
    18321830            ENDIF  
    18331831         END DO  
    18341832      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

     
    116116      !!---------------------------------------------------------------------- 
    117117      INTEGER ::   ji, jj, jk 
    118118      INTEGER ::   ii0, ii1, ij0, ij1 
    119       REAL(wp)::   zcoef 
    120119      !!---------------------------------------------------------------------- 
    121120      ! 
    122121      IF(lwp) WRITE(numout,*) 
     
    161160      DO jk = 2, jpk                               ! vertical sum 
    162161         DO jj = 1,jpj 
    163162            DO ji = 1,jpi 
    164                !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
    165                !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
    166                !                             ! 0.5 where jk = mikt      
    167 !!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ?? 
    168                zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 
    169163               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 
    170                gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  & 
    171                   &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk))  
     164               gdept_n(ji,jj,jk) = gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)  
    172165               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 
    173166               gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1) 
    174                gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  & 
    175                   &                + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk))  
     167               gdept_b(ji,jj,jk) = gdept_b(ji,jj,jk-1) + e3w_b(ji,jj,jk)  
    176168            END DO 
    177169         END DO 
    178170      END DO 
     
    589581      INTEGER, INTENT( in ) ::   kt   ! time step 
    590582      ! 
    591583      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    592       REAL(wp) ::   zcoef        ! local scalar 
    593584      !!---------------------------------------------------------------------- 
    594585      ! 
    595586      IF( ln_linssh )   RETURN      ! No calculation in linear free surface 
     
    645636      DO jk = 2, jpk 
    646637         DO jj = 1,jpj 
    647638            DO ji = 1,jpi 
    648               !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
    649                                                                  ! 1 for jk = mikt 
    650                zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
    651639               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 
    652                gdept_n(ji,jj,jk) =    zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk) )  & 
    653                    &             + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk) )  
     640               gdept_n(ji,jj,jk) = gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk) 
    654641               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 
    655642            END DO 
    656643         END DO 

The correction term is also updated in iscplrst.F90:

  • iscplrst.F90

     
    222222               END DO 
    223223               IF (mikt(ji,jj) > 1) THEN 
    224224                  jk = mikt(ji,jj) 
    225                   gdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * e3w_n(ji,jj,jk) 
     225                  gdept_n(ji,jj,jk) = gdept_0(ji,jj,jk-1) + e3w_n(ji,jj,jk) 
    226226                  gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    227227                  gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn   (ji,jj) 
    228228               END IF 

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

  • dynhpg.F90

     
    581581      REAL(wp), DIMENSION(jpi,jpj,jpk ) ::  zhpi, zhpj 
    582582      REAL(wp), DIMENSION(jpi,jpj,jpts) ::  zts_top 
    583583      REAL(wp), DIMENSION(jpi,jpj)      ::  zrhdtop_oce 
     584      REAL(wp), DIMENSION(jpi,jpj)      ::  ze3w_top 
    584585      !!---------------------------------------------------------------------- 
    585586      ! 
    586587      zcoef0 = - grav * 0.5_wp   ! Local constant initialization 
     
    592593 
    593594      ! compute rhd at the ice/oce interface (ocean side) 
    594595      ! usefull to reduce residual current in the test case ISOMIP with no melting 
     596      ! compute the half wet cell thickness just under the ice sheet. NB this is not 
     597      ! 0.5 * e3w_n(:,:,ikt) due to the partial step in the ice sheet 
    595598      DO ji = 1, jpi 
    596599        DO jj = 1, jpj 
    597600          ikt = mikt(ji,jj) 
    598           zts_top(ji,jj,1) = tsn(ji,jj,ikt,1) 
    599           zts_top(ji,jj,2) = tsn(ji,jj,ikt,2) 
     601          zts_top (ji,jj,1) = tsn(ji,jj,ikt,1) 
     602          zts_top (ji,jj,2) = tsn(ji,jj,ikt,2) 
     603          ze3w_top(ji,jj  ) = gdept_n(ji,jj,ikt) - gdepw_n(ji,jj,ikt)  
    600604        END DO 
    601605      END DO 
    602606      CALL eos( zts_top, risfdep, zrhdtop_oce ) 
    603607 
    604 !==================================================================================      
    605608!===== Compute surface value =====================================================  
    606609!================================================================================== 
    607610      DO jj = 2, jpjm1 
     
    611614            iktp1j = mikt(ji,jj+1) 
    612615            ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 
    613616            ! we assume ISF is in isostatic equilibrium 
    614             zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w_n(ji+1,jj,iktp1i)                                    & 
     617            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( ze3w_top(ji+1,jj)                                                 & 
    615618               &                                    * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) )   & 
    616                &                                  - 0.5_wp * e3w_n(ji,jj,ikt)                                         & 
     619               &                                  - ze3w_top(ji,jj)                                                   & 
    617620               &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          & 
    618621               &                                  + ( riceload(ji+1,jj) - riceload(ji,jj))                            )  
    619             zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w_n(ji,jj+1,iktp1j)                                    & 
     622            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( ze3w_top(ji,jj+1)                                                 & 
    620623               &                                    * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) )   & 
    621                &                                  - 0.5_wp * e3w_n(ji,jj,ikt)                                         &  
     624               &                                  - ze3w_top(ji,jj)                                                   & 
    622625               &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          & 
    623626               &                                  + ( riceload(ji,jj+1) - riceload(ji,jj))                            )  
    624627            ! s-coordinate pressure gradient correction (=0 if z coordinate) 

Commit History (0)

(No commits)

Change History (2)

comment:1 Changed 17 months ago by 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.

comment:2 Changed 2 months ago by clevy

  • Type changed from Enhancement to Defect
Note: See TracTickets for help on using tickets.