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.
domzgr.F90 in branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 9746

Last change on this file since 9746 was 9746, checked in by mathiot, 4 years ago

isf geometry: add comments and minor cleaning

File size: 125.8 KB
Line 
1MODULE domzgr
2   !!==============================================================================
3   !!                       ***  MODULE domzgr   ***
4   !! Ocean initialization : domain initialization
5   !!==============================================================================
6   !! History :  OPA  ! 1995-12  (G. Madec)  Original code : s vertical coordinate
7   !!                 ! 1997-07  (G. Madec)  lbc_lnk call
8   !!                 ! 1997-04  (J.-O. Beismann)
9   !!            8.5  ! 2002-09  (A. Bozec, G. Madec)  F90: Free form and module
10   !!             -   ! 2002-09  (A. de Miranda)  rigid-lid + islands
11   !!  NEMO      1.0  ! 2003-08  (G. Madec)  F90: Free form and module
12   !!             -   ! 2005-10  (A. Beckmann)  modifications for hybrid s-ccordinates & new stretching function
13   !!            2.0  ! 2006-04  (R. Benshila, G. Madec)  add zgr_zco
14   !!            3.0  ! 2008-06  (G. Madec)  insertion of domzgr_zps.h90 & conding style
15   !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option
16   !!            3.3  ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level
17   !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn and Furner stretching function
18   !!            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 
20   !!----------------------------------------------------------------------
21
22   !!----------------------------------------------------------------------
23   !!   dom_zgr          : defined the ocean vertical coordinate system
24   !!       zgr_bat      : bathymetry fields (levels and meters)
25   !!       zgr_bat_zoom : modify the bathymetry field if zoom domain
26   !!       zgr_bat_ctl  : check the bathymetry files
27   !!       zgr_bot_level: deepest ocean level for t-, u, and v-points
28   !!       zgr_z        : reference z-coordinate
29   !!       zgr_zco      : z-coordinate
30   !!       zgr_zps      : z-coordinate with partial steps
31   !!       zgr_sco      : s-coordinate
32   !!       fssig        : tanh stretch function
33   !!       fssig1       : Song and Haidvogel 1994 stretch function
34   !!       fgamma       : Siddorn and Furner 2012 stretching function
35   !!---------------------------------------------------------------------
36   USE oce               ! ocean variables
37   USE dom_oce           ! ocean domain
38   USE closea            ! closed seas
39   USE c1d               ! 1D vertical configuration
40   USE in_out_manager    ! I/O manager
41   USE iom               ! I/O library
42   USE lbclnk            ! ocean lateral boundary conditions (or mpp link)
43   USE lib_mpp           ! distributed memory computing library
44   USE wrk_nemo          ! Memory allocation
45   USE timing            ! Timing
46
47   IMPLICIT NONE
48   PRIVATE
49
50   PUBLIC   dom_zgr        ! called by dom_init.F90
51
52   ! remove subglacial lake under the ice sheet.
53   LOGICAL , PUBLIC  :: ln_isfsubgl
54   REAL(wp), PUBLIC  :: rn_isfsubgllon, rn_isfsubgllat
55   !                              !!* Namelist namzgr_sco *
56   LOGICAL  ::   ln_s_sh94         ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T)
57   LOGICAL  ::   ln_s_sf12         ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T)
58   !
59   REAL(wp) ::   rn_sbot_min       ! minimum depth of s-bottom surface (>0) (m)
60   REAL(wp) ::   rn_sbot_max       ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
61   REAL(wp) ::   rn_rmax           ! maximum cut-off r-value allowed (0<rn_rmax<1)
62   REAL(wp) ::   rn_hc             ! Critical depth for transition from sigma to stretched coordinates
63   ! Song and Haidvogel 1994 stretching parameters
64   REAL(wp) ::   rn_theta          ! surface control parameter (0<=rn_theta<=20)
65   REAL(wp) ::   rn_thetb          ! bottom control parameter  (0<=rn_thetb<= 1)
66   REAL(wp) ::   rn_bb             ! stretching parameter
67   !                                        ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
68   ! Siddorn and Furner stretching parameters
69   LOGICAL  ::   ln_sigcrit        ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
70   REAL(wp) ::   rn_alpha          ! control parameter ( > 1 stretch towards surface, < 1 towards seabed)
71   REAL(wp) ::   rn_efold          !  efold length scale for transition to stretched coord
72   REAL(wp) ::   rn_zs             !  depth of surface grid box
73                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
74   REAL(wp) ::   rn_zb_a           !  bathymetry scaling factor for calculating Zb
75   REAL(wp) ::   rn_zb_b           !  offset for calculating Zb
76
77  !! * Substitutions
78#  include "domzgr_substitute.h90"
79#  include "vectopt_loop_substitute.h90"
80   !!----------------------------------------------------------------------
81   !! NEMO/OPA 3.3.1 , NEMO Consortium (2011)
82   !! $Id$
83   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
84   !!----------------------------------------------------------------------
85CONTAINS       
86
87   SUBROUTINE dom_zgr
88      !!----------------------------------------------------------------------
89      !!                ***  ROUTINE dom_zgr  ***
90      !!                   
91      !! ** Purpose :   set the depth of model levels and the resulting
92      !!              vertical scale factors.
93      !!
94      !! ** Method  : - reference 1D vertical coordinate (gdep._1d, e3._1d)
95      !!              - read/set ocean depth and ocean levels (bathy, mbathy)
96      !!              - vertical coordinate (gdep., e3.) depending on the
97      !!                coordinate chosen :
98      !!                   ln_zco=T   z-coordinate   
99      !!                   ln_zps=T   z-coordinate with partial steps
100      !!                   ln_zco=T   s-coordinate
101      !!
102      !! ** Action  :   define gdep., e3., mbathy and bathy
103      !!----------------------------------------------------------------------
104      INTEGER ::   ioptio, ibat   ! local integer
105      INTEGER ::   ios
106      !
107      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav
108      !!----------------------------------------------------------------------
109      !
110      IF( nn_timing == 1 )   CALL timing_start('dom_zgr')
111      !
112      REWIND( numnam_ref )              ! Namelist namzgr in reference namelist : Vertical coordinate
113      READ  ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901 )
114901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist', lwp )
115
116      REWIND( numnam_cfg )              ! Namelist namzgr in configuration namelist : Vertical coordinate
117      READ  ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 )
118902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp )
119      IF(lwm) WRITE ( numond, namzgr )
120
121      IF(lwp) THEN                     ! Control print
122         WRITE(numout,*)
123         WRITE(numout,*) 'dom_zgr : vertical coordinate'
124         WRITE(numout,*) '~~~~~~~'
125         WRITE(numout,*) '          Namelist namzgr : set vertical coordinate'
126         WRITE(numout,*) '             z-coordinate - full steps      ln_zco    = ', ln_zco
127         WRITE(numout,*) '             z-coordinate - partial steps   ln_zps    = ', ln_zps
128         WRITE(numout,*) '             s- or hybrid z-s-coordinate    ln_sco    = ', ln_sco
129         WRITE(numout,*) '             ice shelf cavities             ln_isfcav = ', ln_isfcav
130      ENDIF
131
132      ioptio = 0                       ! Check Vertical coordinate options
133      IF( ln_zco      )   ioptio = ioptio + 1
134      IF( ln_zps      )   ioptio = ioptio + 1
135      IF( ln_sco      )   ioptio = ioptio + 1
136      IF( ioptio /= 1 )   CALL ctl_stop( ' none or several vertical coordinate options used' )
137      !
138      ! Build the vertical coordinate system
139      ! ------------------------------------
140                          CALL zgr_z            ! Reference z-coordinate system (always called)
141                          CALL zgr_bat          ! Bathymetry fields (levels and meters)
142      IF( lk_c1d      )   CALL lbc_lnk( bathy , 'T', 1._wp )   ! 1D config.: same bathy value over the 3x3 domain
143      IF( ln_zco      )   CALL zgr_zco          ! z-coordinate
144      IF( ln_zps      )   CALL zgr_zps          ! Partial step z-coordinate
145      IF( ln_sco      )   CALL zgr_sco          ! s-coordinate or hybrid z-s coordinate
146      !
147      ! final adjustment of mbathy & check
148      ! -----------------------------------
149      IF( lzoom       )   CALL zgr_bat_zoom     ! correct mbathy in case of zoom subdomain
150      IF( .NOT.lk_c1d )   CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isolated ocean points
151                          CALL zgr_bot_level    ! deepest ocean level for t-, u- and v-points
152                          CALL zgr_top_level    ! shallowest ocean level for T-, U-, V- points
153      !
154      IF( lk_c1d ) THEN                         ! 1D config.: same mbathy value over the 3x3 domain
155         ibat = mbathy(2,2)
156         mbathy(:,:) = ibat
157      END IF
158      !
159      IF( nprint == 1 .AND. lwp )   THEN
160         WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
161         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
162            &                   ' w ',   MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gdep3w_0(:,:,:) )
163         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL( e3f_0(:,:,:) ),  &
164            &                   ' u ',   MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL( e3v_0(:,:,:) ),  &
165            &                   ' uw',   MINVAL( e3uw_0(:,:,:)), ' vw', MINVAL( e3vw_0(:,:,:)),   &
166            &                   ' w ',   MINVAL( e3w_0(:,:,:) )
167
168         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
169            &                   ' w ',   MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gdep3w_0(:,:,:) )
170         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL( e3f_0(:,:,:) ),  &
171            &                   ' u ',   MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL( e3v_0(:,:,:) ),  &
172            &                   ' uw',   MAXVAL( e3uw_0(:,:,:)), ' vw', MAXVAL( e3vw_0(:,:,:)),   &
173            &                   ' w ',   MAXVAL( e3w_0(:,:,:) )
174      ENDIF
175      !
176      IF( nn_timing == 1 )  CALL timing_stop('dom_zgr')
177      !
178   END SUBROUTINE dom_zgr
179
180
181   SUBROUTINE zgr_z
182      !!----------------------------------------------------------------------
183      !!                   ***  ROUTINE zgr_z  ***
184      !!                   
185      !! ** Purpose :   set the depth of model levels and the resulting
186      !!      vertical scale factors.
187      !!
188      !! ** Method  :   z-coordinate system (use in all type of coordinate)
189      !!        The depth of model levels is defined from an analytical
190      !!      function the derivative of which gives the scale factors.
191      !!        both depth and scale factors only depend on k (1d arrays).
192      !!              w-level: gdepw_1d  = gdep(k)
193      !!                       e3w_1d(k) = dk(gdep)(k)     = e3(k)
194      !!              t-level: gdept_1d  = gdep(k+0.5)
195      !!                       e3t_1d(k) = dk(gdep)(k+0.5) = e3(k+0.5)
196      !!
197      !! ** Action  : - gdept_1d, gdepw_1d : depth of T- and W-point (m)
198      !!              - e3t_1d  , e3w_1d   : scale factors at T- and W-levels (m)
199      !!
200      !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
201      !!----------------------------------------------------------------------
202      INTEGER  ::   jk                     ! dummy loop indices
203      REAL(wp) ::   zt, zw                 ! temporary scalars
204      REAL(wp) ::   zsur, za0, za1, zkth   ! Values set from parameters in
205      REAL(wp) ::   zacr, zdzmin, zhmax    ! par_CONFIG_Rxx.h90
206      REAL(wp) ::   zrefdep                ! depth of the reference level (~10m)
207      REAL(wp) ::   za2, zkth2, zacr2      ! Values for optional double tanh function set from parameters
208      !!----------------------------------------------------------------------
209      !
210      IF( nn_timing == 1 )  CALL timing_start('zgr_z')
211      !
212      ! Set variables from parameters
213      ! ------------------------------
214       zkth = ppkth       ;   zacr = ppacr
215       zdzmin = ppdzmin   ;   zhmax = pphmax
216       zkth2 = ppkth2     ;   zacr2 = ppacr2   ! optional (ldbletanh=T) double tanh parameters
217
218      ! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
219      !  za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
220      IF(   ppa1  == pp_to_be_computed  .AND.  &
221         &  ppa0  == pp_to_be_computed  .AND.  &
222         &  ppsur == pp_to_be_computed           ) THEN
223         !
224#if defined key_agrif
225         za1  = (  ppdzmin - pphmax / FLOAT(jpkdta-1)  )                                                   &
226            & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpkdta-1) * (  LOG( COSH( (jpkdta - ppkth) / ppacr) )&
227            &                                                      - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
228#else
229         za1  = (  ppdzmin - pphmax / FLOAT(jpkm1)  )                                                      &
230            & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * (  LOG( COSH( (jpk - ppkth) / ppacr) )      &
231            &                                                   - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
232#endif
233         za0  = ppdzmin - za1 *              TANH( (1-ppkth) / ppacr )
234         zsur =   - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
235      ELSE
236         za1 = ppa1 ;       za0 = ppa0 ;          zsur = ppsur
237         za2 = ppa2                            ! optional (ldbletanh=T) double tanh parameter
238      ENDIF
239
240      IF(lwp) THEN                         ! Parameter print
241         WRITE(numout,*)
242         WRITE(numout,*) '    zgr_z   : Reference vertical z-coordinates'
243         WRITE(numout,*) '    ~~~~~~~'
244         IF(  ppkth == 0._wp ) THEN             
245              WRITE(numout,*) '            Uniform grid with ',jpk-1,' layers'
246              WRITE(numout,*) '            Total depth    :', zhmax
247#if defined key_agrif
248              WRITE(numout,*) '            Layer thickness:', zhmax/(jpkdta-1)
249#else
250              WRITE(numout,*) '            Layer thickness:', zhmax/(jpk-1)
251#endif
252         ELSE
253            IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN
254               WRITE(numout,*) '         zsur, za0, za1 computed from '
255               WRITE(numout,*) '                 zdzmin = ', zdzmin
256               WRITE(numout,*) '                 zhmax  = ', zhmax
257            ENDIF
258            WRITE(numout,*) '           Value of coefficients for vertical mesh:'
259            WRITE(numout,*) '                 zsur = ', zsur
260            WRITE(numout,*) '                 za0  = ', za0
261            WRITE(numout,*) '                 za1  = ', za1
262            WRITE(numout,*) '                 zkth = ', zkth
263            WRITE(numout,*) '                 zacr = ', zacr
264            IF( ldbletanh ) THEN
265               WRITE(numout,*) ' (Double tanh    za2  = ', za2
266               WRITE(numout,*) '  parameters)    zkth2= ', zkth2
267               WRITE(numout,*) '                 zacr2= ', zacr2
268            ENDIF
269         ENDIF
270      ENDIF
271
272
273      ! Reference z-coordinate (depth - scale factor at T- and W-points)
274      ! ======================
275      IF( ppkth == 0._wp ) THEN            !  uniform vertical grid
276#if defined key_agrif
277         za1 = zhmax / FLOAT(jpkdta-1) 
278#else
279         za1 = zhmax / FLOAT(jpk-1) 
280#endif
281         DO jk = 1, jpk
282            zw = FLOAT( jk )
283            zt = FLOAT( jk ) + 0.5_wp
284            gdepw_1d(jk) = ( zw - 1 ) * za1
285            gdept_1d(jk) = ( zt - 1 ) * za1
286            e3w_1d  (jk) =  za1
287            e3t_1d  (jk) =  za1
288         END DO
289      ELSE                                ! Madec & Imbard 1996 function
290         IF( .NOT. ldbletanh ) THEN
291            DO jk = 1, jpk
292               zw = REAL( jk , wp )
293               zt = REAL( jk , wp ) + 0.5_wp
294               gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) )  )
295               gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) )  )
296               e3w_1d  (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   )
297               e3t_1d  (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   )
298            END DO
299         ELSE
300            DO jk = 1, jpk
301               zw = FLOAT( jk )
302               zt = FLOAT( jk ) + 0.5_wp
303               ! Double tanh function
304               gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr  ) )    &
305                  &                             + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) )  )
306               gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr  ) )    &
307                  &                             + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) )  )
308               e3w_1d  (jk) =          za0      + za1        * TANH(       (zw-zkth ) / zacr  )      &
309                  &                             + za2        * TANH(       (zw-zkth2) / zacr2 )
310               e3t_1d  (jk) =          za0      + za1        * TANH(       (zt-zkth ) / zacr  )      &
311                  &                             + za2        * TANH(       (zt-zkth2) / zacr2 )
312            END DO
313         ENDIF
314         gdepw_1d(1) = 0._wp                    ! force first w-level to be exactly at zero
315      ENDIF
316
317      IF ( ln_isfcav ) THEN
318! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth)
319! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively
320         DO jk = 1, jpkm1
321            e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) 
322         END DO
323         e3t_1d(jpk) = e3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO
324
325         DO jk = 2, jpk
326            e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1) 
327         END DO
328         e3w_1d(1  ) = 2._wp * (gdept_1d(1) - gdepw_1d(1)) 
329      END IF
330
331!!gm BUG in s-coordinate this does not work!
332      ! deepest/shallowest W level Above/Below ~10m
333      zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d )                   ! ref. depth with tolerance (10% of minimum layer thickness)
334      nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
335      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m
336!!gm end bug
337
338      IF(lwp) THEN                        ! control print
339         WRITE(numout,*)
340         WRITE(numout,*) '              Reference z-coordinate depth and scale factors:'
341         WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
342         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk )
343      ENDIF
344      DO jk = 1, jpk                      ! control positivity
345         IF( e3w_1d  (jk) <= 0._wp .OR. e3t_1d  (jk) <= 0._wp )   CALL ctl_stop( 'dom:zgr_z: e3w_1d or e3t_1d =< 0 '    )
346         IF( gdepw_1d(jk) <  0._wp .OR. gdept_1d(jk) <  0._wp )   CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' )
347      END DO
348      !
349      IF( nn_timing == 1 )  CALL timing_stop('zgr_z')
350      !
351   END SUBROUTINE zgr_z
352
353
354   SUBROUTINE zgr_bat
355      !!----------------------------------------------------------------------
356      !!                    ***  ROUTINE zgr_bat  ***
357      !!
358      !! ** Purpose :   set bathymetry both in levels and meters
359      !!
360      !! ** Method  :   read or define mbathy and bathy arrays
361      !!       * level bathymetry:
362      !!      The ocean basin geometry is given by a two-dimensional array,
363      !!      mbathy, which is defined as follow :
364      !!            mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
365      !!                              at t-point (ji,jj).
366      !!                            = 0  over the continental t-point.
367      !!      The array mbathy is checked to verified its consistency with
368      !!      model option. in particular:
369      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
370      !!                  along closed boundary.
371      !!            mbathy must be cyclic IF jperio=1.
372      !!            mbathy must be lower or equal to jpk-1.
373      !!            isolated ocean grid points are suppressed from mbathy
374      !!                  since they are only connected to remaining
375      !!                  ocean through vertical diffusion.
376      !!      ntopo=-1 :   rectangular channel or bassin with a bump
377      !!      ntopo= 0 :   flat rectangular channel or basin
378      !!      ntopo= 1 :   mbathy is read in 'bathy_level.nc' NetCDF file
379      !!                   bathy  is read in 'bathy_meter.nc' NetCDF file
380      !!
381      !! ** Action  : - mbathy: level bathymetry (in level index)
382      !!              - bathy : meter bathymetry (in meters)
383      !!----------------------------------------------------------------------
384      INTEGER  ::   ji, jj, jl, jk            ! dummy loop indices
385      INTEGER  ::   inum                      ! temporary logical unit
386      INTEGER  ::   ierror                    ! error flag
387      INTEGER  ::   ii_bump, ij_bump, ih      ! bump center position
388      INTEGER  ::   ii0, ii1, ij0, ij1, ik    ! local indices
389      REAL(wp) ::   r_bump , h_bump , h_oce   ! bump characteristics
390      REAL(wp) ::   zi, zj, zh, zhmin         ! local scalars
391      INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   idta   ! global domain integer data
392      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdta   ! global domain scalar data
393      !!----------------------------------------------------------------------
394      !
395      IF( nn_timing == 1 )  CALL timing_start('zgr_bat')
396      !
397      IF(lwp) WRITE(numout,*)
398      IF(lwp) WRITE(numout,*) '    zgr_bat : defines level and meter bathymetry'
399      IF(lwp) WRITE(numout,*) '    ~~~~~~~'
400      !
401      ! (ISF) initialisation ice shelf draft and top level
402      risfdep(:,:)=0._wp
403      misfdep(:,:)=1
404      !                                               ! ================== !
405      IF( ntopo == 0 .OR. ntopo == -1 ) THEN          !   defined by hand  !
406         !                                            ! ================== !
407         !                                            ! global domain level and meter bathymetry (idta,zdta)
408         !
409         ALLOCATE( idta(jpidta,jpjdta), STAT=ierror )
410         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate idta array' )
411         ALLOCATE( zdta(jpidta,jpjdta), STAT=ierror )
412         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zdta array' )
413         !
414         IF( ntopo == 0 ) THEN                        ! flat basin
415            IF(lwp) WRITE(numout,*)
416            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin'
417            IF( rn_bathy > 0.01 ) THEN
418               IF(lwp) WRITE(numout,*) '         Depth = rn_bathy read in namelist'
419               zdta(:,:) = rn_bathy
420               IF( ln_sco ) THEN                                   ! s-coordinate (zsc       ): idta()=jpk
421                  idta(:,:) = jpkm1
422               ELSE                                                ! z-coordinate (zco or zps): step-like topography
423                  idta(:,:) = jpkm1
424                  DO jk = 1, jpkm1
425                     WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) )   idta(:,:) = jk
426                  END DO
427               ENDIF
428            ELSE
429               IF(lwp) WRITE(numout,*) '         Depth = depthw(jpkm1)'
430               idta(:,:) = jpkm1                            ! before last level
431               zdta(:,:) = gdepw_1d(jpk)                     ! last w-point depth
432               h_oce     = gdepw_1d(jpk)
433            ENDIF
434         ELSE                                         ! bump centered in the basin
435            IF(lwp) WRITE(numout,*)
436            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin with a bump'
437            ii_bump = jpidta / 2                           ! i-index of the bump center
438            ij_bump = jpjdta / 2                           ! j-index of the bump center
439            r_bump  = 50000._wp                            ! bump radius (meters)       
440            h_bump  =  2700._wp                            ! bump height (meters)
441            h_oce   = gdepw_1d(jpk)                        ! background ocean depth (meters)
442            IF(lwp) WRITE(numout,*) '            bump characteristics: '
443            IF(lwp) WRITE(numout,*) '               bump center (i,j)   = ', ii_bump, ii_bump
444            IF(lwp) WRITE(numout,*) '               bump height         = ', h_bump , ' meters'
445            IF(lwp) WRITE(numout,*) '               bump radius         = ', r_bump , ' index'
446            IF(lwp) WRITE(numout,*) '            background ocean depth = ', h_oce  , ' meters'
447            !                                       
448            DO jj = 1, jpjdta                              ! zdta :
449               DO ji = 1, jpidta
450                  zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
451                  zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
452                  zdta(ji,jj) = h_oce - h_bump * EXP( -( zi*zi + zj*zj ) )
453               END DO
454            END DO
455            !                                              ! idta :
456            IF( ln_sco ) THEN                                   ! s-coordinate (zsc       ): idta()=jpk
457               idta(:,:) = jpkm1
458            ELSE                                                ! z-coordinate (zco or zps): step-like topography
459               idta(:,:) = jpkm1
460               DO jk = 1, jpkm1
461                  WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) )   idta(:,:) = jk
462               END DO
463            ENDIF
464         ENDIF
465         !                                            ! set GLOBAL boundary conditions
466         !                                            ! Caution : idta on the global domain: use of jperio, not nperio
467         IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
468            idta( :    , 1    ) = -1                ;      zdta( :    , 1    ) = -1._wp
469            idta( :    ,jpjdta) =  0                ;      zdta( :    ,jpjdta) =  0._wp
470         ELSEIF( jperio == 2 ) THEN
471            idta( :    , 1    ) = idta( : ,  3  )   ;      zdta( :    , 1    ) = zdta( : ,  3  )
472            idta( :    ,jpjdta) = 0                 ;      zdta( :    ,jpjdta) =  0._wp
473            idta( 1    , :    ) = 0                 ;      zdta( 1    , :    ) =  0._wp
474            idta(jpidta, :    ) = 0                 ;      zdta(jpidta, :    ) =  0._wp
475         ELSE
476            ih = 0                                  ;      zh = 0._wp
477            IF( ln_sco )   ih = jpkm1               ;      IF( ln_sco )   zh = h_oce
478            idta( :    , 1    ) = ih                ;      zdta( :    , 1    ) =  zh
479            idta( :    ,jpjdta) = ih                ;      zdta( :    ,jpjdta) =  zh
480            idta( 1    , :    ) = ih                ;      zdta( 1    , :    ) =  zh
481            idta(jpidta, :    ) = ih                ;      zdta(jpidta, :    ) =  zh
482         ENDIF
483
484         !                                            ! local domain level and meter bathymetries (mbathy,bathy)
485         mbathy(:,:) = 0                                   ! set to zero extra halo points
486         bathy (:,:) = 0._wp                               ! (require for mpp case)
487         DO jj = 1, nlcj                                   ! interior values
488            DO ji = 1, nlci
489               mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
490               bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
491            END DO
492         END DO
493         !
494         DEALLOCATE( idta, zdta )
495         !
496         !                                            ! ================ !
497      ELSEIF( ntopo == 1 ) THEN                       !   read in file   ! (over the local domain)
498         !                                            ! ================ !
499         !
500         IF( ln_zco )   THEN                          ! zco : read level bathymetry
501            CALL iom_open ( 'bathy_level.nc', inum ) 
502            CALL iom_get  ( inum, jpdom_data, 'Bathy_level', bathy )
503            CALL iom_close( inum )
504            mbathy(:,:) = INT( bathy(:,:) )
505            !                                                ! =====================
506            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
507               !                                             ! =====================
508               IF( nn_cla == 0 ) THEN
509                  ii0 = 140   ;   ii1 = 140                  ! Gibraltar Strait open
510                  ij0 = 102   ;   ij1 = 102                  ! (Thomson, Ocean Modelling, 1995)
511                  DO ji = mi0(ii0), mi1(ii1)
512                     DO jj = mj0(ij0), mj1(ij1)
513                        mbathy(ji,jj) = 15
514                     END DO
515                  END DO
516                  IF(lwp) WRITE(numout,*)
517                  IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
518                  !
519                  ii0 = 160   ;   ii1 = 160                  ! Bab el mandeb Strait open
520                  ij0 = 88    ;   ij1 = 88                   ! (Thomson, Ocean Modelling, 1995)
521                  DO ji = mi0(ii0), mi1(ii1)
522                     DO jj = mj0(ij0), mj1(ij1)
523                        mbathy(ji,jj) = 12
524                     END DO
525                  END DO
526                  IF(lwp) WRITE(numout,*)
527                  IF(lwp) WRITE(numout,*) '      orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
528               ENDIF
529               !
530            ENDIF
531            !
532         ENDIF
533         IF( ln_zps .OR. ln_sco )   THEN              ! zps or sco : read meter bathymetry
534            CALL iom_open ( 'bathy_meter.nc', inum ) 
535            IF ( ln_isfcav ) THEN
536               CALL iom_get  ( inum, jpdom_data, 'Bathymetry_isf', bathy, lrowattr=.false. )
537            ELSE
538               CALL iom_get  ( inum, jpdom_data, 'Bathymetry'    , bathy, lrowattr=ln_use_jattr  )
539            END IF
540            CALL iom_close( inum )
541            !                                               
542            IF ( ln_isfcav ) THEN
543               CALL iom_open ( 'isf_draft_meter.nc', inum ) 
544               CALL iom_get  ( inum, jpdom_data, 'isf_draft', risfdep )
545               CALL iom_close( inum )
546               WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp
547            END IF
548            !       
549            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
550               !
551              IF( nn_cla == 0 ) THEN
552                 ii0 = 140   ;   ii1 = 140                   ! Gibraltar Strait open
553                 ij0 = 102   ;   ij1 = 102                   ! (Thomson, Ocean Modelling, 1995)
554                 DO ji = mi0(ii0), mi1(ii1)
555                    DO jj = mj0(ij0), mj1(ij1)
556                       bathy(ji,jj) = 284._wp
557                    END DO
558                 END DO
559                 IF(lwp) WRITE(numout,*)     
560                 IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
561                 !
562                 ii0 = 160   ;   ii1 = 160                   ! Bab el mandeb Strait open
563                 ij0 = 88    ;   ij1 = 88                    ! (Thomson, Ocean Modelling, 1995)
564                 DO ji = mi0(ii0), mi1(ii1)
565                    DO jj = mj0(ij0), mj1(ij1)
566                       bathy(ji,jj) = 137._wp
567                    END DO
568                 END DO
569                 IF(lwp) WRITE(numout,*)
570                 IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
571              ENDIF
572              !
573           ENDIF
574            !
575        ENDIF
576         !                                            ! =============== !
577      ELSE                                            !      error      !
578         !                                            ! =============== !
579         WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
580         CALL ctl_stop( '    zgr_bat : '//trim(ctmp1) )
581      ENDIF
582      !
583      IF( nn_closea == 0 )   CALL clo_bat( bathy, mbathy )    !==  NO closed seas or lakes  ==!
584      !                       
585      IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==!
586         ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin
587         IF ( ln_isfcav ) THEN
588            WHERE (bathy == risfdep)
589               bathy   = 0.0_wp ; risfdep = 0.0_wp
590            END WHERE
591         END IF
592         ! end patch
593         IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level
594         ELSE                          ;   ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 )  ! from a depth
595         ENDIF
596         zhmin = gdepw_1d(ik+1)                                                         ! minimum depth = ik+1 w-levels
597         WHERE( bathy(:,:) <= 0._wp )   ;   bathy(:,:) = 0._wp                         ! min=0     over the lands
598         ELSE WHERE                     ;   bathy(:,:) = MAX(  zhmin , bathy(:,:)  )   ! min=zhmin over the oceans
599         END WHERE
600         IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
601      ENDIF
602      !
603      IF( nn_timing == 1 )  CALL timing_stop('zgr_bat')
604      !
605   END SUBROUTINE zgr_bat
606
607
608   SUBROUTINE zgr_bat_zoom
609      !!----------------------------------------------------------------------
610      !!                    ***  ROUTINE zgr_bat_zoom  ***
611      !!
612      !! ** Purpose : - Close zoom domain boundary if necessary
613      !!              - Suppress Med Sea from ORCA R2 and R05 arctic zoom
614      !!
615      !! ** Method  :
616      !!
617      !! ** Action  : - update mbathy: level bathymetry (in level index)
618      !!----------------------------------------------------------------------
619      INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integers
620      !!----------------------------------------------------------------------
621      !
622      IF(lwp) WRITE(numout,*)
623      IF(lwp) WRITE(numout,*) '    zgr_bat_zoom : modify the level bathymetry for zoom domain'
624      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~'
625      !
626      ! Zoom domain
627      ! ===========
628      !
629      ! Forced closed boundary if required
630      IF( lzoom_s )   mbathy(  : , mj0(jpjzoom):mj1(jpjzoom) )      = 0
631      IF( lzoom_w )   mbathy(      mi0(jpizoom):mi1(jpizoom) , :  ) = 0
632      IF( lzoom_e )   mbathy(      mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , :  ) = 0
633      IF( lzoom_n )   mbathy(  : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) )      = 0
634      !
635      ! Configuration specific domain modifications
636      ! (here, ORCA arctic configuration: suppress Med Sea)
637      IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN
638         SELECT CASE ( jp_cfg )
639         !                                        ! =======================
640         CASE ( 2 )                               !  ORCA_R2 configuration
641            !                                     ! =======================
642            IF(lwp) WRITE(numout,*) '                   ORCA R2 arctic zoom: suppress the Med Sea'
643            ii0 = 141   ;   ii1 = 162      ! Sea box i,j indices
644            ij0 =  98   ;   ij1 = 110
645            !                                     ! =======================
646         CASE ( 05 )                              !  ORCA_R05 configuration
647            !                                     ! =======================
648            IF(lwp) WRITE(numout,*) '                   ORCA R05 arctic zoom: suppress the Med Sea'
649            ii0 = 563   ;   ii1 = 642      ! zero over the Med Sea boxe
650            ij0 = 314   ;   ij1 = 370 
651         END SELECT
652         !
653         mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0   ! zero over the Med Sea boxe
654         !
655      ENDIF
656      !
657   END SUBROUTINE zgr_bat_zoom
658
659
660   SUBROUTINE zgr_bat_ctl
661      !!----------------------------------------------------------------------
662      !!                    ***  ROUTINE zgr_bat_ctl  ***
663      !!
664      !! ** Purpose :   check the bathymetry in levels
665      !!
666      !! ** Method  :   The array mbathy is checked to verified its consistency
667      !!      with the model options. in particular:
668      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
669      !!                  along closed boundary.
670      !!            mbathy must be cyclic IF jperio=1.
671      !!            mbathy must be lower or equal to jpk-1.
672      !!            isolated ocean grid points are suppressed from mbathy
673      !!                  since they are only connected to remaining
674      !!                  ocean through vertical diffusion.
675      !!      C A U T I O N : mbathy will be modified during the initializa-
676      !!      tion phase to become the number of non-zero w-levels of a water
677      !!      column, with a minimum value of 1.
678      !!
679      !! ** Action  : - update mbathy: level bathymetry (in level index)
680      !!              - update bathy : meter bathymetry (in meters)
681      !!----------------------------------------------------------------------
682      !!
683      INTEGER ::   ji, jj, jl                    ! dummy loop indices
684      INTEGER ::   icompt, ibtest, ikmax         ! temporary integers
685      REAL(wp), POINTER, DIMENSION(:,:) ::  zbathy
686
687      !!----------------------------------------------------------------------
688      !
689      IF( nn_timing == 1 )  CALL timing_start('zgr_bat_ctl')
690      !
691      CALL wrk_alloc( jpi, jpj, zbathy )
692      !
693      IF(lwp) WRITE(numout,*)
694      IF(lwp) WRITE(numout,*) '    zgr_bat_ctl : check the bathymetry'
695      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
696      !                                          ! Suppress isolated ocean grid points
697      IF(lwp) WRITE(numout,*)
698      IF(lwp) WRITE(numout,*)'                   suppress isolated ocean grid points'
699      IF(lwp) WRITE(numout,*)'                   -----------------------------------'
700      icompt = 0
701      DO jl = 1, 2
702         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
703            mbathy( 1 ,:) = mbathy(jpim1,:)           ! local domain is cyclic east-west
704            mbathy(jpi,:) = mbathy(  2  ,:)
705         ENDIF
706         DO jj = 2, jpjm1
707            DO ji = 2, jpim1
708               ibtest = MAX(  mbathy(ji-1,jj), mbathy(ji+1,jj),   &
709                  &           mbathy(ji,jj-1), mbathy(ji,jj+1)  )
710               IF( ibtest < mbathy(ji,jj) ) THEN
711                  IF(lwp) WRITE(numout,*) ' the number of ocean level at ',   &
712                     &   'grid-point (i,j) =  ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
713                  mbathy(ji,jj) = ibtest
714                  icompt = icompt + 1
715               ENDIF
716            END DO
717         END DO
718      END DO
719      IF( lk_mpp )   CALL mpp_sum( icompt )
720      IF( icompt == 0 ) THEN
721         IF(lwp) WRITE(numout,*)'     no isolated ocean grid points'
722      ELSE
723         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed'
724      ENDIF
725      IF( lk_mpp ) THEN
726         zbathy(:,:) = FLOAT( mbathy(:,:) )
727         CALL lbc_lnk( zbathy, 'T', 1._wp )
728         mbathy(:,:) = INT( zbathy(:,:) )
729      ENDIF
730      !                                          ! East-west cyclic boundary conditions
731      IF( nperio == 0 ) THEN
732         IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio
733         IF( lk_mpp ) THEN
734            IF( nbondi == -1 .OR. nbondi == 2 ) THEN
735               IF( jperio /= 1 )   mbathy(1,:) = 0
736            ENDIF
737            IF( nbondi == 1 .OR. nbondi == 2 ) THEN
738               IF( jperio /= 1 )   mbathy(nlci,:) = 0
739            ENDIF
740         ELSE
741            IF( ln_zco .OR. ln_zps ) THEN
742               mbathy( 1 ,:) = 0
743               mbathy(jpi,:) = 0
744            ELSE
745               mbathy( 1 ,:) = jpkm1
746               mbathy(jpi,:) = jpkm1
747            ENDIF
748         ENDIF
749      ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio ==  6 ) THEN
750         IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio
751         mbathy( 1 ,:) = mbathy(jpim1,:)
752         mbathy(jpi,:) = mbathy(  2  ,:)
753      ELSEIF( nperio == 2 ) THEN
754         IF(lwp) WRITE(numout,*) '   equatorial boundary conditions on mbathy: nperio = ', nperio
755      ELSE
756         IF(lwp) WRITE(numout,*) '    e r r o r'
757         IF(lwp) WRITE(numout,*) '    parameter , nperio = ', nperio
758         !         STOP 'dom_mba'
759      ENDIF
760      !  Boundary condition on mbathy
761      IF( .NOT.lk_mpp ) THEN 
762!!gm     !!bug ???  think about it !
763         !   ... mono- or macro-tasking: T-point, >0, 2D array, no slab
764         zbathy(:,:) = FLOAT( mbathy(:,:) )
765         CALL lbc_lnk( zbathy, 'T', 1._wp )
766         mbathy(:,:) = INT( zbathy(:,:) )
767      ENDIF
768      ! Number of ocean level inferior or equal to jpkm1
769      ikmax = 0
770      DO jj = 1, jpj
771         DO ji = 1, jpi
772            ikmax = MAX( ikmax, mbathy(ji,jj) )
773         END DO
774      END DO
775!!gm  !!! test to do:   ikmax = MAX( mbathy(:,:) )   ???
776      IF( ikmax > jpkm1 ) THEN
777         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' >  jpk-1'
778         IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
779      ELSE IF( ikmax < jpkm1 ) THEN
780         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1' 
781         IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
782      ENDIF
783
784      IF( lwp .AND. nprint == 1 ) THEN      ! control print
785         WRITE(numout,*)
786         WRITE(numout,*) ' bathymetric field :   number of non-zero T-levels '
787         WRITE(numout,*) ' ------------------'
788         CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )
789         WRITE(numout,*)
790      ENDIF
791      !
792      CALL wrk_dealloc( jpi, jpj, zbathy )
793      !
794      IF( nn_timing == 1 )  CALL timing_stop('zgr_bat_ctl')
795      !
796   END SUBROUTINE zgr_bat_ctl
797
798
799   SUBROUTINE zgr_bot_level
800      !!----------------------------------------------------------------------
801      !!                    ***  ROUTINE zgr_bot_level  ***
802      !!
803      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
804      !!
805      !! ** Method  :   computes from mbathy with a minimum value of 1 over land
806      !!
807      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest
808      !!                                     ocean level at t-, u- & v-points
809      !!                                     (min value = 1 over land)
810      !!----------------------------------------------------------------------
811      !!
812      INTEGER ::   ji, jj   ! dummy loop indices
813      REAL(wp), POINTER, DIMENSION(:,:) ::  zmbk
814      !!----------------------------------------------------------------------
815      !
816      IF( nn_timing == 1 )  CALL timing_start('zgr_bot_level')
817      !
818      CALL wrk_alloc( jpi, jpj, zmbk )
819      !
820      IF(lwp) WRITE(numout,*)
821      IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
822      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
823      !
824      mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
825 
826      !                                     ! bottom k-index of W-level = mbkt+1
827      DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level
828         DO ji = 1, jpim1
829            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
830            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
831         END DO
832      END DO
833      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
834      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
835      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
836      !
837      CALL wrk_dealloc( jpi, jpj, zmbk )
838      !
839      IF( nn_timing == 1 )  CALL timing_stop('zgr_bot_level')
840      !
841   END SUBROUTINE zgr_bot_level
842
843      SUBROUTINE zgr_top_level
844      !!----------------------------------------------------------------------
845      !!                    ***  ROUTINE zgr_bot_level  ***
846      !!
847      !! ** Purpose :   defines the vertical index of ocean top (mik. arrays)
848      !!
849      !! ** Method  :   computes from misfdep with a minimum value of 1
850      !!
851      !! ** Action  :   mikt, miku, mikv :   vertical indices of the shallowest
852      !!                                     ocean level at t-, u- & v-points
853      !!                                     (min value = 1)
854      !!----------------------------------------------------------------------
855      !!
856      INTEGER ::   ji, jj   ! dummy loop indices
857      REAL(wp), POINTER, DIMENSION(:,:) ::  zmik
858      !!----------------------------------------------------------------------
859      !
860      IF( nn_timing == 1 )  CALL timing_start('zgr_top_level')
861      !
862      CALL wrk_alloc( jpi, jpj, zmik )
863      !
864      IF(lwp) WRITE(numout,*)
865      IF(lwp) WRITE(numout,*) '    zgr_top_level : ocean top k-index of T-, U-, V- and W-levels '
866      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
867      !
868      mikt(:,:) = MAX( misfdep(:,:) , 1 )    ! top k-index of T-level (=1)
869      !                                      ! top k-index of W-level (=mikt)
870      DO jj = 1, jpjm1                       ! top k-index of U- (U-) level
871         DO ji = 1, jpim1
872            miku(ji,jj) = MAX(  mikt(ji+1,jj  ) , mikt(ji,jj)  )
873            mikv(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj)  )
874            mikf(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj), mikt(ji+1,jj  ), mikt(ji+1,jj+1)  )
875         END DO
876      END DO
877
878      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
879      zmik(:,:) = REAL( miku(:,:), wp )   ;   CALL lbc_lnk(zmik,'U',1.)   ;   miku  (:,:) = MAX( INT( zmik(:,:) ), 1 )
880      zmik(:,:) = REAL( mikv(:,:), wp )   ;   CALL lbc_lnk(zmik,'V',1.)   ;   mikv  (:,:) = MAX( INT( zmik(:,:) ), 1 )
881      zmik(:,:) = REAL( mikf(:,:), wp )   ;   CALL lbc_lnk(zmik,'F',1.)   ;   mikf  (:,:) = MAX( INT( zmik(:,:) ), 1 )
882      !
883      CALL wrk_dealloc( jpi, jpj, zmik )
884      !
885      IF( nn_timing == 1 )  CALL timing_stop('zgr_top_level')
886      !
887   END SUBROUTINE zgr_top_level
888
889   SUBROUTINE zgr_zco
890      !!----------------------------------------------------------------------
891      !!                  ***  ROUTINE zgr_zco  ***
892      !!
893      !! ** Purpose :   define the z-coordinate system
894      !!
895      !! ** Method  :   set 3D coord. arrays to reference 1D array
896      !!----------------------------------------------------------------------
897      INTEGER  ::   jk
898      !!----------------------------------------------------------------------
899      !
900      IF( nn_timing == 1 )  CALL timing_start('zgr_zco')
901      !
902      DO jk = 1, jpk
903         gdept_0 (:,:,jk) = gdept_1d(jk)
904         gdepw_0 (:,:,jk) = gdepw_1d(jk)
905         gdep3w_0(:,:,jk) = gdepw_1d(jk)
906         e3t_0   (:,:,jk) = e3t_1d  (jk)
907         e3u_0   (:,:,jk) = e3t_1d  (jk)
908         e3v_0   (:,:,jk) = e3t_1d  (jk)
909         e3f_0   (:,:,jk) = e3t_1d  (jk)
910         e3w_0   (:,:,jk) = e3w_1d  (jk)
911         e3uw_0  (:,:,jk) = e3w_1d  (jk)
912         e3vw_0  (:,:,jk) = e3w_1d  (jk)
913      END DO
914      !
915      IF( nn_timing == 1 )  CALL timing_stop('zgr_zco')
916      !
917   END SUBROUTINE zgr_zco
918
919
920   SUBROUTINE zgr_zps
921      !!----------------------------------------------------------------------
922      !!                  ***  ROUTINE zgr_zps  ***
923      !!                     
924      !! ** Purpose :   the depth and vertical scale factor in partial step
925      !!      z-coordinate case
926      !!
927      !! ** Method  :   Partial steps : computes the 3D vertical scale factors
928      !!      of T-, U-, V-, W-, UW-, VW and F-points that are associated with
929      !!      a partial step representation of bottom topography.
930      !!
931      !!        The reference depth of model levels is defined from an analytical
932      !!      function the derivative of which gives the reference vertical
933      !!      scale factors.
934      !!        From  depth and scale factors reference, we compute there new value
935      !!      with partial steps  on 3d arrays ( i, j, k ).
936      !!
937      !!              w-level: gdepw_0(i,j,k)  = gdep(k)
938      !!                       e3w_0(i,j,k) = dk(gdep)(k)     = e3(i,j,k)
939      !!              t-level: gdept_0(i,j,k)  = gdep(k+0.5)
940      !!                       e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5)
941      !!
942      !!        With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
943      !!      we find the mbathy index of the depth at each grid point.
944      !!      This leads us to three cases:
945      !!
946      !!              - bathy = 0 => mbathy = 0
947      !!              - 1 < mbathy < jpkm1   
948      !!              - bathy > gdepw_0(jpk) => mbathy = jpkm1 
949      !!
950      !!        Then, for each case, we find the new depth at t- and w- levels
951      !!      and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
952      !!      and f-points.
953      !!
954      !!        This routine is given as an example, it must be modified
955      !!      following the user s desiderata. nevertheless, the output as
956      !!      well as the way to compute the model levels and scale factors
957      !!      must be respected in order to insure second order accuracy
958      !!      schemes.
959      !!
960      !!         c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives
961      !!         - - - - - - -   gdept_0, gdepw_0 and e3. are positives
962      !!     
963      !!  Reference :   Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
964      !!----------------------------------------------------------------------
965      !!
966      INTEGER  ::   ji, jj, jk       ! dummy loop indices
967      INTEGER  ::   ik, it, ikb, ikt ! temporary integers
968      LOGICAL  ::   ll_print         ! Allow  control print for debugging
969      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
970      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
971      REAL(wp) ::   zmax             ! Maximum depth
972      REAL(wp) ::   zdiff            ! temporary scalar
973      REAL(wp) ::   zrefdep          ! temporary scalar
974      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt
975      !!---------------------------------------------------------------------
976      !
977      IF( nn_timing == 1 )  CALL timing_start('zgr_zps')
978      !
979      CALL wrk_alloc( jpi, jpj, jpk, zprt )
980      !
981      IF(lwp) WRITE(numout,*)
982      IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps'
983      IF(lwp) WRITE(numout,*) '    ~~~~~~~ '
984      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used'
985
986      ll_print = .FALSE.                   ! Local variable for debugging
987     
988      IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth
989         WRITE(numout,*)
990         WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)'
991         CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
992      ENDIF
993
994
995      ! bathymetry in level (from bathy_meter)
996      ! ===================
997      zmax = gdepw_1d(jpk) + e3t_1d(jpk)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) )
998      bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat)
999      WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0
1000      ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level
1001      END WHERE
1002
1003      ! Compute mbathy for ocean points (i.e. the number of ocean levels)
1004      ! find the number of ocean levels such that the last level thickness
1005      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where
1006      ! e3t_1d is the reference level thickness
1007      DO jk = jpkm1, 1, -1
1008         zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1009         WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1
1010      END DO
1011
1012      IF ( ln_isfcav ) CALL zgr_isf
1013
1014      ! Scale factors and depth at T- and W-points
1015      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate
1016         gdept_0(:,:,jk) = gdept_1d(jk)
1017         gdepw_0(:,:,jk) = gdepw_1d(jk)
1018         e3t_0  (:,:,jk) = e3t_1d  (jk)
1019         e3w_0  (:,:,jk) = e3w_1d  (jk)
1020      END DO
1021      !
1022      DO jj = 1, jpj
1023         DO ji = 1, jpi
1024            ik = mbathy(ji,jj)
1025            IF( ik > 0 ) THEN               ! ocean point only
1026               ! max ocean level case
1027               IF( ik == jpkm1 ) THEN
1028                  zdepwp = bathy(ji,jj)
1029                  ze3tp  = bathy(ji,jj) - gdepw_1d(ik)
1030                  ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
1031                  e3t_0(ji,jj,ik  ) = ze3tp
1032                  e3t_0(ji,jj,ik+1) = ze3tp
1033                  e3w_0(ji,jj,ik  ) = ze3wp
1034                  e3w_0(ji,jj,ik+1) = ze3tp
1035                  gdepw_0(ji,jj,ik+1) = zdepwp
1036                  gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp
1037                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
1038                  !
1039               ELSE                         ! standard case
1040                  IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
1041                  ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
1042                  ENDIF
1043!gm Bug?  check the gdepw_1d
1044                  !       ... on ik
1045                  gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   &
1046                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   &
1047                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) ))
1048                  e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   & 
1049                     &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ) 
1050                  e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   &
1051                     &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
1052                  !       ... on ik+1
1053                  e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
1054                  e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
1055                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
1056               ENDIF
1057            ENDIF
1058         END DO
1059      END DO
1060      !
1061      it = 0
1062      DO jj = 1, jpj
1063         DO ji = 1, jpi
1064            ik = mbathy(ji,jj)
1065            IF( ik > 0 ) THEN               ! ocean point only
1066               e3tp (ji,jj) = e3t_0(ji,jj,ik)
1067               e3wp (ji,jj) = e3w_0(ji,jj,ik)
1068               ! test
1069               zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  )
1070               IF( zdiff <= 0._wp .AND. lwp ) THEN
1071                  it = it + 1
1072                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj
1073                  WRITE(numout,*) ' bathy = ', bathy(ji,jj)
1074                  WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
1075                  WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  )
1076               ENDIF
1077            ENDIF
1078         END DO
1079      END DO
1080      !
1081      IF ( ln_isfcav ) THEN
1082      ! (ISF) Definition of e3t, u, v, w for ISF case
1083         DO jj = 1, jpj 
1084            DO ji = 1, jpi 
1085               ik = misfdep(ji,jj) 
1086               IF( ik > 1 ) THEN               ! ice shelf point only
1087                  IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik) 
1088                  gdepw_0(ji,jj,ik) = risfdep(ji,jj) 
1089!gm Bug?  check the gdepw_0
1090               !       ... on ik
1091                  gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   & 
1092                     &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   & 
1093                     &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      ) 
1094                  e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 
1095                  e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik)
1096
1097                  IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)
1098                     e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 
1099                  ENDIF 
1100               !       ... on ik / ik-1
1101                  e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 
1102                  e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
1103! The next line isn't required and doesn't affect results - included for consistency with bathymetry code
1104                  gdept_0(ji,jj,ik-1) = gdept_1d(ik-1)
1105               ENDIF
1106            END DO
1107         END DO 
1108      !
1109         it = 0 
1110         DO jj = 1, jpj 
1111            DO ji = 1, jpi 
1112               ik = misfdep(ji,jj) 
1113               IF( ik > 1 ) THEN               ! ice shelf point only
1114                  e3tp (ji,jj) = e3t_0(ji,jj,ik  ) 
1115                  e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 
1116               ! test
1117                  zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  ) 
1118                  IF( zdiff <= 0. .AND. lwp ) THEN 
1119                     it = it + 1 
1120                     WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
1121                     WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 
1122                     WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
1123                     WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj) 
1124                  ENDIF
1125               ENDIF
1126            END DO
1127         END DO
1128      END IF
1129      ! END (ISF)
1130
1131      ! Scale factors and depth at U-, V-, UW and VW-points
1132      DO jk = 1, jpk                        ! initialisation to z-scale factors
1133         e3u_0 (:,:,jk) = e3t_1d(jk)
1134         e3v_0 (:,:,jk) = e3t_1d(jk)
1135         e3uw_0(:,:,jk) = e3w_1d(jk)
1136         e3vw_0(:,:,jk) = e3w_1d(jk)
1137      END DO
1138      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
1139         DO jj = 1, jpjm1
1140            DO ji = 1, fs_jpim1   ! vector opt.
1141               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
1142               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
1143               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
1144               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
1145            END DO
1146         END DO
1147      END DO
1148      IF ( ln_isfcav ) THEN
1149      ! (ISF) define e3uw (adapted for 2 cells in the water column)
1150         DO jj = 2, jpjm1 
1151            DO ji = 2, fs_jpim1   ! vector opt.
1152               ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))
1153               ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))
1154               IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) &
1155                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) )
1156               ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))
1157               ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))
1158               IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) &
1159                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) )
1160            END DO
1161         END DO
1162      END IF
1163
1164      CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions
1165      CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp )
1166      !
1167      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
1168         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk)
1169         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk)
1170         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk)
1171         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk)
1172      END DO
1173     
1174      ! Scale factor at F-point
1175      DO jk = 1, jpk                        ! initialisation to z-scale factors
1176         e3f_0(:,:,jk) = e3t_1d(jk)
1177      END DO
1178      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors
1179         DO jj = 1, jpjm1
1180            DO ji = 1, fs_jpim1   ! vector opt.
1181               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
1182            END DO
1183         END DO
1184      END DO
1185      CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions
1186      !
1187      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
1188         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk)
1189      END DO
1190!!gm  bug ? :  must be a do loop with mj0,mj1
1191      !
1192      e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2
1193      e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 
1194      e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 
1195      e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 
1196      e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 
1197
1198      ! Control of the sign
1199      IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' )
1200      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' )
1201      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' )
1202      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' )
1203     
1204      ! Compute gdep3w_0 (vertical sum of e3w)
1205      IF ( ln_isfcav ) THEN ! if cavity
1206         WHERE (misfdep == 0) misfdep = 1
1207         DO jj = 1,jpj
1208            DO ji = 1,jpi
1209               gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)
1210               DO jk = 2, misfdep(ji,jj)
1211                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
1212               END DO
1213               IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))
1214               DO jk = misfdep(ji,jj) + 1, jpk
1215                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
1216               END DO
1217            END DO
1218         END DO
1219      ELSE ! no cavity
1220         gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
1221         DO jk = 2, jpk
1222            gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk)
1223         END DO
1224      END IF
1225      !                                               ! ================= !
1226      IF(lwp .AND. ll_print) THEN                     !   Control print   !
1227         !                                            ! ================= !
1228         DO jj = 1,jpj
1229            DO ji = 1, jpi
1230               ik = MAX( mbathy(ji,jj), 1 )
1231               zprt(ji,jj,1) = e3t_0   (ji,jj,ik)
1232               zprt(ji,jj,2) = e3w_0   (ji,jj,ik)
1233               zprt(ji,jj,3) = e3u_0   (ji,jj,ik)
1234               zprt(ji,jj,4) = e3v_0   (ji,jj,ik)
1235               zprt(ji,jj,5) = e3f_0   (ji,jj,ik)
1236               zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)
1237            END DO
1238         END DO
1239         WRITE(numout,*)
1240         WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1241         WRITE(numout,*)
1242         WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1243         WRITE(numout,*)
1244         WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1245         WRITE(numout,*)
1246         WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1247         WRITE(numout,*)
1248         WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1249         WRITE(numout,*)
1250         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1251      ENDIF 
1252      !
1253      CALL wrk_dealloc( jpi, jpj, jpk, zprt )
1254      !
1255      IF( nn_timing == 1 )  CALL timing_stop('zgr_zps')
1256      !
1257   END SUBROUTINE zgr_zps
1258
1259   SUBROUTINE zgr_isf
1260      !!----------------------------------------------------------------------
1261      !!                    ***  ROUTINE zgr_isf  ***
1262      !!   
1263      !! ** Purpose :   check the bathymetry in levels
1264      !!   
1265      !! ** Method  :   The water column have to contained at least 2 cells
1266      !!                Bathymetry and isfdraft are modified (dig/close) to respect
1267      !!                this criterion.
1268      !!                          0.0 = definition of minal ice shelf draft
1269      !!                digging and filling base on depth criterion only
1270      !!                          1.0 = set iceshelf to the minimum depth allowed
1271      !!                          1.1 = ground ice shelf if water column less than X m
1272      !!                          1.2 = ensure a minimum thickness for iceshelf cavity
1273      !!                          1.3 = remove channels and single point 'bay'
1274      !!                          1.4 = close single isolated point
1275      !!                compute level
1276      !!                          2.0 = compute level
1277      !!                          2.1 = ensure misfdep is not below bathymetry after step 2.0
1278      !!                digging and filling base on level criterion only               
1279      !!                          3.0 = dig to fit the 2 water level criterion (closed pool possible after this step)
1280      !!                          3.1 = dig enough to ensure connectivity of all the cell beneath an ice shelf
1281      !!                                (most of the close pool should be remove after this step)
1282      !!                          3.2 = fill chimney
1283      !!   
1284      !! ** Action  : - test compatibility between isfdraft and bathy
1285      !!              - bathy and isfdraft are modified
1286      !!----------------------------------------------------------------------
1287      !! 
1288      INTEGER  ::   ji, jj, jk, jn
1289      INTEGER  ::   ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1   ! (ISF)
1290      INTEGER  ::   zmbathyij, zmbathyip1, zmbathyim1, zmbathyjp1, zmbathyjm1
1291      INTEGER , POINTER, DIMENSION(:,:)   ::   zmisfdep, zmbathy         ! 2D workspace (ISH)
1292      !
1293      REAL(wp) ::   zdepth, vmskjp1, vmskjm1, umskip1, umskim1, umskip1_r, umskim1_r, vmskjp1_r, vmskjm1_r
1294      REAL(wp) :: zisfdep_min
1295      REAL(wp), POINTER, DIMENSION(:,:)   ::   zdummy, zmask, zrisfdep   ! 2D workspace (ISH)
1296      !
1297      !!-----------------
1298      ! NAMELIST
1299      INTEGER  ::   ios
1300      INTEGER  :: nn_kisfmax    = 999.                             
1301      REAL(wp) :: rn_isfdep_min = 10.0_wp                           ! ice shelf minimal thickness
1302      REAL(wp) :: rn_glhw_min   = 1.0e-3                             ! threshold on hw to define grounding line water
1303      REAL(wp) :: rn_isfhw_min  = 1.0e-3                            ! threshold on hw to define isf draft into the cavity
1304      REAL(wp) :: rn_isfshallow = 0.0_wp                            ! threshold to define shallow ice shelf cavity
1305      REAL(wp) :: rn_zisfmax    = 6000.0_wp                         ! maximun meter of ice we are allowed to dig to assure connectivity
1306      LOGICAL  :: ln_isfcheminey= .FALSE. , ln_isfconnect= .FALSE. , ln_isfchannel= .FALSE. !  remove cheminey, assure connectivity, remove channel in water column (on z space not level space)
1307      !!---------------------------------------------------------------------
1308      NAMELIST/namzgr_isf/nn_kisfmax, rn_zisfmax, &
1309              &           rn_isfdep_min, rn_glhw_min, rn_isfhw_min, &
1310              &           ln_isfcheminey, ln_isfconnect, ln_isfchannel, ln_isfsubgl, rn_isfsubgllon, rn_isfsubgllat
1311      !
1312      IF( nn_timing == 1 )  CALL timing_start('zgr_isf')
1313      !
1314      CALL wrk_alloc( jpi, jpj, zdummy, zmask, zrisfdep)
1315      CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy )
1316      !
1317      !
1318      REWIND( numnam_ref )              ! Namelist namzgr_isf in reference namelist : ice shelf geometry definition
1319      READ  ( numnam_ref, namzgr_isf, IOSTAT = ios, ERR = 901)
1320901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_isf in reference namelist', lwp )
1321
1322      REWIND( numnam_cfg )              ! Namelist namzgr_sco in configuration namelist : ice shelf geometry definition
1323      READ  ( numnam_cfg, namzgr_isf, IOSTAT = ios, ERR = 902 )
1324902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_isf in configuration namelist', lwp )
1325      IF(lwm) WRITE ( numond, namzgr_isf )
1326      !
1327      IF (lwp) THEN
1328         WRITE(numout,*) ' nn_kisfmax       = ',nn_kisfmax           !
1329         WRITE(numout,*) ' rn_zisfmax       = ',rn_zisfmax           !
1330         WRITE(numout,*) ' rn_isfdep_min    = ',rn_isfdep_min        !
1331         WRITE(numout,*) ' rn_isfhw_min     = ',rn_isfhw_min         !
1332         WRITE(numout,*) ' rn_glhw_min      = ',rn_glhw_min          !
1333         WRITE(numout,*) ' ln_isfcheminey   = ',ln_isfcheminey       !
1334         WRITE(numout,*) ' ln_isfconnect    = ',ln_isfconnect        !
1335         WRITE(numout,*) ' ln_isfchannel    = ',ln_isfchannel        !
1336      END IF
1337      !
1338      ! 0.0 variable definition of minimal ice shelf draft allowed
1339      ! an ice shelf draft have to be larger than the first level thickness
1340      zisfdep_min = MAX(rn_isfdep_min,e3t_1d(1))
1341
1342      ! 1.0 set iceshelf to the minimum depth allowed
1343      WHERE(risfdep(:,:) > 0.0_wp .AND. risfdep(:,:) < zisfdep_min)
1344         risfdep(:,:)=zisfdep_min
1345      END WHERE
1346     
1347      ! 1.1 ground ice shelf if water column less than X m => set the grounding
1348      ! line position
1349      WHERE( bathy(:,:) - risfdep(:,:) < rn_glhw_min ) 
1350         risfdep(:,:)=0.0 ; misfdep(:,:) = 1
1351         bathy  (:,:)=0.0 ; mbathy (:,:) = 0
1352      END WHERE
1353      !
1354      ! 1.2 ensure a minimum thickness for iceshelf cavity => avoid to negative
1355      ! e3t if ssh + sum(e3t*tmask) < 0
1356      DO jj = 1, jpj
1357         DO ji = 1, jpi
1358            IF (bathy(ji,jj)-risfdep(ji,jj) < rn_isfhw_min) THEN
1359               risfdep(ji,jj) = bathy(ji,jj) - rn_isfhw_min
1360               ! sanity check on risfdep (if < zisfdep_min) => we ground it
1361               IF (risfdep(ji,jj) < zisfdep_min) THEN
1362                  risfdep(ji,jj)=0.0 ; misfdep(ji,jj) = 1
1363                  bathy  (ji,jj)=0.0 ; mbathy (ji,jj) = 0
1364               END IF
1365            END IF
1366         END DO
1367      END DO
1368      !
1369      ! 1.3 Remove channels and single point 'bay'.
1370      ! channel could be created if connectivity is enforced later.
1371      IF (ln_isfchannel) THEN
1372         zmask(:,:)=1._wp
1373         WHERE (mbathy(:,:)==0.0)
1374            zmask(:,:)=0._wp
1375         END WHERE
1376         DO jj = 2, jpjm1
1377            DO ji = 2, jpim1
1378               IF(  (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN
1379                  umskip1=zmask(ji,jj)*zmask(ji+1,jj  )  ! 1 = connexion
1380                  vmskjp1=zmask(ji,jj)*zmask(ji  ,jj+1)  ! 1 = connexion
1381                  umskim1=zmask(ji,jj)*zmask(ji-1,jj  )  ! 1 = connexion
1382                  vmskjm1=zmask(ji,jj)*zmask(ji  ,jj-1)  ! 1 = connexion
1383               ! zonal channel and single bay
1384                  IF ((umskip1+umskim1>=1) .AND. (vmskjp1+vmskjm1==0)) THEN
1385                     risfdep(ji,jj)=0.0 ; misfdep(ji,jj) = 1
1386                     bathy  (ji,jj)=0.0 ; mbathy (ji,jj) = 0
1387                  END IF
1388               ! meridionnal channel and single bay
1389                  IF ((vmskjp1+vmskjm1>=1) .AND. (umskip1+umskim1==0)) THEN
1390                     risfdep(ji,jj)=0.0 ; misfdep(ji,jj) = 1
1391                     bathy  (ji,jj)=0.0 ; mbathy (ji,jj) = 0
1392                  END IF
1393               END IF
1394            END DO
1395         END DO
1396         ! ensure halo correct
1397         IF( lk_mpp ) THEN
1398            zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1399            zdummy(:,:) = FLOAT( mbathy(:,:)  ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:)  = INT( zdummy(:,:) )
1400            CALL lbc_lnk( risfdep, 'T', 1. )
1401            CALL lbc_lnk( bathy  , 'T', 1. )
1402         ENDIF       
1403      END IF
1404      !
1405      ! 1.4 Closed single isolated point
1406      zmask(:,:)=1._wp
1407      WHERE (mbathy(:,:)==0.0)
1408         zmask(:,:)=0._wp
1409      END WHERE
1410      DO jj = 2, jpjm1
1411         DO ji = 2, jpim1
1412            IF(  (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN
1413               umskip1=zmask(ji,jj)*zmask(ji+1,jj  )  ! 1 = connexion
1414               vmskjp1=zmask(ji,jj)*zmask(ji  ,jj+1)  ! 1 = connexion
1415               umskim1=zmask(ji,jj)*zmask(ji-1,jj  )  ! 1 = connexion
1416               vmskjm1=zmask(ji,jj)*zmask(ji  ,jj-1)  ! 1 = connexion
1417               ! remove single point
1418               IF (umskip1+umskim1+vmskjp1+vmskjm1==0.0_wp) THEN
1419                  risfdep(ji,jj)=0.0 ; misfdep(ji,jj) = 1
1420                  bathy  (ji,jj)=0.0 ; mbathy (ji,jj) = 0
1421               END IF
1422            END IF
1423         END DO
1424      END DO
1425      ! ensure halo correct
1426      IF( lk_mpp ) THEN
1427         zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1428         zdummy(:,:) = FLOAT( mbathy(:,:)  ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:)  = INT( zdummy(:,:) )
1429         CALL lbc_lnk( risfdep, 'T', 1. )
1430         CALL lbc_lnk( bathy  , 'T', 1. )
1431      ENDIF
1432      !
1433      ! 2 Compute misfdep for ocean points (i.e. first wet level)
1434      ! find the first ocean level such that the first level thickness
1435      ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where
1436      ! e3t_0 is the reference level thickness
1437      !
1438      ! 2.0 set misfdep to 1 on open water and initialisation beneath ice shelf
1439      WHERE( risfdep(:,:) == 0.0 ) ;   misfdep(:,:) = 1   ! open water or grounded : set misfdep to 1 
1440      ELSEWHERE                    ;   misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level
1441      END WHERE 
1442      !
1443      DO jk = 2, jpkm1 
1444         zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
1445         WHERE( risfdep(:,:) > 0.0 .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1 
1446      END DO 
1447      !
1448      ! 2.1 be sure misfdep not below mbathy
1449      ! warning because of condition 4 we could have 'wet cell with misfdep below mbathy
1450      ! risfdep of these cells will be fix later on (see 5)
1451      WHERE( misfdep > mbathy .AND. mbathy > 0 ) misfdep=mbathy
1452      !
1453      ! 3.0 Assure 2 wet cells in the column at T point and along the edge.
1454      ! first do T, then loop 4 times on U-1, U, V-1, V and then at T point
1455      ! If un-luky, digging cell U-1 can trigger case for U+1, then V-1, then V+1
1456      !
1457      ! find the deepest isfdep level that fit the 2 wet cell on the water column
1458      ! on all the sides (still need 4 pass)
1459      DO jn = 1, 4
1460         zmisfdep = misfdep
1461         DO jj = 2, jpjm1
1462            DO ji = 2, jpim1
1463               ! ISF cell only
1464               IF(  (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN
1465                  zmbathyij  = mbathy(ji  ,jj)
1466                  ! set ground edge value to jpk to skip it later on
1467                  zmbathyip1 = mbathy(ji+1,jj) ; IF (zmbathyip1 < misfdep(ji,jj)) zmbathyip1 = jpk ! not wet cell in common on this edge
1468                  zmbathyim1 = mbathy(ji-1,jj) ; IF (zmbathyim1 < misfdep(ji,jj)) zmbathyim1 = jpk ! not wet cell in common on this edge
1469                  zmbathyjp1 = mbathy(ji,jj+1) ; IF (zmbathyjp1 < misfdep(ji,jj)) zmbathyjp1 = jpk ! not wet cell in common on this edge
1470                  zmbathyjm1 = mbathy(ji,jj-1) ; IF (zmbathyjm1 < misfdep(ji,jj)) zmbathyjm1 = jpk ! not wet cell in common on this edge
1471                  ! shallowest bathy level
1472                  zmbathyij = MIN(zmbathyij,zmbathyip1,zmbathyim1,zmbathyjp1,zmbathyjm1)
1473                  ! misfdep need to be <= zmbathyij-1 to fit 2 wet cell on the
1474                  ! water column
1475                  jk = MIN(misfdep(ji,jj),zmbathyij-1)
1476                  ! update isf draft if needed
1477                  IF ( jk < misfdep(ji,jj) ) THEN
1478                     zmisfdep(ji,jj) = jk
1479                     risfdep(ji,jj)  = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1480                  ELSE
1481                     ! isf depth not modify, nothing to do
1482                  END IF
1483               ENDIF
1484            END DO
1485         END DO
1486         misfdep=zmisfdep
1487         ! ensure halo correct before new pass
1488         IF( lk_mpp ) THEN
1489            zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1490            zdummy(:,:) = FLOAT( mbathy(:,:)  ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:)  = INT( zdummy(:,:) )
1491            CALL lbc_lnk( risfdep, 'T', 1. )
1492            CALL lbc_lnk( bathy  , 'T', 1. )
1493         ENDIF
1494      END DO ! jn
1495      !
1496      ! 3.1 condition block to assure connectivity every where beneath an ice shelf
1497      IF (ln_isfconnect) THEN
1498         zmask(:,:)=1.0
1499         WHERE (mbathy==0)
1500            zmask(:,:)=jpk
1501         END WHERE
1502
1503         zmbathy  = mbathy
1504         zmisfdep = misfdep
1505         zrisfdep = risfdep
1506         WHERE (mbathy == 0) zmbathy=jpk
1507         DO jj = 2, jpjm1
1508            DO ji = 2, jpim1
1509               IF(  (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN
1510! what it should be
1511                  umskip1=zmask(ji,jj)*zmask(ji+1,jj  )  ! 1 = should be connected
1512                  umskim1=zmask(ji,jj)*zmask(ji-1,jj  )  ! 1 = should be connected
1513                  vmskjp1=zmask(ji,jj)*zmask(ji  ,jj+1)  ! 1 = should be connected
1514                  vmskjm1=zmask(ji,jj)*zmask(ji  ,jj-1)  ! 1 = should be connected
1515! what it is
1516                  umskip1_r=jpk ; umskim1_r=jpk; vmskjp1_r=jpk; vmskjm1_r=jpk
1517                  IF (misfdep(ji,jj) > zmbathy(ji+1,jj  )) umskip1_r=1.0 ! 1 = no effective connection
1518                  IF (misfdep(ji,jj) > zmbathy(ji-1,jj  )) umskim1_r=1.0
1519                  IF (misfdep(ji,jj) > zmbathy(ji  ,jj+1)) vmskjp1_r=1.0
1520                  IF (misfdep(ji,jj) > zmbathy(ji  ,jj-1)) vmskjm1_r=1.0
1521                  ! defining level need for connectivity
1522                  jk=NINT(MIN(zmbathy(ji+1,jj  ) * umskip1_r * umskip1, &
1523                     &        zmbathy(ji-1,jj  ) * umskim1_r * umskim1, &
1524                     &        zmbathy(ji  ,jj+1) * vmskjp1_r * vmskjp1, &
1525                     &        zmbathy(ji  ,jj-1) * vmskjm1_r * vmskjm1, &
1526                     &        jpk+1.0 )) ! add jpk in the MIN to avoid out of boundary later on
1527! if connectivity is OK or no connection needed (grounding line) or grounded, zmisfdep=misfdep
1528                  zmisfdep(ji,jj)=MIN(misfdep(ji,jj),jk-1)
1529                  !
1530                  ! update isf draft if needed (need to be done here because next test is in meter and level)
1531                  IF ( zmisfdep(ji,jj) < misfdep(ji,jj) ) THEN
1532                     jk = zmisfdep(ji,jj)
1533                     zrisfdep(ji,jj)  = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1534                  END IF
1535   
1536                  ! sanity check
1537                  ! check if we dig more than nn_kisfmax level or reach the surface
1538                  ! check if we dig more than rn_zisfmax meter
1539                  ! => if this is the case, undo what has been done before
1540                  IF (      (misfdep(ji,jj)-zmisfdep(ji,jj) > MIN(nn_kisfmax,misfdep(ji,jj)-2)) &
1541                     & .OR. (risfdep(ji,jj)-zrisfdep(ji,jj) > MIN(rn_zisfmax,risfdep(ji,jj)  )) ) THEN
1542                     zmisfdep(ji,jj)=misfdep(ji,jj) 
1543                     zrisfdep(ji,jj)=risfdep(ji,jj)
1544                  END IF
1545               END IF
1546            END DO
1547         END DO
1548         misfdep=zmisfdep
1549         risfdep=zrisfdep
1550         ! ensure halo correct
1551         IF( lk_mpp ) THEN
1552            zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1553            zdummy(:,:) = FLOAT( mbathy(:,:)  ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:)  = INT( zdummy(:,:) )
1554            CALL lbc_lnk( risfdep, 'T', 1. )
1555            CALL lbc_lnk( bathy  , 'T', 1. )
1556         ENDIF
1557      END IF
1558      !
1559      IF (ln_isfcheminey) THEN
1560      ! 3.2 fill hole in ice shelf (ie cell with no velocity point)
1561      !      => misfdep = MIN(misfdep at U, U-1, V, V-1)
1562      !         risfdep = gdepw(misfdep) (ie full cell)
1563         zmisfdep = misfdep
1564         WHERE (mbathy == 0) zmisfdep=jpk   ! grounded
1565         DO jj = 2, jpjm1
1566            DO ji = 2, jpim1
1567               ibtest    = zmisfdep(ji  ,jj  )
1568               ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  )
1569               ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1)
1570               ! case unconected wet cell (T point) where isfdraft(ii,jj) deeper
1571               ! than bathy U-1/U+1/V-1/V+1
1572               IF( zmisfdep(ji,jj) > mbathy(ji-1,jj  ) ) ibtestim1 = jpk ! mask at U-1
1573               IF( zmisfdep(ji,jj) > mbathy(ji+1,jj  ) ) ibtestip1 = jpk ! mask at U+1
1574               IF( zmisfdep(ji,jj) > mbathy(ji  ,jj-1) ) ibtestjm1 = jpk ! mask at V-1
1575               IF( zmisfdep(ji,jj) > mbathy(ji  ,jj+1) ) ibtestjp1 = jpk ! mask at V+1
1576               ! case unconected wet cell (T point) where bathy(ii,jj) shallower
1577               ! than isfdraft U-1/U+1/V-1/V+1
1578               IF( mbathy(ji,jj) < zmisfdep(ji-1,jj  ) ) ibtestim1 = jpk ! mask at U-1
1579               IF( mbathy(ji,jj) < zmisfdep(ji+1,jj  ) ) ibtestip1 = jpk ! mask at U+1
1580               IF( mbathy(ji,jj) < zmisfdep(ji  ,jj-1) ) ibtestjm1 = jpk ! mask at V-1
1581               IF( mbathy(ji,jj) < zmisfdep(ji  ,jj+1) ) ibtestjp1 = jpk ! mask at V+1
1582               ! if surround in fully mask => mask cell  (1)
1583               ! if hole in the ice shelf, set to the min of surrounding (the MIN is doing the job)
1584               ! if misfdep is not changed, nothing is done (the MAX is doing the
1585               ! job)
1586               ibtest = MAX(ibtest,MIN(ibtestim1, ibtestip1, ibtestjm1,ibtestjp1))
1587               IF (misfdep(ji,jj) < ibtest) THEN
1588                  misfdep(ji,jj) = ibtest
1589                  risfdep(ji,jj) = gdepw_1d(ibtest)
1590               END IF
1591            ENDDO
1592         ENDDO
1593         WHERE( misfdep == jpk)   ! ground case (1)
1594            mbathy = 0 ; bathy = 0.0_wp ; misfdep = 1 ; risfdep = 0.0_wp
1595         END WHERE
1596         !
1597         ! ensure halo correct
1598         IF( lk_mpp ) THEN
1599            zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1600            zdummy(:,:) = FLOAT( mbathy(:,:)  ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:)  = INT( zdummy(:,:) )
1601            CALL lbc_lnk( risfdep, 'T', 1. )
1602            CALL lbc_lnk( bathy  , 'T', 1. )
1603         END IF
1604      END IF
1605      !
1606      CALL wrk_dealloc( jpi, jpj, zmask, zdummy, zrisfdep)
1607      CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy)
1608      !
1609      IF( nn_timing == 1 )  CALL timing_stop('zgr_isf')
1610     
1611   END SUBROUTINE
1612
1613   SUBROUTINE zgr_sco
1614      !!----------------------------------------------------------------------
1615      !!                  ***  ROUTINE zgr_sco  ***
1616      !!                     
1617      !! ** Purpose :   define the s-coordinate system
1618      !!
1619      !! ** Method  :   s-coordinate
1620      !!         The depth of model levels is defined as the product of an
1621      !!      analytical function by the local bathymetry, while the vertical
1622      !!      scale factors are defined as the product of the first derivative
1623      !!      of the analytical function by the bathymetry.
1624      !!      (this solution save memory as depth and scale factors are not
1625      !!      3d fields)
1626      !!          - Read bathymetry (in meters) at t-point and compute the
1627      !!         bathymetry at u-, v-, and f-points.
1628      !!            hbatu = mi( hbatt )
1629      !!            hbatv = mj( hbatt )
1630      !!            hbatf = mi( mj( hbatt ) )
1631      !!          - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
1632      !!         function and its derivative given as function.
1633      !!            z_gsigt(k) = fssig (k    )
1634      !!            z_gsigw(k) = fssig (k-0.5)
1635      !!            z_esigt(k) = fsdsig(k    )
1636      !!            z_esigw(k) = fsdsig(k-0.5)
1637      !!      Three options for stretching are give, and they can be modified
1638      !!      following the users requirements. Nevertheless, the output as
1639      !!      well as the way to compute the model levels and scale factors
1640      !!      must be respected in order to insure second order accuracy
1641      !!      schemes.
1642      !!
1643      !!      The three methods for stretching available are:
1644      !!
1645      !!           s_sh94 (Song and Haidvogel 1994)
1646      !!                a sinh/tanh function that allows sigma and stretched sigma
1647      !!
1648      !!           s_sf12 (Siddorn and Furner 2012?)
1649      !!                allows the maintenance of fixed surface and or
1650      !!                bottom cell resolutions (cf. geopotential coordinates)
1651      !!                within an analytically derived stretched S-coordinate framework.
1652      !!
1653      !!          s_tanh  (Madec et al 1996)
1654      !!                a cosh/tanh function that gives stretched coordinates       
1655      !!
1656      !!----------------------------------------------------------------------
1657      !
1658      INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
1659      INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
1660      INTEGER  ::   ios                      ! Local integer output status for namelist read
1661      REAL(wp) ::   zrmax, ztaper   ! temporary scalars
1662      REAL(wp) ::   zrfact
1663      !
1664      REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
1665      REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
1666
1667      NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
1668                           rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
1669     !!----------------------------------------------------------------------
1670      !
1671      IF( nn_timing == 1 )  CALL timing_start('zgr_sco')
1672      !
1673      CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
1674      !
1675      REWIND( numnam_ref )              ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
1676      READ  ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901)
1677901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist', lwp )
1678
1679      REWIND( numnam_cfg )              ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters
1680      READ  ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 )
1681902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist', lwp )
1682      IF(lwm) WRITE ( numond, namzgr_sco )
1683
1684      IF(lwp) THEN                           ! control print
1685         WRITE(numout,*)
1686         WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
1687         WRITE(numout,*) '~~~~~~~~~~~'
1688         WRITE(numout,*) '   Namelist namzgr_sco'
1689         WRITE(numout,*) '     stretching coeffs '
1690         WRITE(numout,*) '        maximum depth of s-bottom surface (>0)       rn_sbot_max   = ',rn_sbot_max
1691         WRITE(numout,*) '        minimum depth of s-bottom surface (>0)       rn_sbot_min   = ',rn_sbot_min
1692         WRITE(numout,*) '        Critical depth                               rn_hc         = ',rn_hc
1693         WRITE(numout,*) '        maximum cut-off r-value allowed              rn_rmax       = ',rn_rmax
1694         WRITE(numout,*) '     Song and Haidvogel 1994 stretching              ln_s_sh94     = ',ln_s_sh94
1695         WRITE(numout,*) '        Song and Haidvogel 1994 stretching coefficients'
1696         WRITE(numout,*) '        surface control parameter (0<=rn_theta<=20)  rn_theta      = ',rn_theta
1697         WRITE(numout,*) '        bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ',rn_thetb
1698         WRITE(numout,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb
1699         WRITE(numout,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12
1700         WRITE(numout,*) '        switching to sigma (T) or Z (F) at H<Hc      ln_sigcrit    = ',ln_sigcrit
1701         WRITE(numout,*) '        Siddorn and Furner 2012 stretching coefficients'
1702         WRITE(numout,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha
1703         WRITE(numout,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold
1704         WRITE(numout,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs
1705         WRITE(numout,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a
1706         WRITE(numout,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b
1707         WRITE(numout,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b'
1708      ENDIF
1709
1710      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate
1711      hifu(:,:) = rn_sbot_min
1712      hifv(:,:) = rn_sbot_min
1713      hiff(:,:) = rn_sbot_min
1714
1715      !                                        ! set maximum ocean depth
1716      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1717
1718      DO jj = 1, jpj
1719         DO ji = 1, jpi
1720           IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1721         END DO
1722      END DO
1723      !                                        ! =============================
1724      !                                        ! Define the envelop bathymetry   (hbatt)
1725      !                                        ! =============================
1726      ! use r-value to create hybrid coordinates
1727      zenv(:,:) = bathy(:,:)
1728      !
1729      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
1730      DO jj = 1, jpj
1731         DO ji = 1, jpi
1732           IF( bathy(ji,jj) == 0._wp ) THEN
1733             iip1 = MIN( ji+1, jpi )
1734             ijp1 = MIN( jj+1, jpj )
1735             iim1 = MAX( ji-1, 1 )
1736             ijm1 = MAX( jj-1, 1 )
1737             IF( ( + bathy(iim1,ijp1) + bathy(ji,ijp1) + bathy(iip1,ijp1)  &
1738                &  + bathy(iim1,jj  )                  + bathy(iip1,jj  )  &
1739                &  + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1)  ) > 0._wp ) THEN
1740                zenv(ji,jj) = rn_sbot_min
1741             ENDIF
1742           ENDIF
1743         END DO
1744      END DO
1745      ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1746      CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1747      !
1748      ! smooth the bathymetry (if required)
1749      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
1750      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
1751      !
1752      jl = 0
1753      zrmax = 1._wp
1754      !   
1755      !     
1756      ! set scaling factor used in reducing vertical gradients
1757      zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
1758      !
1759      ! initialise temporary evelope depth arrays
1760      ztmpi1(:,:) = zenv(:,:)
1761      ztmpi2(:,:) = zenv(:,:)
1762      ztmpj1(:,:) = zenv(:,:)
1763      ztmpj2(:,:) = zenv(:,:)
1764      !
1765      ! initialise temporary r-value arrays
1766      zri(:,:) = 1._wp
1767      zrj(:,:) = 1._wp
1768      !                                                            ! ================ !
1769      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) !  Iterative loop  !
1770         !                                                         ! ================ !
1771         jl = jl + 1
1772         zrmax = 0._wp
1773         ! we set zrmax from previous r-values (zri and zrj) first
1774         ! if set after current r-value calculation (as previously)
1775         ! we could exit DO WHILE prematurely before checking r-value
1776         ! of current zenv
1777         DO jj = 1, nlcj
1778            DO ji = 1, nlci
1779               zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
1780            END DO
1781         END DO
1782         zri(:,:) = 0._wp
1783         zrj(:,:) = 0._wp
1784         DO jj = 1, nlcj
1785            DO ji = 1, nlci
1786               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
1787               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
1788               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
1789                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
1790               END IF
1791               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
1792                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
1793               END IF
1794               IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact
1795               IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact
1796               IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact
1797               IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact
1798            END DO
1799         END DO
1800         IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
1801         !
1802         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax
1803         !
1804         DO jj = 1, nlcj
1805            DO ji = 1, nlci
1806               zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
1807            END DO
1808         END DO
1809         ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1810         CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1811         !                                                  ! ================ !
1812      END DO                                                !     End loop     !
1813      !                                                     ! ================ !
1814      DO jj = 1, jpj
1815         DO ji = 1, jpi
1816            zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
1817         END DO
1818      END DO
1819      !
1820      ! Envelope bathymetry saved in hbatt
1821      hbatt(:,:) = zenv(:,:) 
1822      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
1823         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1824         DO jj = 1, jpj
1825            DO ji = 1, jpi
1826               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
1827               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
1828            END DO
1829         END DO
1830      ENDIF
1831      !
1832      IF(lwp) THEN                             ! Control print
1833         WRITE(numout,*)
1834         WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
1835         WRITE(numout,*)
1836         CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )
1837         IF( nprint == 1 )   THEN       
1838            WRITE(numout,*) ' bathy  MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
1839            WRITE(numout,*) ' hbatt  MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
1840         ENDIF
1841      ENDIF
1842
1843      !                                        ! ==============================
1844      !                                        !   hbatu, hbatv, hbatf fields
1845      !                                        ! ==============================
1846      IF(lwp) THEN
1847         WRITE(numout,*)
1848         WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
1849      ENDIF
1850      hbatu(:,:) = rn_sbot_min
1851      hbatv(:,:) = rn_sbot_min
1852      hbatf(:,:) = rn_sbot_min
1853      DO jj = 1, jpjm1
1854        DO ji = 1, jpim1   ! NO vector opt.
1855           hbatu(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
1856           hbatv(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
1857           hbatf(ji,jj) = 0.25_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
1858              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
1859        END DO
1860      END DO
1861      !
1862      ! Apply lateral boundary condition
1863!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
1864      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk( hbatu, 'U', 1._wp )
1865      DO jj = 1, jpj
1866         DO ji = 1, jpi
1867            IF( hbatu(ji,jj) == 0._wp ) THEN
1868               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
1869               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
1870            ENDIF
1871         END DO
1872      END DO
1873      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk( hbatv, 'V', 1._wp )
1874      DO jj = 1, jpj
1875         DO ji = 1, jpi
1876            IF( hbatv(ji,jj) == 0._wp ) THEN
1877               IF( zhbat(ji,jj) == 0._wp )   hbatv(ji,jj) = rn_sbot_min
1878               IF( zhbat(ji,jj) /= 0._wp )   hbatv(ji,jj) = zhbat(ji,jj)
1879            ENDIF
1880         END DO
1881      END DO
1882      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk( hbatf, 'F', 1._wp )
1883      DO jj = 1, jpj
1884         DO ji = 1, jpi
1885            IF( hbatf(ji,jj) == 0._wp ) THEN
1886               IF( zhbat(ji,jj) == 0._wp )   hbatf(ji,jj) = rn_sbot_min
1887               IF( zhbat(ji,jj) /= 0._wp )   hbatf(ji,jj) = zhbat(ji,jj)
1888            ENDIF
1889         END DO
1890      END DO
1891
1892!!bug:  key_helsinki a verifer
1893      hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1894      hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1895      hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1896      hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1897
1898      IF( nprint == 1 .AND. lwp )   THEN
1899         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
1900            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
1901         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
1902            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
1903         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
1904            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1905         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
1906            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1907      ENDIF
1908!! helsinki
1909
1910      !                                            ! =======================
1911      !                                            !   s-ccordinate fields     (gdep., e3.)
1912      !                                            ! =======================
1913      !
1914      ! non-dimensional "sigma" for model level depth at w- and t-levels
1915
1916
1917!========================================================================
1918! Song and Haidvogel  1994 (ln_s_sh94=T)
1919! Siddorn and Furner 2012 (ln_sf12=T)
1920! or  tanh function       (both false)                   
1921!========================================================================
1922      IF      ( ln_s_sh94 ) THEN
1923                           CALL s_sh94()
1924      ELSE IF ( ln_s_sf12 ) THEN
1925                           CALL s_sf12()
1926      ELSE                 
1927                           CALL s_tanh()
1928      ENDIF
1929
1930      CALL lbc_lnk( e3t_0 , 'T', 1._wp )
1931      CALL lbc_lnk( e3u_0 , 'U', 1._wp )
1932      CALL lbc_lnk( e3v_0 , 'V', 1._wp )
1933      CALL lbc_lnk( e3f_0 , 'F', 1._wp )
1934      CALL lbc_lnk( e3w_0 , 'W', 1._wp )
1935      CALL lbc_lnk( e3uw_0, 'U', 1._wp )
1936      CALL lbc_lnk( e3vw_0, 'V', 1._wp )
1937
1938      fsdepw(:,:,:) = gdepw_0 (:,:,:)
1939      fsde3w(:,:,:) = gdep3w_0(:,:,:)
1940      !
1941      where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0
1942      where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0
1943      where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0
1944      where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0
1945      where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0
1946      where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0
1947      where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0
1948
1949#if defined key_agrif
1950      ! Ensure meaningful vertical scale factors in ghost lines/columns
1951      IF( .NOT. Agrif_Root() ) THEN
1952         
1953         IF((nbondi == -1).OR.(nbondi == 2)) THEN
1954            e3u_0(1,:,:) = e3u_0(2,:,:)
1955         ENDIF
1956         !
1957         IF((nbondi ==  1).OR.(nbondi == 2)) THEN
1958            e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:)
1959         ENDIF
1960         !
1961         IF((nbondj == -1).OR.(nbondj == 2)) THEN
1962            e3v_0(:,1,:) = e3v_0(:,2,:)
1963         ENDIF
1964         !
1965         IF((nbondj ==  1).OR.(nbondj == 2)) THEN
1966            e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:)
1967         ENDIF
1968         !
1969      ENDIF
1970#endif
1971
1972      fsdept(:,:,:) = gdept_0 (:,:,:)
1973      fsdepw(:,:,:) = gdepw_0 (:,:,:)
1974      fsde3w(:,:,:) = gdep3w_0(:,:,:)
1975      fse3t (:,:,:) = e3t_0   (:,:,:)
1976      fse3u (:,:,:) = e3u_0   (:,:,:)
1977      fse3v (:,:,:) = e3v_0   (:,:,:)
1978      fse3f (:,:,:) = e3f_0   (:,:,:)
1979      fse3w (:,:,:) = e3w_0   (:,:,:)
1980      fse3uw(:,:,:) = e3uw_0  (:,:,:)
1981      fse3vw(:,:,:) = e3vw_0  (:,:,:)
1982!!
1983      ! HYBRID :
1984      DO jj = 1, jpj
1985         DO ji = 1, jpi
1986            DO jk = 1, jpkm1
1987               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
1988            END DO
1989            IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0
1990         END DO
1991      END DO
1992      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
1993         &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
1994
1995      IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain
1996         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) )
1997         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
1998            &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w '  , MINVAL( gdep3w_0(:,:,:) )
1999         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0   (:,:,:) ),   &
2000            &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0   (:,:,:) ),   &
2001            &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0  (:,:,:) ),   &
2002            &                          ' w ', MINVAL( e3w_0  (:,:,:) )
2003
2004         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
2005            &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w '  , MAXVAL( gdep3w_0(:,:,:) )
2006         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0   (:,:,:) ),   &
2007            &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0   (:,:,:) ),   &
2008            &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0  (:,:,:) ),   &
2009            &                          ' w ', MAXVAL( e3w_0  (:,:,:) )
2010      ENDIF
2011      !  END DO
2012      IF(lwp) THEN                                  ! selected vertical profiles
2013         WRITE(numout,*)
2014         WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
2015         WRITE(numout,*) ' ~~~~~~  --------------------'
2016         WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2017         WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     &
2018            &                                 e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
2019         DO jj = mj0(20), mj1(20)
2020            DO ji = mi0(20), mi1(20)
2021               WRITE(numout,*)
2022               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
2023               WRITE(numout,*) ' ~~~~~~  --------------------'
2024               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2025               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
2026                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2027            END DO
2028         END DO
2029         DO jj = mj0(74), mj1(74)
2030            DO ji = mi0(100), mi1(100)
2031               WRITE(numout,*)
2032               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
2033               WRITE(numout,*) ' ~~~~~~  --------------------'
2034               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2035               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
2036                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2037            END DO
2038         END DO
2039      ENDIF
2040
2041!================================================================================
2042! check the coordinate makes sense
2043!================================================================================
2044      DO ji = 1, jpi
2045         DO jj = 1, jpj
2046
2047            IF( hbatt(ji,jj) > 0._wp) THEN
2048               DO jk = 1, mbathy(ji,jj)
2049                 ! check coordinate is monotonically increasing
2050                 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
2051                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
2052                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
2053                    WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
2054                    WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
2055                    CALL ctl_stop( ctmp1 )
2056                 ENDIF
2057                 ! and check it has never gone negative
2058                 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
2059                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
2060                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk
2061                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2062                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2063                    CALL ctl_stop( ctmp1 )
2064                 ENDIF
2065                 ! and check it never exceeds the total depth
2066                 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
2067                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
2068                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
2069                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2070                    CALL ctl_stop( ctmp1 )
2071                 ENDIF
2072               END DO
2073
2074               DO jk = 1, mbathy(ji,jj)-1
2075                 ! and check it never exceeds the total depth
2076                IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
2077                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
2078                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
2079                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2080                    CALL ctl_stop( ctmp1 )
2081                 ENDIF
2082               END DO
2083
2084            ENDIF
2085
2086         END DO
2087      END DO
2088      !
2089      CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
2090      !
2091      IF( nn_timing == 1 )  CALL timing_stop('zgr_sco')
2092      !
2093   END SUBROUTINE zgr_sco
2094
2095!!======================================================================
2096   SUBROUTINE s_sh94()
2097
2098      !!----------------------------------------------------------------------
2099      !!                  ***  ROUTINE s_sh94  ***
2100      !!                     
2101      !! ** Purpose :   stretch the s-coordinate system
2102      !!
2103      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994
2104      !!                mixed S/sigma coordinate
2105      !!
2106      !! Reference : Song and Haidvogel 1994.
2107      !!----------------------------------------------------------------------
2108      !
2109      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2110      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
2111      !
2112      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2113      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
2114
2115      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2116      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2117
2118      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
2119      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
2120      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
2121      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
2122
2123      DO ji = 1, jpi
2124         DO jj = 1, jpj
2125
2126            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
2127               DO jk = 1, jpk
2128                  z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
2129                  z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
2130               END DO
2131            ELSE ! shallow water, uniform sigma
2132               DO jk = 1, jpk
2133                  z_gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
2134                  z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
2135                  END DO
2136            ENDIF
2137            !
2138            DO jk = 1, jpkm1
2139               z_esigt3(ji,jj,jk  ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2140               z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2141            END DO
2142            z_esigw3(ji,jj,1  ) = 2._wp * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) )
2143            z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
2144            !
2145            ! Coefficients for vertical depth as the sum of e3w scale factors
2146            z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
2147            DO jk = 2, jpk
2148               z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2149            END DO
2150            !
2151            DO jk = 1, jpk
2152               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2153               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2154               gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
2155               gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
2156               gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
2157            END DO
2158           !
2159         END DO   ! for all jj's
2160      END DO    ! for all ji's
2161
2162      DO ji = 1, jpim1
2163         DO jj = 1, jpjm1
2164            DO jk = 1, jpk
2165               z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
2166                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2167               z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
2168                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2169               z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     &
2170                  &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   &
2171                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2172               z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
2173                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2174               z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
2175                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2176               !
2177               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2178               e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2179               e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2180               e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2181               !
2182               e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2183               e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2184               e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2185            END DO
2186        END DO
2187      END DO
2188
2189      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2190      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2191
2192   END SUBROUTINE s_sh94
2193
2194   SUBROUTINE s_sf12
2195
2196      !!----------------------------------------------------------------------
2197      !!                  ***  ROUTINE s_sf12 ***
2198      !!                     
2199      !! ** Purpose :   stretch the s-coordinate system
2200      !!
2201      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
2202      !!                mixed S/sigma/Z coordinate
2203      !!
2204      !!                This method allows the maintenance of fixed surface and or
2205      !!                bottom cell resolutions (cf. geopotential coordinates)
2206      !!                within an analytically derived stretched S-coordinate framework.
2207      !!
2208      !!
2209      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
2210      !!----------------------------------------------------------------------
2211      !
2212      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2213      REAL(wp) ::   zsmth               ! smoothing around critical depth
2214      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
2215      !
2216      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2217      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
2218
2219      !
2220      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2221      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2222
2223      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
2224      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
2225      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
2226      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
2227
2228      DO ji = 1, jpi
2229         DO jj = 1, jpj
2230
2231          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
2232             
2233              zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
2234                                                     ! could be changed by users but care must be taken to do so carefully
2235              zzb = 1.0_wp-(zzb/hbatt(ji,jj))
2236           
2237              zzs = rn_zs / hbatt(ji,jj) 
2238             
2239              IF (rn_efold /= 0.0_wp) THEN
2240                zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
2241              ELSE
2242                zsmth = 1.0_wp 
2243              ENDIF
2244               
2245              DO jk = 1, jpk
2246                z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
2247                z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
2248              ENDDO
2249              z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth  )
2250              z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth  )
2251 
2252          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
2253
2254            DO jk = 1, jpk
2255              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
2256              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
2257            END DO
2258
2259          ELSE  ! shallow water, z coordinates
2260
2261            DO jk = 1, jpk
2262              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2263              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2264            END DO
2265
2266          ENDIF
2267
2268          DO jk = 1, jpkm1
2269             z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2270             z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2271          END DO
2272          z_esigw3(ji,jj,1  ) = 2.0_wp * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ))
2273          z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
2274
2275          ! Coefficients for vertical depth as the sum of e3w scale factors
2276          z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
2277          DO jk = 2, jpk
2278             z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2279          END DO
2280
2281          DO jk = 1, jpk
2282             gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
2283             gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
2284             gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
2285          END DO
2286
2287        ENDDO   ! for all jj's
2288      ENDDO    ! for all ji's
2289
2290      DO ji=1,jpi-1
2291        DO jj=1,jpj-1
2292
2293          DO jk = 1, jpk
2294                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / &
2295                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2296                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / &
2297                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2298                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  &
2299                                      hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / &
2300                                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2301                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / &
2302                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2303                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / &
2304                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2305
2306             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
2307             e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
2308             e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
2309             e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
2310             !
2311             e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
2312             e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
2313             e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
2314          END DO
2315
2316        ENDDO
2317      ENDDO
2318      !
2319      CALL lbc_lnk(e3t_0 ,'T',1.) ; CALL lbc_lnk(e3u_0 ,'T',1.)
2320      CALL lbc_lnk(e3v_0 ,'T',1.) ; CALL lbc_lnk(e3f_0 ,'T',1.)
2321      CALL lbc_lnk(e3w_0 ,'T',1.)
2322      CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.)
2323      !
2324      !                                               ! =============
2325
2326      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2327      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2328
2329   END SUBROUTINE s_sf12
2330
2331   SUBROUTINE s_tanh()
2332
2333      !!----------------------------------------------------------------------
2334      !!                  ***  ROUTINE s_tanh***
2335      !!                     
2336      !! ** Purpose :   stretch the s-coordinate system
2337      !!
2338      !! ** Method  :   s-coordinate stretch
2339      !!
2340      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
2341      !!----------------------------------------------------------------------
2342
2343      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2344      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
2345
2346      REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
2347      REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw
2348
2349      CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      )
2350      CALL wrk_alloc( jpk, z_esigt, z_esigw                                               )
2351
2352      z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp   ;   z_gsi3w  = 0._wp
2353      z_esigt  = 0._wp   ;   z_esigw  = 0._wp 
2354
2355      DO jk = 1, jpk
2356        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
2357        z_gsigt(jk) = -fssig( REAL(jk,wp)        )
2358      END DO
2359      IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk)
2360      !
2361      ! Coefficients for vertical scale factors at w-, t- levels
2362!!gm bug :  define it from analytical function, not like juste bellow....
2363!!gm        or betteroffer the 2 possibilities....
2364      DO jk = 1, jpkm1
2365         z_esigt(jk  ) = z_gsigw(jk+1) - z_gsigw(jk)
2366         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
2367      END DO
2368      z_esigw( 1 ) = 2._wp * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
2369      z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
2370      !
2371      ! Coefficients for vertical depth as the sum of e3w scale factors
2372      z_gsi3w(1) = 0.5_wp * z_esigw(1)
2373      DO jk = 2, jpk
2374         z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
2375      END DO
2376!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
2377      DO jk = 1, jpk
2378         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2379         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2380         gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
2381         gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
2382         gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )
2383      END DO
2384!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
2385      DO jj = 1, jpj
2386         DO ji = 1, jpi
2387            DO jk = 1, jpk
2388              e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2389              e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2390              e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2391              e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
2392              !
2393              e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2394              e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2395              e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2396            END DO
2397         END DO
2398      END DO
2399
2400      CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      )
2401      CALL wrk_dealloc( jpk, z_esigt, z_esigw                                               )
2402
2403   END SUBROUTINE s_tanh
2404
2405   FUNCTION fssig( pk ) RESULT( pf )
2406      !!----------------------------------------------------------------------
2407      !!                 ***  ROUTINE fssig ***
2408      !!       
2409      !! ** Purpose :   provide the analytical function in s-coordinate
2410      !!         
2411      !! ** Method  :   the function provide the non-dimensional position of
2412      !!                T and W (i.e. between 0 and 1)
2413      !!                T-points at integer values (between 1 and jpk)
2414      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2415      !!----------------------------------------------------------------------
2416      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
2417      REAL(wp)             ::   pf   ! sigma value
2418      !!----------------------------------------------------------------------
2419      !
2420      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   &
2421         &     - TANH( rn_thetb * rn_theta                                )  )   &
2422         & * (   COSH( rn_theta                           )                      &
2423         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              &
2424         & / ( 2._wp * SINH( rn_theta ) )
2425      !
2426   END FUNCTION fssig
2427
2428
2429   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
2430      !!----------------------------------------------------------------------
2431      !!                 ***  ROUTINE fssig1 ***
2432      !!
2433      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
2434      !!
2435      !! ** Method  :   the function provides the non-dimensional position of
2436      !!                T and W (i.e. between 0 and 1)
2437      !!                T-points at integer values (between 1 and jpk)
2438      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2439      !!----------------------------------------------------------------------
2440      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
2441      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
2442      REAL(wp)             ::   pf1   ! sigma value
2443      !!----------------------------------------------------------------------
2444      !
2445      IF ( rn_theta == 0 ) then      ! uniform sigma
2446         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
2447      ELSE                        ! stretched sigma
2448         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
2449            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
2450            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
2451      ENDIF
2452      !
2453   END FUNCTION fssig1
2454
2455
2456   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
2457      !!----------------------------------------------------------------------
2458      !!                 ***  ROUTINE fgamma  ***
2459      !!
2460      !! ** Purpose :   provide analytical function for the s-coordinate
2461      !!
2462      !! ** Method  :   the function provides the non-dimensional position of
2463      !!                T and W (i.e. between 0 and 1)
2464      !!                T-points at integer values (between 1 and jpk)
2465      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2466      !!
2467      !!                This method allows the maintenance of fixed surface and or
2468      !!                bottom cell resolutions (cf. geopotential coordinates)
2469      !!                within an analytically derived stretched S-coordinate framework.
2470      !!
2471      !! Reference  :   Siddorn and Furner, in prep
2472      !!----------------------------------------------------------------------
2473      REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate
2474      REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate
2475      REAL(wp), INTENT(in   ) ::   pzb           ! Bottom box depth
2476      REAL(wp), INTENT(in   ) ::   pzs           ! surface box depth
2477      REAL(wp), INTENT(in   ) ::   psmth       ! Smoothing parameter
2478      REAL(wp)                ::   za1,za2,za3    ! local variables
2479      REAL(wp)                ::   zn1,zn2        ! local variables
2480      REAL(wp)                ::   za,zb,zx       ! local variables
2481      integer                 ::   jk
2482      !!----------------------------------------------------------------------
2483      !
2484
2485      zn1  =  1./(jpk-1.)
2486      zn2  =  1. -  zn1
2487
2488      za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 
2489      za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp)
2490      za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
2491     
2492      za = pzb - za3*(pzs-za1)-za2
2493      za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) )
2494      zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
2495      zx = 1.0_wp-za/2.0_wp-zb
2496 
2497      DO jk = 1, jpk
2498        p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp +  &
2499                    & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &
2500                    &      (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
2501        p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth)
2502      ENDDO 
2503
2504      !
2505   END FUNCTION fgamma
2506
2507   !!======================================================================
2508END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.