New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 11308 – NEMO

Changeset 11308


Ignore:
Timestamp:
2019-07-19T15:04:36+02:00 (5 years ago)
Author:
mathiot
Message:

ENHANCE-02_ISF_domcfg: create domisf to isolate isf related work and add new isf geometry computation used by UKESM (ticket #2142)

Location:
NEMO/branches/2019/ENHANCE-02_ISF_domcfg
Files:
1 added
6 edited
1 copied

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/ENHANCE-02_ISF_domcfg/namelist_cfg

    r10727 r11308  
    11!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    2 !! NEMO/OPA  Configuration namelist : used to overwrite defaults values defined in SHARED/namelist_ref 
     2!! NEMO/OCE :   NEMO/OPA  Configuration namelist : used to overwrite defaults values defined in SHARED/namelist_ref 
    33!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    4 ! 
     4!! NEMO/OCE  :  1 - Domain & run manager (namrun, namcfg, namdom, namzgr, namzgr_isf, namzgr_sco ) 
     5!!              8 - diagnostics      (namnc4) 
     6!!             10 - miscellaneous    (nammpp, namctl) 
     7!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    58!----------------------------------------------------------------------- 
    69&namrun        !   parameters of the run 
    710!----------------------------------------------------------------------- 
    8    nn_no       =       0   !  job number (no more used...) 
    9    cn_exp      =  "domaincfg"  !  experience name 
    10    nn_it000    =       1   !  first time step 
    11    nn_itend    =      75   !  last  time step (std 5475) 
    1211/ 
    1312!----------------------------------------------------------------------- 
    1413&namcfg        !   parameters of the configuration 
    1514!----------------------------------------------------------------------- 
    16    ! 
    17    ln_e3_dep   = .true.    ! =T : e3=dk[depth] in discret sens.  
    18    !                       !      ===>>> will become the only possibility in v4.0 
    19    !                       ! =F : e3 analytical derivative of depth function 
    20    !                       !      only there for backward compatibility test with v3.6 
    21    !                       !       
    22    cp_cfg      =  "orca"   !  name of the configuration 
    23    jp_cfg      =       2   !  resolution of the configuration 
    24    jpidta      =     182   !  1st lateral dimension ( >= jpi ) 
    25    jpjdta      =     149   !  2nd    "         "    ( >= jpj ) 
    26    jpkdta      =      31   !  number of levels      ( >= jpk ) 
    27    jpiglo      =     182   !  1st dimension of global domain --> i =jpidta 
    28    jpjglo      =     149   !  2nd    -                  -    --> j  =jpjdta 
    29    jpizoom     =       1   !  left bottom (i,j) indices of the zoom 
    30    jpjzoom     =       1   !  in data domain indices 
    31    jperio      =       4   !  lateral cond. type (between 0 and 6) 
    3215/ 
    33 !----------------------------------------------------------------------- 
    34 &namzgr        !   vertical coordinate 
    35 !----------------------------------------------------------------------- 
    36    ln_zps      = .true.    !  z-coordinate - partial steps 
    37    ln_linssh   = .true.    !  linear free surface 
    38 / 
    39 &namzgr_sco 
    40 / 
    41 &namlbc 
    42 / 
    43 &namagrif 
    44 / 
    45 &nambdy 
    46 / 
    47 &nam_vvl 
    48 / 
    49  
    5016!----------------------------------------------------------------------- 
    5117&namdom        !   space and time domain (bathymetry, mesh, timestep) 
    5218!----------------------------------------------------------------------- 
    53    jphgr_msh   =       0               !  type of horizontal mesh 
    54    ppglam0     =  999999.0             !  longitude of first raw and column T-point (jphgr_msh = 1) 
    55    ppgphi0     =  999999.0             ! latitude  of first raw and column T-point (jphgr_msh = 1) 
    56    ppe1_deg    =  999999.0             !  zonal      grid-spacing (degrees) 
    57    ppe2_deg    =  999999.0             !  meridional grid-spacing (degrees) 
    58    ppe1_m      =  999999.0             !  zonal      grid-spacing (degrees) 
    59    ppe2_m      =  999999.0             !  meridional grid-spacing (degrees) 
    60    ppsur       =   -4762.96143546300   !  ORCA r4, r2 and r05 coefficients 
    61    ppa0        =     255.58049070440   ! (default coefficients) 
    62    ppa1        =     245.58132232490   ! 
    63    ppkth       =      21.43336197938   ! 
    64    ppacr       =       3.0             ! 
    65    ppdzmin     =  999999.              !  Minimum vertical spacing 
    66    pphmax      =  999999.              !  Maximum depth 
    67    ldbletanh   =  .FALSE.              !  Use/do not use double tanf function for vertical coordinates 
    68    ppa2        =  999999.              !  Double tanh function parameters 
    69    ppkth2      =  999999.              ! 
    70    ppacr2      =  999999.              ! 
    7119/ 
    7220!----------------------------------------------------------------------- 
    73 &nammpp        !   Massively Parallel Processing                        ("key_mpp_mpi) 
     21&namzgr        !   vertical coordinate                                  (default: NO selection) 
    7422!----------------------------------------------------------------------- 
    7523/ 
    7624!----------------------------------------------------------------------- 
    77 &namctl        !   Control prints & Benchmark 
     25&namzgr_isf    !   isf cavity geometry definition 
    7826!----------------------------------------------------------------------- 
    7927/ 
     28!----------------------------------------------------------------------- 
     29&namzgr_sco    !   s-coordinate or hybrid z-s-coordinate                (default F) 
     30!----------------------------------------------------------------------- 
     31/ 
     32!----------------------------------------------------------------------- 
     33&namlbc        !   lateral momentum boundary condition                  (default: NO selection) 
     34!----------------------------------------------------------------------- 
     35/ 
     36!----------------------------------------------------------------------- 
     37&namagrif      !  AGRIF zoom                                            ("key_agrif") 
     38!----------------------------------------------------------------------- 
     39/ 
     40!----------------------------------------------------------------------- 
     41&nambdy        !  unstructured open boundaries                          (default: OFF) 
     42!----------------------------------------------------------------------- 
     43/ 
     44!----------------------------------------------------------------------- 
     45&namnc4        !   netcdf4 chunking and compression settings            ("key_netcdf4") 
     46!----------------------------------------------------------------------- 
     47/ 
     48!----------------------------------------------------------------------- 
     49&nammpp        !   Massively Parallel Processing                        ("key_mpp_mpi") 
     50!----------------------------------------------------------------------- 
     51/ 
     52!----------------------------------------------------------------------- 
     53&namctl        !   Control prints                                       (default: OFF) 
     54!----------------------------------------------------------------------- 
     55/ 

  • NEMO/branches/2019/ENHANCE-02_ISF_domcfg/namelist_ref

    r11201 r11308  
    4444   nn_msh      =    0      !  create (=1) a mesh file or not (=0) 
    4545   rn_hmin     =   -3.     !  min depth of the ocean (>0) or min number of ocean level (<0) 
    46    rn_isfhmin  =    1.00   !  treshold (m) to discriminate grounding ice to floating ice 
    4746   rn_e3zps_min=   20.     !  partial step thickness is set larger than the minimum of 
    4847   rn_e3zps_rat=    0.1    !  rn_e3zps_min and rn_e3zps_rat*e3t, with 0<rn_e3zps_rat<1 
    … …  
    8786   cp_cfz      = "no zoom" !  name of the zoom of configuration 
    8887   jp_cfg      =      0    !  resolution of the configuration 
    89    jpidta      =     10    !  1st lateral dimension ( >= jpi ) 
    90    jpjdta      =     12    !  2nd    "         "    ( >= jpj ) 
    9188   jpkdta      =     31    !  number of levels      ( >= jpk ) 
    9289   jpiglo      =     10    !  1st dimension of global domain --> i =jpidta 
    … …  
    10097   ln_use_jattr = .false.  !  use (T) the file attribute: open_ocean_jstart, if present 
    10198                           !  in netcdf input files, as the start j-row for reading 
    102    ln_domclo = .true.      ! computation of closed sea masks (see namclo) 
     99   ln_domclo = .false.     ! computation of closed sea masks (see namclo) 
    103100/ 
    104101!----------------------------------------------------------------------- 
    … …  
    108105   ln_zps      = .false.   !  z-coordinate - partial steps 
    109106   ln_sco      = .false.   !  s- or hybrid z-s-coordinate 
    110    ln_isfcav   = .false.   !  ice shelf cavity 
     107   ln_isfcav   = .false.   !  ice shelf cavity             (T: see namzgr_isf) 
    111108   ln_linssh   = .false.   !  linear free surface 
    112109/ 
    113110!----------------------------------------------------------------------- 
    114 &namzgr_sco    !   s-coordinate or hybrid z-s-coordinate                (default F) 
     111&namzgr_isf    !   isf cavity geometry definition                       (default: OFF) 
     112!----------------------------------------------------------------------- 
     113   rn_isfdep_min    = 10.         ! minimum isf draft tickness (if lower, isf draft set to this value) 
     114   rn_glhw_min      = 1.e-3       ! minimum water column thickness to define the grounding line 
     115   rn_isfhw_min     = 10          ! minimum water column thickness in the cavity once the grounding line defined. 
     116   ln_isfchannel    = .false.     ! remove channel (based on 2d mask build from isfdraft-bathy) 
     117   ln_isfconnect    = .false.     ! force connection under the ice shelf (based on 2d mask build from isfdraft-bathy) 
     118      nn_kisfmax       = 999         ! limiter in level on the previous condition. (if change larger than this number, get back to value before we enforce the connection) 
     119      rn_zisfmax       = 7000.       ! limiter in m     on the previous condition. (if change larger than this number, get back to value before we enforce the connection) 
     120   ln_isfcheminey   = .false.     ! close cheminey 
     121   ln_isfsubgl      = .false.     ! remove subglacial lake created by the remapping process 
     122      rn_isfsubgllon   =    0.0      !  longitude of the seed to determine the open ocean 
     123      rn_isfsubgllat   =    0.0      !  latitude  of the seed to determine the open ocean 
     124/ 
     125!----------------------------------------------------------------------- 
     126&namzgr_sco    !   s-coordinate or hybrid z-s-coordinate                (default: OFF) 
    115127!----------------------------------------------------------------------- 
    116128   ln_s_sh94   = .false.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)| 
    … …  
    137149/ 
    138150!----------------------------------------------------------------------- 
    139 &namclo ! (closed sea : need ln_domclo = .true. in namcfg) 
     151&namclo ! (closed sea : need ln_domclo = .true. in namcfg)              (default: OFF) 
    140152!----------------------------------------------------------------------- 
    141153   rn_lon_opnsea = -2.0     ! longitude seed of open ocean 
    … …  
    231243&namctl        !   Control prints                                       (default: OFF) 
    232244!----------------------------------------------------------------------- 
    233    ln_ctl = .FALSE.                 ! Toggle all report printing on/off (T/F); Ignored if sn_cfctl%l_config is T 
     245   ln_ctl = .TRUE.                 ! Toggle all report printing on/off (T/F); Ignored if sn_cfctl%l_config is T 
    234246     sn_cfctl%l_config = .TRUE.     ! IF .true. then control which reports are written with the following 
    235247       sn_cfctl%l_runstat = .FALSE. ! switches and which areas produce reports with the proc integer settings. 

  • NEMO/branches/2019/ENHANCE-02_ISF_domcfg/src/dom_oce.F90

    r11201 r11308  
    3636   REAL(wp), PUBLIC ::   rn_e3zps_rat    !: minimum thickness ration for partial steps 
    3737   INTEGER , PUBLIC ::   nn_msh          !: = 1 create a mesh-mask file 
     38 
     39   INTEGER , PUBLIC ::   ntopo           !: = 0/1 ,compute/read the bathymetry file 
     40   INTEGER, PUBLIC  ::   nperio          !: type of lateral boundary condition 
     41   REAL(wp), PUBLIC ::   e3zps_min       !: miminum thickness for partial steps (meters) 
     42   REAL(wp), PUBLIC ::   e3zps_rat       !: minimum thickness ration for partial steps 
    3843 
    3944   INTEGER, PUBLIC :: nn_interp 

  • NEMO/branches/2019/ENHANCE-02_ISF_domcfg/src/domain.F90

    r11201 r11308  
    105105         &             ln_cfmeta, ln_iscpl 
    106106      NAMELIST/namdom/ nn_bathy, cn_topo, cn_bath, cn_lon, cn_lat, nn_interp,                        & 
    107          &             rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, rn_isfhmin,           & 
    108          &             rn_atfp , rn_rdt   , ln_crs      , jphgr_msh ,                  & 
     107         &             rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin,                       & 
     108         &             rn_atfp , rn_rdt   , ln_crs      , jphgr_msh ,                                & 
    109109         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m,                         & 
    110110         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh,                  & 
    … …  
    204204         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin 
    205205         WRITE(numout,*) '      min number of ocean level (<0)       ' 
    206          WRITE(numout,*) '      treshold to open the isf cavity   rn_isfhmin   = ', rn_isfhmin, ' (m)' 
    207206         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)' 
    208207         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat 
    … …  
    415414         END DO 
    416415      END DO 
    417       CALL iom_rstput( 0, 0, inum, 'bathy_metry'   , z2d , ktype = jp_r4 ) 
     416      CALL iom_rstput( 0, 0, inum, 'bathy_metry_e3'   , z2d , ktype = jp_r4 ) 
     417      DO jj = 1,jpj 
     418         DO ji = 1,jpi 
     419            z2d (ji,jj) = SUM ( e3t_0(ji,jj, 1:mikt(ji,jj)-1 ) ) * ssmask(ji,jj)  
     420         END DO 
     421      END DO 
     422      CALL iom_rstput( 0, 0, inum, 'isf_draft_e3'   , z2d , ktype = jp_r4 ) 
     423      CALL iom_rstput( 0, 0, inum, 'isf_draft'   , risfdep , ktype = jp_r4 ) 
     424      CALL iom_rstput( 0, 0, inum, 'bathy_metry'   , bathy , ktype = jp_r4 ) 
     425      CALL iom_rstput( 0, 0, inum, 'hw',bathy-risfdep, ktype = jp_r4 ) 
     426      CALL iom_rstput( 0, 0, inum, 'mhw',mbkt*ssmask-mikt*ssmask, ktype = jp_i4 ) 
    418427      ! 
    419428      !                              !== closed sea ==! 

  • NEMO/branches/2019/ENHANCE-02_ISF_domcfg/src/dommsk.F90

    r11201 r11308  
    2424   !!---------------------------------------------------------------------- 
    2525   USE dom_oce        ! ocean space and time domain 
     26   USE domisf         ! domain: ice shelf 
    2627   USE domwri         ! domain: write the meshmask file 
    2728   USE bdy_oce        ! open boundary 
    … …  
    142143            END DO 
    143144         END DO 
    144       END DO 
    145       
     145      END DO     
     146  
     147      IF ( ln_isfsubgl ) CALL zgr_isf_subgl 
     148 
    146149!SF  add here lbc_lnk: bug not still understood : cause now domain configuration is read ! 
    147150!!gm I don't understand why...   

  • NEMO/branches/2019/ENHANCE-02_ISF_domcfg/src/domzgr.F90

    r11201 r11308  
    1717   !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn and Furner stretching function 
    1818   !!            3.4  ! 2012-12  (R. Bourdalle-Badie and G. Reffray)  modify C1D case   
    19    !!            3.6  ! 2014-11  (P. Mathiot and C. Harris) add ice shelf capabilitye 
     19   !!            3.6  ! 2014-11  (P. Mathiot and C. Harris) add ice shelf capability 
    2020   !!            3.?  ! 2015-11  (H. Liu) Modifications for Wetting/Drying 
    2121   !!---------------------------------------------------------------------- 
    … …  
    4343   USE lib_fortran 
    4444   USE dombat 
     45   USE domisf 
    4546 
    4647   IMPLICIT NONE 
    … …  
    4849 
    4950   PUBLIC   dom_zgr        ! called by dom_init.F90 
    50  
     51   ! 
    5152   !                              !!* Namelist namzgr_sco * 
    5253   LOGICAL  ::   ln_s_sh94         ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T) 
    … …  
    5758   REAL(wp) ::   rn_rmax           ! maximum cut-off r-value allowed (0<rn_rmax<1) 
    5859   REAL(wp) ::   rn_hc             ! Critical depth for transition from sigma to stretched coordinates 
    59    INTEGER , PUBLIC ::   ntopo           !: = 0/1 ,compute/read the bathymetry file 
    60    REAL(wp), PUBLIC ::   e3zps_min       !: miminum thickness for partial steps (meters) 
    61    REAL(wp), PUBLIC ::   e3zps_rat       !: minimum thickness ration for partial steps 
    62    INTEGER, PUBLIC  ::   nperio          !: type of lateral boundary condition 
    6360 
    6461   ! Song and Haidvogel 1994 stretching parameters 
    … …  
    148145      IF( ioptio /= 1 )   CALL ctl_stop( ' none or several vertical coordinate options used' ) 
    149146      ! 
     147      IF ( ln_isfcav ) CALL zgr_isf_nam 
    150148      ioptio = 0 
    151149      IF ( ln_zco .AND. ln_isfcav ) ioptio = ioptio + 1 
    … …  
    160158      IF( ln_zps      )   CALL zgr_zps          ! Partial step z-coordinate 
    161159      IF( ln_sco      )   CALL zgr_sco          ! s-coordinate or hybrid z-s coordinate 
     160                          CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isolated ocean points 
    162161      ! 
    163162      ! final adjustment of mbathy & check  
    164163      ! ----------------------------------- 
    165                           CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isolated ocean points 
    166164                          CALL zgr_bot_level    ! deepest ocean level for t-, u- and v-points 
    167165                          CALL zgr_top_level    ! shallowest ocean level for T-, U-, V- points 
    … …  
    313311      IF ( ln_isfcav .OR. ln_e3_dep ) THEN      ! e3. = dk[gdep]    
    314312         ! 
    315 !==>>>   need to be like this to compute the pressure gradient with ISF.  
    316 !        If not, level beneath the ISF are not aligned (sum(e3t) /= depth) 
    317 !        define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively 
    318 ! 
    319313         DO jk = 1, jpkm1 
    320314            e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk)  
    … …  
    578572               CALL iom_get  ( inum, jpdom_data, 'isf_draft', risfdep ) 
    579573               CALL iom_close( inum ) 
    580                WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp 
    581  
    582                ! set grounded point to 0  
    583                ! (a treshold could be set here if needed, or set it offline based on the grounded fraction) 
    584                WHERE ( bathy(:,:) <= risfdep(:,:) + rn_isfhmin ) 
    585                   misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 
    586                   mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp 
    587                END WHERE 
    588574            END IF 
    589575            !        
    … …  
    624610         ELSE                          ;   ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 )  ! from a depth 
    625611         ENDIF 
    626          zhmin = gdepw_1d(ik+1)                                                         ! minimum depth = ik+1 w-levels  
     612         zhmin = gdepw_1d(ik+1)                                                        ! minimum depth = ik+1 w-levels  
    627613         WHERE( bathy(:,:) <= 0._wp )   ;   bathy(:,:) = 0._wp                         ! min=0     over the lands 
    628614         ELSE WHERE                     ;   bathy(:,:) = MAX(  zhmin , bathy(:,:)  )   ! min=zhmin over the oceans 
    … …  
    931917      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used' 
    932918 
     919      ! compute position of the ice shelf grounding line 
     920      ! set bathy and isfdraft to 0 where grounded 
     921      IF ( ln_isfcav ) CALL zgr_isf_zspace 
     922 
    933923      ! bathymetry in level (from bathy_meter) 
    934924      ! =================== 
    … …  
    948938      END DO 
    949939 
     940      ! Check compatibility between bathy and iceshelf draft 
     941      ! insure at least 2 wet level on the vertical under an ice shelf 
     942      ! compute misfdep and adjust isf draft if needed 
     943      IF ( ln_isfcav ) CALL zgr_isf_kspace 
     944 
    950945      ! Scale factors and depth at T- and W-points 
    951946      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate 
    … …  
    956951      END DO 
    957952       
    958       ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf 
    959       IF ( ln_isfcav ) CALL zgr_isf 
    960  
    961953      ! Scale factors and depth at T- and W-points 
    962       IF ( .NOT. ln_isfcav ) THEN 
    963          DO jj = 1, jpj 
    964             DO ji = 1, jpi 
    965                ik = mbathy(ji,jj) 
    966                IF( ik > 0 ) THEN               ! ocean point only 
    967                   ! max ocean level case 
    968                   IF( ik == jpkm1 ) THEN 
    969                      zdepwp = bathy(ji,jj) 
    970                      ze3tp  = bathy(ji,jj) - gdepw_1d(ik) 
    971                      ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 
    972                      e3t_0(ji,jj,ik  ) = ze3tp 
    973                      e3t_0(ji,jj,ik+1) = ze3tp 
    974                      e3w_0(ji,jj,ik  ) = ze3wp 
    975                      e3w_0(ji,jj,ik+1) = ze3tp 
    976                      gdepw_0(ji,jj,ik+1) = zdepwp 
    977                      gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp 
    978                      gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 
    979                      ! 
    980                   ELSE                         ! standard case 
    981                      IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 
    982                      ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
    983                      ENDIF 
    984    !gm Bug?  check the gdepw_1d 
    985                      !       ... on ik 
    986                      gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   & 
    987                         &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
    988                         &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
    989                      e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
    990                         &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
    991                      e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
    992                         &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
    993                      !       ... on ik+1 
    994                      e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    995                      e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    996                      gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
    997                   ENDIF 
    998                ENDIF 
    999             END DO 
    1000          END DO 
    1001          ! 
    1002          it = 0 
    1003          DO jj = 1, jpj 
    1004             DO ji = 1, jpi 
    1005                ik = mbathy(ji,jj) 
    1006                IF( ik > 0 ) THEN               ! ocean point only 
    1007                   e3tp (ji,jj) = e3t_0(ji,jj,ik) 
    1008                   e3wp (ji,jj) = e3w_0(ji,jj,ik) 
    1009                   ! test 
    1010                   zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  ) 
    1011                   IF( zdiff <= 0._wp .AND. lwp ) THEN  
    1012                      it = it + 1 
    1013                      WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
    1014                      WRITE(numout,*) ' bathy = ', bathy(ji,jj) 
    1015                      WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
    1016                      WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  ) 
    1017                   ENDIF 
    1018                ENDIF 
    1019             END DO 
    1020          END DO 
    1021       END IF 
    1022       ! 
    1023       ! Scale factors and depth at U-, V-, UW and VW-points 
    1024       DO jk = 1, jpk                        ! initialisation to z-scale factors 
    1025          e3u_0 (:,:,jk) = e3t_1d(jk) 
    1026          e3v_0 (:,:,jk) = e3t_1d(jk) 
    1027          e3uw_0(:,:,jk) = e3w_1d(jk) 
    1028          e3vw_0(:,:,jk) = e3w_1d(jk) 
    1029       END DO 
    1030  
    1031       DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors 
    1032          DO jj = 1, jpjm1 
    1033             DO ji = 1, jpim1   ! vector opt. 
    1034                e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 
    1035                e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 
    1036                e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 
    1037                e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 
    1038             END DO 
    1039          END DO 
    1040       END DO 
    1041       IF ( ln_isfcav ) THEN 
    1042       ! (ISF) define e3uw (adapted for 2 cells in the water column) 
    1043          DO jj = 2, jpjm1  
    1044             DO ji = 2, jpim1   ! vector opt.  
    1045                ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) 
    1046                ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) 
    1047                IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) & 
    1048                                        &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) ) 
    1049                ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) 
    1050                ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) 
    1051                IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) & 
    1052                                        &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) ) 
    1053             END DO 
    1054          END DO 
    1055       END IF 
    1056  
    1057       CALL lbc_lnk('domzgr', e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk('domzgr', e3uw_0, 'U', 1._wp )   ! lateral boundary conditions 
    1058       CALL lbc_lnk('domzgr', e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp ) 
    1059       ! 
    1060  
    1061       DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
    1062          WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk) 
    1063          WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk) 
    1064          WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk) 
    1065          WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk) 
    1066       END DO 
    1067        
    1068       ! Scale factor at F-point 
    1069       DO jk = 1, jpk                        ! initialisation to z-scale factors 
    1070          e3f_0(:,:,jk) = e3t_1d(jk) 
    1071       END DO 
    1072       DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors 
    1073          DO jj = 1, jpjm1 
    1074             DO ji = 1, jpim1   ! vector opt. 
    1075                e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 
    1076             END DO 
    1077          END DO 
    1078       END DO 
    1079       CALL lbc_lnk('domzgr', e3f_0, 'F', 1._wp )       ! Lateral boundary conditions 
    1080       ! 
    1081       DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
    1082          WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk) 
    1083       END DO 
    1084 !!gm  bug ? :  must be a do loop with mj0,mj1 
    1085       !  
    1086       e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2 
    1087       e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)  
    1088       e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)  
    1089       e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)  
    1090       e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)  
    1091  
    1092       ! Control of the sign 
    1093       IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' ) 
    1094       IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' ) 
    1095       IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' ) 
    1096       IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' ) 
    1097       
    1098       ! Compute gde3w_0 (vertical sum of e3w) 
    1099       IF ( ln_isfcav ) THEN ! if cavity 
    1100          WHERE( misfdep == 0 )   misfdep = 1 
    1101          DO jj = 1,jpj 
    1102             DO ji = 1,jpi 
    1103                gde3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 
    1104                DO jk = 2, misfdep(ji,jj) 
    1105                   gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
    1106                END DO 
    1107                IF( misfdep(ji,jj) >= 2 )   gde3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 
    1108                DO jk = misfdep(ji,jj) + 1, jpk 
    1109                   gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
    1110                END DO 
    1111             END DO 
    1112          END DO 
    1113       ELSE ! no cavity 
    1114          gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 
    1115          DO jk = 2, jpk 
    1116             gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 
    1117          END DO 
    1118       END IF 
    1119       ! 
    1120       DEALLOCATE( zprt ) 
    1121       ! 
    1122    END SUBROUTINE zgr_zps 
    1123  
    1124  
    1125    SUBROUTINE zgr_isf 
    1126       !!---------------------------------------------------------------------- 
    1127       !!                    ***  ROUTINE zgr_isf  *** 
    1128       !!    
    1129       !! ** Purpose :   check the bathymetry in levels 
    1130       !!    
    1131       !! ** Method  :   THe water column have to contained at least 2 cells 
    1132       !!                Bathymetry and isfdraft are modified (dig/close) to respect 
    1133       !!                this criterion. 
    1134       !!    
    1135       !! ** Action  : - test compatibility between isfdraft and bathy  
    1136       !!              - bathy and isfdraft are modified 
    1137       !!---------------------------------------------------------------------- 
    1138       INTEGER  ::   ji, jj, jl, jk       ! dummy loop indices 
    1139       INTEGER  ::   ik, it               ! temporary integers 
    1140       INTEGER  ::   icompt, ibtest       ! (ISF) 
    1141       INTEGER  ::   ibtestim1, ibtestip1 ! (ISF) 
    1142       INTEGER  ::   ibtestjm1, ibtestjp1 ! (ISF) 
    1143       REAL(wp) ::   zdepth           ! Ajusted ocean depth to avoid too small e3t 
    1144       REAL(wp) ::   zmax             ! Maximum and minimum depth 
    1145       REAL(wp) ::   zbathydiff       ! isf temporary scalar 
    1146       REAL(wp) ::   zrisfdepdiff     ! isf temporary scalar 
    1147       REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
    1148       REAL(wp) ::   zdepwp           ! Ajusted ocean depth to avoid too small e3t 
    1149       REAL(wp) ::   zdiff            ! temporary scalar 
    1150       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 2D workspace (ISH) 
    1151       INTEGER , ALLOCATABLE, DIMENSION(:,:)   ::   zmbathy, zmisfdep         ! 2D workspace (ISH) 
    1152       !!--------------------------------------------------------------------- 
    1153       ! 
    1154       ALLOCATE( zbathy(jpi,jpj), zmask(jpi,jpj), zrisfdep(jpi,jpj) ) 
    1155       ALLOCATE( zmisfdep(jpi,jpj), zmbathy(jpi,jpj) ) 
    1156       ! 
    1157       ! (ISF) compute misfdep 
    1158       WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) /= 0 ) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1   
    1159       ELSEWHERE                      ;                         misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level  
    1160       END WHERE   
    1161  
    1162       ! Compute misfdep for ocean points (i.e. first wet level)  
    1163       ! find the first ocean level such that the first level thickness  
    1164       ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where  
    1165       ! e3t_0 is the reference level thickness  
    1166       DO jk = 2, jpkm1  
    1167          zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )  
    1168          WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1  
    1169       END DO  
    1170       WHERE ( 0._wp < risfdep(:,:) .AND. risfdep(:,:) <= e3t_1d(1) ) 
    1171          risfdep(:,:) = 0. ; misfdep(:,:) = 1 
    1172       END WHERE 
    1173  
    1174       ! remove very shallow ice shelf (less than ~ 10m if 75L) 
    1175       WHERE (risfdep(:,:) <= 10._wp .AND. misfdep(:,:) > 1) 
    1176          misfdep = 0; risfdep = 0.0_wp; 
    1177          mbathy  = 0; bathy   = 0.0_wp; 
    1178       END WHERE 
    1179       WHERE (bathy(:,:) <= 30.0_wp .AND. gphit < -60._wp) 
    1180          misfdep = 0; risfdep = 0.0_wp; 
    1181          mbathy  = 0; bathy   = 0.0_wp; 
    1182       END WHERE 
    1183   
    1184 ! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation 
    1185       icompt = 0  
    1186 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 
    1187       DO jl = 1, 10      
    1188          ! check at each iteration if isf is grounded or not (1cm treshold have to be update after first coupling experiments) 
    1189          WHERE (bathy(:,:) <= risfdep(:,:) + rn_isfhmin) 
    1190             misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 
    1191             mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp 
    1192          END WHERE 
    1193          WHERE (mbathy(:,:) <= 0)  
    1194             misfdep(:,:) = 0; risfdep(:,:) = 0._wp  
    1195             mbathy (:,:) = 0; bathy  (:,:) = 0._wp 
    1196          END WHERE 
    1197          IF( lk_mpp ) THEN 
    1198             zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
    1199             CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. ) 
    1200             misfdep(:,:) = INT( zbathy(:,:) ) 
    1201  
    1202             CALL lbc_lnk( 'domzgr',risfdep,'T', 1. ) 
    1203             CALL lbc_lnk( 'domzgr',bathy,  'T', 1. ) 
    1204  
    1205             zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
    1206             CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. ) 
    1207             mbathy(:,:)  = INT( zbathy(:,:) ) 
    1208          ENDIF 
    1209          IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN  
    1210             misfdep( 1 ,:) = misfdep(jpim1,:)            ! local domain is cyclic east-west  
    1211             misfdep(jpi,:) = misfdep(  2  ,:)  
    1212             mbathy( 1 ,:)  = mbathy(jpim1,:)             ! local domain is cyclic east-west 
    1213             mbathy(jpi,:)  = mbathy(  2  ,:) 
    1214          ENDIF 
    1215  
    1216          ! split last cell if possible (only where water column is 2 cell or less) 
    1217          ! if coupled to ice sheet, we do not modify the bathymetry (can be discuss). 
    1218          IF ( .NOT. ln_iscpl) THEN 
    1219             DO jk = jpkm1, 1, -1 
    1220                zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
    1221                WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 
    1222                   mbathy(:,:) = jk 
    1223                   bathy(:,:)  = zmax 
    1224                END WHERE 
    1225             END DO 
    1226          END IF 
    1227   
    1228          ! split top cell if possible (only where water column is 2 cell or less) 
    1229          DO jk = 2, jpkm1 
    1230             zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
    1231             WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy) 
    1232                misfdep(:,:) = jk 
    1233                risfdep(:,:) = zmax 
    1234             END WHERE 
    1235          END DO 
    1236  
    1237   
    1238  ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 
    1239          DO jj = 1, jpj 
    1240             DO ji = 1, jpi 
    1241                ! find the minimum change option: 
    1242                ! test bathy 
    1243                IF (risfdep(ji,jj) > 1) THEN 
    1244                   IF ( .NOT. ln_iscpl ) THEN 
    1245                      zbathydiff  =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) & 
    1246                          &            + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
    1247                      zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) & 
    1248                          &            - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
    1249                      IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) <  misfdep(ji,jj)) THEN 
    1250                         IF (zbathydiff <= zrisfdepdiff) THEN 
    1251                            bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 
    1252                            mbathy(ji,jj)= mbathy(ji,jj) + 1 
    1253                         ELSE 
    1254                            risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 
    1255                            misfdep(ji,jj) = misfdep(ji,jj) - 1 
    1256                         END IF 
    1257                      ENDIF 
    1258                   ELSE 
    1259                      IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) <  misfdep(ji,jj)) THEN 
    1260                         risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 
    1261                         misfdep(ji,jj) = misfdep(ji,jj) - 1 
    1262                      END IF 
    1263                   END IF 
    1264                END IF 
    1265             END DO 
    1266          END DO 
    1267   
    1268          ! At least 2 levels for water thickness at T, U, and V point. 
    1269          DO jj = 1, jpj 
    1270             DO ji = 1, jpi 
    1271                ! find the minimum change option: 
    1272                ! test bathy 
    1273                IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 
    1274                   IF ( .NOT. ln_iscpl ) THEN  
    1275                      zbathydiff  =ABS(bathy(ji,jj)   - ( gdepw_1d(mbathy (ji,jj)+1) & 
    1276                          &                             + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
    1277                      zrisfdepdiff=ABS(risfdep(ji,jj) - ( gdepw_1d(misfdep(ji,jj)  ) &  
    1278                          &                             - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
    1279                      IF (zbathydiff <= zrisfdepdiff) THEN 
    1280                         mbathy(ji,jj) = mbathy(ji,jj) + 1 
    1281                         bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
    1282                      ELSE 
    1283                         misfdep(ji,jj)= misfdep(ji,jj) - 1 
    1284                         risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 
    1285                      END IF 
    1286                   ELSE 
    1287                      misfdep(ji,jj)= misfdep(ji,jj) - 1 
    1288                      risfdep(ji,jj)= gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 
    1289                   END IF 
    1290                ENDIF 
    1291             END DO 
    1292          END DO 
    1293   
    1294  ! point V mbathy(ji,jj) == misfdep(ji,jj+1)  
    1295          DO jj = 1, jpjm1 
    1296             DO ji = 1, jpim1 
    1297                IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 
    1298                   IF ( .NOT. ln_iscpl ) THEN  
    1299                      zbathydiff  =ABS(bathy(ji,jj    ) - ( gdepw_1d(mbathy (ji,jj)+1) & 
    1300                           &                              + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat ))) 
    1301                      zrisfdepdiff=ABS(risfdep(ji,jj+1) - ( gdepw_1d(misfdep(ji,jj+1)) & 
    1302                           &                              - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 
    1303                      IF (zbathydiff <= zrisfdepdiff) THEN 
    1304                         mbathy(ji,jj) = mbathy(ji,jj) + 1 
    1305                         bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj   )+1)*e3zps_rat ) 
    1306                      ELSE 
    1307                         misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1 
    1308                         risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 
    1309                      END IF 
    1310                   ELSE 
    1311                      misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1 
    1312                      risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 
    1313                   END IF 
    1314                ENDIF 
    1315             END DO 
    1316          END DO 
    1317   
    1318          IF( lk_mpp ) THEN 
    1319             zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
    1320             CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. ) 
    1321             misfdep(:,:) = INT( zbathy(:,:) ) 
    1322  
    1323             CALL lbc_lnk( 'domzgr',risfdep,'T', 1. ) 
    1324             CALL lbc_lnk( 'domzgr',bathy,  'T', 1. ) 
    1325  
    1326             zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
    1327             CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. ) 
    1328             mbathy(:,:)  = INT( zbathy(:,:) ) 
    1329          ENDIF 
    1330  ! point V misdep(ji,jj) == mbathy(ji,jj+1)  
    1331          DO jj = 1, jpjm1 
    1332             DO ji = 1, jpim1 
    1333                IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) > 1) THEN 
    1334                   IF ( .NOT. ln_iscpl ) THEN  
    1335                      zbathydiff  =ABS(  bathy(ji,jj+1) - ( gdepw_1d(mbathy (ji,jj+1)+1) & 
    1336                            &                             + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 
    1337                      zrisfdepdiff=ABS(risfdep(ji,jj  ) - ( gdepw_1d(misfdep(ji,jj  )  ) & 
    1338                            &                             - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat ))) 
    1339                      IF (zbathydiff <= zrisfdepdiff) THEN 
    1340                         mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 
    1341                         bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 
    1342                      ELSE 
    1343                         misfdep(ji,jj)   = misfdep(ji,jj) - 1 
    1344                         risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat ) 
    1345                      END IF 
    1346                   ELSE 
    1347                      misfdep(ji,jj)   = misfdep(ji,jj) - 1 
    1348                      risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat ) 
    1349                   END IF 
    1350                ENDIF 
    1351             END DO 
    1352          END DO 
    1353   
    1354   
    1355          IF( lk_mpp ) THEN  
    1356             zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
    1357             CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. ) 
    1358             misfdep(:,:) = INT( zbathy(:,:) ) 
    1359  
    1360             CALL lbc_lnk( 'domzgr',risfdep,'T', 1. ) 
    1361             CALL lbc_lnk( 'domzgr',bathy,  'T', 1. ) 
    1362  
    1363             zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
    1364             CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. ) 
    1365             mbathy(:,:)  = INT( zbathy(:,:) ) 
    1366          ENDIF  
    1367   
    1368  ! point U mbathy(ji,jj) == misfdep(ji,jj+1)  
    1369          DO jj = 1, jpjm1 
    1370             DO ji = 1, jpim1 
    1371                IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 
    1372                   IF ( .NOT. ln_iscpl ) THEN  
    1373                   zbathydiff  =ABS(  bathy(ji  ,jj) - ( gdepw_1d(mbathy (ji,jj)+1) & 
    1374                        &                              + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat ))) 
    1375                   zrisfdepdiff=ABS(risfdep(ji+1,jj) - ( gdepw_1d(misfdep(ji+1,jj)) & 
    1376                        &                              - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 
    1377                   IF (zbathydiff <= zrisfdepdiff) THEN 
    1378                      mbathy(ji,jj) = mbathy(ji,jj) + 1 
    1379                      bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
    1380                   ELSE 
    1381                      misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 
    1382                      risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 
    1383                   END IF 
    1384                   ELSE 
    1385                      misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 
    1386                      risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 
    1387                   ENDIF 
    1388                ENDIF 
    1389             ENDDO 
    1390          ENDDO 
    1391   
    1392          IF( lk_mpp ) THEN  
    1393             zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
    1394             CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. ) 
    1395             misfdep(:,:) = INT( zbathy(:,:) ) 
    1396  
    1397             CALL lbc_lnk( 'domzgr',risfdep,'T', 1. ) 
    1398             CALL lbc_lnk( 'domzgr',bathy,  'T', 1. ) 
    1399  
    1400             zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
    1401             CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. ) 
    1402             mbathy(:,:)  = INT( zbathy(:,:) ) 
    1403          ENDIF  
    1404   
    1405  ! point U misfdep(ji,jj) == bathy(ji,jj+1)  
    1406          DO jj = 1, jpjm1 
    1407             DO ji = 1, jpim1 
    1408                IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) > 1) THEN 
    1409                   IF ( .NOT. ln_iscpl ) THEN  
    1410                      zbathydiff  =ABS(  bathy(ji+1,jj) - ( gdepw_1d(mbathy (ji+1,jj)+1) & 
    1411                           &                              + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 
    1412                      zrisfdepdiff=ABS(risfdep(ji  ,jj) - ( gdepw_1d(misfdep(ji  ,jj)  ) & 
    1413                           &                              - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat ))) 
    1414                      IF (zbathydiff <= zrisfdepdiff) THEN 
    1415                         mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1 
    1416                         bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 
    1417                      ELSE 
    1418                         misfdep(ji,jj)   = misfdep(ji  ,jj) - 1 
    1419                         risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat ) 
    1420                      END IF 
    1421                   ELSE 
    1422                      misfdep(ji,jj)   = misfdep(ji  ,jj) - 1 
    1423                      risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat ) 
    1424                   ENDIF 
    1425                ENDIF 
    1426             ENDDO 
    1427          ENDDO 
    1428   
    1429          IF( lk_mpp ) THEN 
    1430             zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
    1431             CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. ) 
    1432             misfdep(:,:) = INT( zbathy(:,:) ) 
    1433  
    1434             CALL lbc_lnk( 'domzgr',risfdep,'T', 1. ) 
    1435             CALL lbc_lnk( 'domzgr',bathy,  'T', 1. ) 
    1436  
    1437             zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
    1438             CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. ) 
    1439             mbathy(:,:)  = INT( zbathy(:,:) ) 
    1440          ENDIF 
    1441       END DO 
    1442       ! end dig bathy/ice shelf to be compatible 
    1443       ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness 
    1444       DO jl = 1,20 
    1445   
    1446  ! remove single point "bay" on isf coast line in the ice shelf draft' 
    1447          DO jk = 2, jpk 
    1448             WHERE (misfdep==0) misfdep=jpk 
    1449             zmask=0._wp 
    1450             WHERE (misfdep <= jk) zmask=1 
    1451             DO jj = 2, jpjm1 
    1452                DO ji = 2, jpim1 
    1453                   IF (misfdep(ji,jj) == jk) THEN 
    1454                      ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 
    1455                      IF (ibtest <= 1) THEN 
    1456                         risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1 
    1457                         IF (misfdep(ji,jj) > mbathy(ji,jj)) misfdep(ji,jj) = jpk 
    1458                      END IF 
    1459                   END IF 
    1460                END DO 
    1461             END DO 
    1462          END DO 
    1463          WHERE (misfdep==jpk) 
    1464              misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 
    1465          END WHERE 
    1466          IF( lk_mpp ) THEN 
    1467             zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
    1468             CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. ) 
    1469             misfdep(:,:) = INT( zbathy(:,:) ) 
    1470  
    1471             CALL lbc_lnk( 'domzgr',risfdep,'T', 1. ) 
    1472             CALL lbc_lnk('domzgr', bathy,  'T', 1. ) 
    1473  
    1474             zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
    1475             CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. ) 
    1476             mbathy(:,:)  = INT( zbathy(:,:) ) 
    1477          ENDIF 
    1478   
    1479  ! remove single point "bay" on bathy coast line beneath an ice shelf' 
    1480          DO jk = jpk,1,-1 
    1481             zmask=0._wp 
    1482             WHERE (mbathy >= jk ) zmask=1 
    1483             DO jj = 2, jpjm1 
    1484                DO ji = 2, jpim1 
    1485                   IF (mbathy(ji,jj) == jk .AND. misfdep(ji,jj) >= 2) THEN 
    1486                      ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 
    1487                      IF (ibtest <= 1) THEN 
    1488                         bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1 
    1489                         IF (misfdep(ji,jj) > mbathy(ji,jj)) mbathy(ji,jj) = 0 
    1490                      END IF 
    1491                   END IF 
    1492                END DO 
    1493             END DO 
    1494          END DO 
    1495          WHERE (mbathy==0) 
    1496              misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 
    1497          END WHERE 
    1498          IF( lk_mpp ) THEN 
    1499             zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
    1500             CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. ) 
    1501             misfdep(:,:) = INT( zbathy(:,:) ) 
    1502  
    1503             CALL lbc_lnk( 'domzgr',risfdep,'T', 1. ) 
    1504             CALL lbc_lnk( 'domzgr',bathy,  'T', 1. ) 
    1505  
    1506             zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
    1507             CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. ) 
    1508             mbathy(:,:)  = INT( zbathy(:,:) ) 
    1509          ENDIF 
    1510   
    1511  ! fill hole in ice shelf 
    1512          zmisfdep = misfdep 
    1513          zrisfdep = risfdep 
    1514          WHERE (zmisfdep <= 1._wp) zmisfdep=jpk 
    1515          DO jj = 2, jpjm1 
    1516             DO ji = 2, jpim1 
    1517                ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  ) 
    1518                ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1) 
    1519                IF( zmisfdep(ji,jj) >= mbathy(ji-1,jj  ) ) ibtestim1 = jpk 
    1520                IF( zmisfdep(ji,jj) >= mbathy(ji+1,jj  ) ) ibtestip1 = jpk 
    1521                IF( zmisfdep(ji,jj) >= mbathy(ji  ,jj-1) ) ibtestjm1 = jpk 
    1522                IF( zmisfdep(ji,jj) >= mbathy(ji  ,jj+1) ) ibtestjp1 = jpk 
    1523                ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
    1524                IF( ibtest == jpk .AND. misfdep(ji,jj) >= 2) THEN 
    1525                   mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp 
    1526                END IF 
    1527                IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) >= 2) THEN 
    1528                   misfdep(ji,jj) = ibtest 
    1529                   risfdep(ji,jj) = gdepw_1d(ibtest) 
    1530                ENDIF 
    1531             ENDDO 
    1532          ENDDO 
    1533   
    1534          IF( lk_mpp ) THEN  
    1535             zbathy(:,:)  = FLOAT( misfdep(:,:) )  
    1536             CALL lbc_lnk( 'domzgr',zbathy,  'T', 1. )  
    1537             misfdep(:,:) = INT( zbathy(:,:) )  
    1538  
    1539             CALL lbc_lnk( 'domzgr',risfdep, 'T', 1. )  
    1540             CALL lbc_lnk( 'domzgr',bathy,   'T', 1. ) 
    1541  
    1542             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1543             CALL lbc_lnk( 'domzgr',zbathy,  'T', 1. ) 
    1544             mbathy(:,:) = INT( zbathy(:,:) ) 
    1545          ENDIF  
    1546  ! 
    1547  !! fill hole in bathymetry 
    1548          zmbathy (:,:)=mbathy (:,:) 
    1549          DO jj = 2, jpjm1 
    1550             DO ji = 2, jpim1 
    1551                ibtestim1 = zmbathy(ji-1,jj  ) ; ibtestip1 = zmbathy(ji+1,jj  ) 
    1552                ibtestjm1 = zmbathy(ji  ,jj-1) ; ibtestjp1 = zmbathy(ji  ,jj+1) 
    1553                IF( zmbathy(ji,jj) <  misfdep(ji-1,jj  ) ) ibtestim1 = 0 
    1554                IF( zmbathy(ji,jj) <  misfdep(ji+1,jj  ) ) ibtestip1 = 0 
    1555                IF( zmbathy(ji,jj) <  misfdep(ji  ,jj-1) ) ibtestjm1 = 0 
    1556                IF( zmbathy(ji,jj) <  misfdep(ji  ,jj+1) ) ibtestjp1 = 0 
    1557                ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
    1558                IF( ibtest == 0 .AND. misfdep(ji,jj) >= 2) THEN 
    1559                   mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 
    1560                END IF 
    1561                IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) >= 2) THEN 
    1562                   mbathy(ji,jj) = ibtest 
    1563                   bathy(ji,jj)  = gdepw_1d(ibtest+1)  
    1564                ENDIF 
    1565             END DO 
    1566          END DO 
    1567          IF( lk_mpp ) THEN  
    1568             zbathy(:,:)  = FLOAT( misfdep(:,:) )  
    1569             CALL lbc_lnk( 'domzgr',zbathy,  'T', 1. )  
    1570             misfdep(:,:) = INT( zbathy(:,:) )  
    1571  
    1572             CALL lbc_lnk( 'domzgr',risfdep, 'T', 1. )  
    1573             CALL lbc_lnk( 'domzgr',bathy,   'T', 1. ) 
    1574  
    1575             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1576             CALL lbc_lnk( 'domzgr',zbathy,  'T', 1. ) 
    1577             mbathy(:,:) = INT( zbathy(:,:) ) 
    1578          ENDIF  
    1579  ! if not compatible after all check (ie U point water column less than 2 cells), mask U 
    1580          DO jj = 1, jpjm1 
    1581             DO ji = 1, jpim1 
    1582                IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 
    1583                   mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ; 
    1584                END IF 
    1585             END DO 
    1586          END DO 
    1587          IF( lk_mpp ) THEN  
    1588             zbathy(:,:)  = FLOAT( misfdep(:,:) )  
    1589             CALL lbc_lnk( 'domzgr',zbathy,  'T', 1. )  
    1590             misfdep(:,:) = INT( zbathy(:,:) )  
    1591  
    1592             CALL lbc_lnk('domzgr', risfdep, 'T', 1. )  
    1593             CALL lbc_lnk('domzgr', bathy,   'T', 1. ) 
    1594  
    1595             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1596             CALL lbc_lnk( 'domzgr',zbathy,  'T', 1. ) 
    1597             mbathy(:,:) = INT( zbathy(:,:) ) 
    1598          ENDIF  
    1599  ! if not compatible after all check (ie U point water column less than 2 cells), mask U 
    1600          DO jj = 1, jpjm1 
    1601             DO ji = 1, jpim1 
    1602                IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 
    1603                   mbathy(ji+1,jj)  = mbathy(ji+1,jj) - 1;   bathy(ji+1,jj)   = gdepw_1d(mbathy(ji+1,jj)+1) ; 
    1604                END IF 
    1605             END DO 
    1606          END DO 
    1607          IF( lk_mpp ) THEN  
    1608             zbathy(:,:)  = FLOAT( misfdep(:,:) )  
    1609             CALL lbc_lnk('domzgr', zbathy, 'T', 1. )  
    1610             misfdep(:,:) = INT( zbathy(:,:) )  
    1611  
    1612             CALL lbc_lnk('domzgr', risfdep,'T', 1. )  
    1613             CALL lbc_lnk( 'domzgr',bathy,  'T', 1. ) 
    1614  
    1615             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1616             CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. ) 
    1617             mbathy(:,:) = INT( zbathy(:,:) ) 
    1618          ENDIF  
    1619  ! if not compatible after all check (ie V point water column less than 2 cells), mask V 
    1620          DO jj = 1, jpjm1 
    1621             DO ji = 1, jpi 
    1622                IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 
    1623                   mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ; 
    1624                END IF 
    1625             END DO 
    1626          END DO 
    1627          IF( lk_mpp ) THEN  
    1628             zbathy(:,:)  = FLOAT( misfdep(:,:) )  
    1629             CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )  
    1630             misfdep(:,:) = INT( zbathy(:,:) )  
    1631  
    1632             CALL lbc_lnk( 'domzgr',risfdep,'T', 1. )  
    1633             CALL lbc_lnk('domzgr', bathy,  'T', 1. ) 
    1634  
    1635             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1636             CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. ) 
    1637             mbathy(:,:) = INT( zbathy(:,:) ) 
    1638          ENDIF  
    1639  ! if not compatible after all check (ie V point water column less than 2 cells), mask V 
    1640          DO jj = 1, jpjm1 
    1641             DO ji = 1, jpi 
    1642                IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 
    1643                   mbathy(ji,jj+1)  = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ; 
    1644                END IF 
    1645             END DO 
    1646          END DO 
    1647          IF( lk_mpp ) THEN  
    1648             zbathy(:,:)  = FLOAT( misfdep(:,:) )  
    1649             CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )  
    1650             misfdep(:,:) = INT( zbathy(:,:) )  
    1651  
    1652             CALL lbc_lnk( 'domzgr',risfdep,'T', 1. )  
    1653             CALL lbc_lnk( 'domzgr',bathy,  'T', 1. ) 
    1654  
    1655             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1656             CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. ) 
    1657             mbathy(:,:) = INT( zbathy(:,:) ) 
    1658          ENDIF  
    1659  ! if not compatible after all check, mask T 
    1660          DO jj = 1, jpj 
    1661             DO ji = 1, jpi 
    1662                IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN 
    1663                   misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj)  = 0 ; bathy(ji,jj)   = 0._wp ; 
    1664                END IF 
    1665             END DO 
    1666          END DO 
    1667   
    1668          WHERE (mbathy(:,:) == 1) 
    1669             mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp 
    1670          END WHERE 
    1671       END DO  
    1672 ! end check compatibility ice shelf/bathy 
    1673       ! remove very shallow ice shelf (less than ~ 10m if 75L) 
    1674       WHERE (risfdep(:,:) <= 10._wp) 
    1675          misfdep = 1; risfdep = 0.0_wp; 
    1676       END WHERE 
    1677  
    1678       IF( icompt == 0 ) THEN  
    1679          IF(lwp) WRITE(numout,*)'     no points with ice shelf too close to bathymetry'  
    1680       ELSE  
    1681          IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry'  
    1682       ENDIF  
    1683  
    1684       ! compute scale factor and depth at T- and W- points 
    1685954      DO jj = 1, jpj 
    1686955         DO ji = 1, jpi 
    … …  
    1704973                  ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
    1705974                  ENDIF 
    1706       !            gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
    1707975!gm Bug?  check the gdepw_1d 
    1708976                  !       ... on ik 
    … …  
    1710978                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
    1711979                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
    1712                   e3t_0  (ji,jj,ik  ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik  ) 
    1713                   e3w_0  (ji,jj,ik  ) = gdept_0(ji,jj,ik  ) - gdept_1d(ik-1) 
     980                  e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
     981                     &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
     982                  e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
     983                     &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
    1714984                  !       ... on ik+1 
    1715985                  e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    1716986                  e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     987                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
    1717988               ENDIF 
    1718989            ENDIF 
    … …  
    17401011      END DO 
    17411012      ! 
    1742       ! (ISF) Definition of e3t, u, v, w for ISF case 
    1743       DO jj = 1, jpj  
    1744          DO ji = 1, jpi  
    1745             ik = misfdep(ji,jj)  
    1746             IF( ik > 1 ) THEN               ! ice shelf point only  
    1747                IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik)  
    1748                gdepw_0(ji,jj,ik) = risfdep(ji,jj)  
    1749 !gm Bug?  check the gdepw_0  
    1750             !       ... on ik  
    1751                gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   &  
    1752                   &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   &  
    1753                   &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      )  
    1754                e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
    1755                e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
    1756  
    1757                IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)  
    1758                   e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
    1759                ENDIF  
    1760             !       ... on ik / ik-1  
    1761                e3w_0  (ji,jj,ik  ) = e3t_0  (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
    1762                gdept_0(ji,jj,ik-1) = gdept_0(ji,jj,ik) - e3w_0(ji,jj,ik) 
    1763                e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
    1764                e3w_0  (ji,jj,ik-1) = gdept_0(ji,jj,ik-1) - gdept_1d(ik-2) 
    1765                gdepw_0(ji,jj,ik-1) = gdepw_0(ji,jj,ik) - e3t_0(ji,jj,ik-1) 
    1766             ENDIF  
    1767          END DO  
    1768       END DO  
    1769     
    1770       it = 0  
    1771       DO jj = 1, jpj  
    1772          DO ji = 1, jpi  
    1773             ik = misfdep(ji,jj)  
    1774             IF( ik > 1 ) THEN               ! ice shelf point only  
    1775                e3tp (ji,jj) = e3t_0(ji,jj,ik  )  
    1776                e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )  
    1777             ! test  
    1778                zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  )  
    1779                IF( zdiff <= 0. .AND. lwp ) THEN   
    1780                   it = it + 1  
    1781                   WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj  
    1782                   WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)  
    1783                   WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff  
    1784                   WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj)  
    1785                ENDIF  
    1786             ENDIF  
    1787          END DO  
    1788       END DO  
    1789  
    1790       DEALLOCATE( zbathy, zmask, zrisfdep ) 
    1791       DEALLOCATE( zmisfdep, zmbathy ) 
    1792       ! 
    1793    END SUBROUTINE zgr_isf 
    1794  
     1013      ! compute top scale factor if ice shelf 
     1014      IF (ln_isfcav) CALL zps_isf 
     1015      ! 
     1016      ! Scale factors and depth at U-, V-, UW and VW-points 
     1017      DO jk = 1, jpk                        ! initialisation to z-scale factors 
     1018         e3u_0 (:,:,jk) = e3t_1d(jk) 
     1019         e3v_0 (:,:,jk) = e3t_1d(jk) 
     1020         e3uw_0(:,:,jk) = e3w_1d(jk) 
     1021         e3vw_0(:,:,jk) = e3w_1d(jk) 
     1022      END DO 
     1023 
     1024      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors 
     1025         DO jj = 1, jpjm1 
     1026            DO ji = 1, jpim1   ! vector opt. 
     1027               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 
     1028               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 
     1029               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 
     1030               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 
     1031            END DO 
     1032         END DO 
     1033      END DO 
     1034 
     1035      ! update e3uw in case only 2 cells in the water column 
     1036      IF ( ln_isfcav ) CALL zps_isf_e3uv_w 
     1037      ! 
     1038      CALL lbc_lnk('domzgr', e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk('domzgr', e3uw_0, 'U', 1._wp )   ! lateral boundary conditions 
     1039      CALL lbc_lnk('domzgr', e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp ) 
     1040      ! 
     1041      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
     1042         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk) 
     1043         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk) 
     1044         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk) 
     1045         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk) 
     1046      END DO 
     1047       
     1048      ! Scale factor at F-point 
     1049      DO jk = 1, jpk                        ! initialisation to z-scale factors 
     1050         e3f_0(:,:,jk) = e3t_1d(jk) 
     1051      END DO 
     1052      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors 
     1053         DO jj = 1, jpjm1 
     1054            DO ji = 1, jpim1   ! vector opt. 
     1055               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 
     1056            END DO 
     1057         END DO 
     1058      END DO 
     1059      CALL lbc_lnk('domzgr', e3f_0, 'F', 1._wp )       ! Lateral boundary conditions 
     1060      ! 
     1061      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
     1062         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk) 
     1063      END DO 
     1064!!gm  bug ? :  must be a do loop with mj0,mj1 
     1065      !  
     1066      e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2 
     1067      e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)  
     1068      e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)  
     1069      e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)  
     1070      e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)  
     1071 
     1072      ! Control of the sign 
     1073      IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' ) 
     1074      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' ) 
     1075      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' ) 
     1076      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' ) 
     1077      
     1078      ! Compute gde3w_0 (vertical sum of e3w) 
     1079      IF ( ln_isfcav ) THEN ! if cavity   ! to check if it can be simplify 
     1080         WHERE( misfdep == 0 )   misfdep = 1 
     1081         DO jj = 1,jpj 
     1082            DO ji = 1,jpi 
     1083               gde3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 
     1084               DO jk = 2, misfdep(ji,jj) 
     1085                  gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     1086               END DO 
     1087               IF( misfdep(ji,jj) >= 2 )   gde3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 
     1088               DO jk = misfdep(ji,jj) + 1, jpk 
     1089                  gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     1090               END DO 
     1091            END DO 
     1092         END DO 
     1093      ELSE ! no cavity 
     1094         gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 
     1095         DO jk = 2, jpk 
     1096            gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 
     1097         END DO 
     1098      END IF 
     1099      ! 
     1100      DEALLOCATE( zprt ) 
     1101      ! 
     1102   END SUBROUTINE zgr_zps 
    17951103 
    17961104   SUBROUTINE zgr_sco 
Note: See TracChangeset for help on using the changeset viewer.