source: branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5111

Last change on this file since 5111 was 5111, checked in by mathiot, 6 years ago

add some missing if ln_isfcav, test of option compatibility with ln_isfcav + small documentation

  • Property svn:keywords set to Id
File size: 130.2 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   !                              !!* Namelist namzgr_sco *
53   LOGICAL  ::   ln_s_sh94         ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T)
54   LOGICAL  ::   ln_s_sf12         ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T)
55   !
56   REAL(wp) ::   rn_sbot_min       ! minimum depth of s-bottom surface (>0) (m)
57   REAL(wp) ::   rn_sbot_max       ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
58   REAL(wp) ::   rn_rmax           ! maximum cut-off r-value allowed (0<rn_rmax<1)
59   REAL(wp) ::   rn_hc             ! Critical depth for transition from sigma to stretched coordinates
60   ! Song and Haidvogel 1994 stretching parameters
61   REAL(wp) ::   rn_theta          ! surface control parameter (0<=rn_theta<=20)
62   REAL(wp) ::   rn_thetb          ! bottom control parameter  (0<=rn_thetb<= 1)
63   REAL(wp) ::   rn_bb             ! stretching parameter
64   !                                        ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
65   ! Siddorn and Furner stretching parameters
66   LOGICAL  ::   ln_sigcrit        ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
67   REAL(wp) ::   rn_alpha          ! control parameter ( > 1 stretch towards surface, < 1 towards seabed)
68   REAL(wp) ::   rn_efold          !  efold length scale for transition to stretched coord
69   REAL(wp) ::   rn_zs             !  depth of surface grid box
70                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
71   REAL(wp) ::   rn_zb_a           !  bathymetry scaling factor for calculating Zb
72   REAL(wp) ::   rn_zb_b           !  offset for calculating Zb
73
74  !! * Substitutions
75#  include "domzgr_substitute.h90"
76#  include "vectopt_loop_substitute.h90"
77   !!----------------------------------------------------------------------
78   !! NEMO/OPA 3.3.1 , NEMO Consortium (2011)
79   !! $Id$
80   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
81   !!----------------------------------------------------------------------
82CONTAINS       
83
84   SUBROUTINE dom_zgr
85      !!----------------------------------------------------------------------
86      !!                ***  ROUTINE dom_zgr  ***
87      !!                   
88      !! ** Purpose :   set the depth of model levels and the resulting
89      !!              vertical scale factors.
90      !!
91      !! ** Method  : - reference 1D vertical coordinate (gdep._1d, e3._1d)
92      !!              - read/set ocean depth and ocean levels (bathy, mbathy)
93      !!              - vertical coordinate (gdep., e3.) depending on the
94      !!                coordinate chosen :
95      !!                   ln_zco=T   z-coordinate   
96      !!                   ln_zps=T   z-coordinate with partial steps
97      !!                   ln_zco=T   s-coordinate
98      !!
99      !! ** Action  :   define gdep., e3., mbathy and bathy
100      !!----------------------------------------------------------------------
101      INTEGER ::   ioptio, ibat   ! local integer
102      INTEGER ::   ios
103      !
104      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav
105      !!----------------------------------------------------------------------
106      !
107      IF( nn_timing == 1 )   CALL timing_start('dom_zgr')
108      !
109      REWIND( numnam_ref )              ! Namelist namzgr in reference namelist : Vertical coordinate
110      READ  ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901 )
111901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist', lwp )
112
113      REWIND( numnam_cfg )              ! Namelist namzgr in configuration namelist : Vertical coordinate
114      READ  ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 )
115902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp )
116      IF(lwm) WRITE ( numond, namzgr )
117
118      IF(lwp) THEN                     ! Control print
119         WRITE(numout,*)
120         WRITE(numout,*) 'dom_zgr : vertical coordinate'
121         WRITE(numout,*) '~~~~~~~'
122         WRITE(numout,*) '          Namelist namzgr : set vertical coordinate'
123         WRITE(numout,*) '             z-coordinate - full steps      ln_zco    = ', ln_zco
124         WRITE(numout,*) '             z-coordinate - partial steps   ln_zps    = ', ln_zps
125         WRITE(numout,*) '             s- or hybrid z-s-coordinate    ln_sco    = ', ln_sco
126         WRITE(numout,*) '             ice shelf cavities             ln_isfcav = ', ln_isfcav
127      ENDIF
128
129      ioptio = 0                       ! Check Vertical coordinate options
130      IF( ln_zco      )   ioptio = ioptio + 1
131      IF( ln_zps      )   ioptio = ioptio + 1
132      IF( ln_sco      )   ioptio = ioptio + 1
133      IF( ioptio /= 1 )   CALL ctl_stop( ' none or several vertical coordinate options used' )
134      !
135      ! Build the vertical coordinate system
136      ! ------------------------------------
137                          CALL zgr_z            ! Reference z-coordinate system (always called)
138                          CALL zgr_bat          ! Bathymetry fields (levels and meters)
139      IF( lk_c1d      )   CALL lbc_lnk( bathy , 'T', 1._wp )   ! 1D config.: same bathy value over the 3x3 domain
140      IF( ln_zco      )   CALL zgr_zco          ! z-coordinate
141      IF( ln_zps      )   CALL zgr_zps          ! Partial step z-coordinate
142      IF( ln_sco      )   CALL zgr_sco          ! s-coordinate or hybrid z-s coordinate
143      !
144      ! final adjustment of mbathy & check
145      ! -----------------------------------
146      IF( lzoom       )   CALL zgr_bat_zoom     ! correct mbathy in case of zoom subdomain
147      IF( .NOT.lk_c1d )   CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isolated ocean points
148                          CALL zgr_bot_level    ! deepest ocean level for t-, u- and v-points
149                          CALL zgr_top_level    ! shallowest ocean level for T-, U-, V- points
150      !
151      IF( lk_c1d ) THEN                         ! 1D config.: same mbathy value over the 3x3 domain
152         ibat = mbathy(2,2)
153         mbathy(:,:) = ibat
154      END IF
155      !
156      IF( nprint == 1 .AND. lwp )   THEN
157         WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
158         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
159            &                   ' w ',   MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gdep3w_0(:,:,:) )
160         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL( e3f_0(:,:,:) ),  &
161            &                   ' u ',   MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL( e3v_0(:,:,:) ),  &
162            &                   ' uw',   MINVAL( e3uw_0(:,:,:)), ' vw', MINVAL( e3vw_0(:,:,:)),   &
163            &                   ' w ',   MINVAL( e3w_0(:,:,:) )
164
165         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
166            &                   ' w ',   MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gdep3w_0(:,:,:) )
167         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL( e3f_0(:,:,:) ),  &
168            &                   ' u ',   MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL( e3v_0(:,:,:) ),  &
169            &                   ' uw',   MAXVAL( e3uw_0(:,:,:)), ' vw', MAXVAL( e3vw_0(:,:,:)),   &
170            &                   ' w ',   MAXVAL( e3w_0(:,:,:) )
171      ENDIF
172      !
173      IF( nn_timing == 1 )  CALL timing_stop('dom_zgr')
174      !
175   END SUBROUTINE dom_zgr
176
177
178   SUBROUTINE zgr_z
179      !!----------------------------------------------------------------------
180      !!                   ***  ROUTINE zgr_z  ***
181      !!                   
182      !! ** Purpose :   set the depth of model levels and the resulting
183      !!      vertical scale factors.
184      !!
185      !! ** Method  :   z-coordinate system (use in all type of coordinate)
186      !!        The depth of model levels is defined from an analytical
187      !!      function the derivative of which gives the scale factors.
188      !!        both depth and scale factors only depend on k (1d arrays).
189      !!              w-level: gdepw_1d  = gdep(k)
190      !!                       e3w_1d(k) = dk(gdep)(k)     = e3(k)
191      !!              t-level: gdept_1d  = gdep(k+0.5)
192      !!                       e3t_1d(k) = dk(gdep)(k+0.5) = e3(k+0.5)
193      !!
194      !! ** Action  : - gdept_1d, gdepw_1d : depth of T- and W-point (m)
195      !!              - e3t_1d  , e3w_1d   : scale factors at T- and W-levels (m)
196      !!
197      !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
198      !!----------------------------------------------------------------------
199      INTEGER  ::   jk                     ! dummy loop indices
200      REAL(wp) ::   zt, zw                 ! temporary scalars
201      REAL(wp) ::   zsur, za0, za1, zkth   ! Values set from parameters in
202      REAL(wp) ::   zacr, zdzmin, zhmax    ! par_CONFIG_Rxx.h90
203      REAL(wp) ::   zrefdep                ! depth of the reference level (~10m)
204      REAL(wp) ::   za2, zkth2, zacr2      ! Values for optional double tanh function set from parameters
205      !!----------------------------------------------------------------------
206      !
207      IF( nn_timing == 1 )  CALL timing_start('zgr_z')
208      !
209      ! Set variables from parameters
210      ! ------------------------------
211       zkth = ppkth       ;   zacr = ppacr
212       zdzmin = ppdzmin   ;   zhmax = pphmax
213       zkth2 = ppkth2     ;   zacr2 = ppacr2   ! optional (ldbletanh=T) double tanh parameters
214
215      ! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
216      !  za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
217      IF(   ppa1  == pp_to_be_computed  .AND.  &
218         &  ppa0  == pp_to_be_computed  .AND.  &
219         &  ppsur == pp_to_be_computed           ) THEN
220         !
221         za1  = (  ppdzmin - pphmax / FLOAT(jpkm1)  )                                                      &
222            & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * (  LOG( COSH( (jpk - ppkth) / ppacr) )      &
223            &                                                   - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
224         za0  = ppdzmin - za1 *              TANH( (1-ppkth) / ppacr )
225         zsur =   - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
226      ELSE
227         za1 = ppa1 ;       za0 = ppa0 ;          zsur = ppsur
228         za2 = ppa2                            ! optional (ldbletanh=T) double tanh parameter
229      ENDIF
230
231      IF(lwp) THEN                         ! Parameter print
232         WRITE(numout,*)
233         WRITE(numout,*) '    zgr_z   : Reference vertical z-coordinates'
234         WRITE(numout,*) '    ~~~~~~~'
235         IF(  ppkth == 0._wp ) THEN             
236              WRITE(numout,*) '            Uniform grid with ',jpk-1,' layers'
237              WRITE(numout,*) '            Total depth    :', zhmax
238              WRITE(numout,*) '            Layer thickness:', zhmax/(jpk-1)
239         ELSE
240            IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN
241               WRITE(numout,*) '         zsur, za0, za1 computed from '
242               WRITE(numout,*) '                 zdzmin = ', zdzmin
243               WRITE(numout,*) '                 zhmax  = ', zhmax
244            ENDIF
245            WRITE(numout,*) '           Value of coefficients for vertical mesh:'
246            WRITE(numout,*) '                 zsur = ', zsur
247            WRITE(numout,*) '                 za0  = ', za0
248            WRITE(numout,*) '                 za1  = ', za1
249            WRITE(numout,*) '                 zkth = ', zkth
250            WRITE(numout,*) '                 zacr = ', zacr
251            IF( ldbletanh ) THEN
252               WRITE(numout,*) ' (Double tanh    za2  = ', za2
253               WRITE(numout,*) '  parameters)    zkth2= ', zkth2
254               WRITE(numout,*) '                 zacr2= ', zacr2
255            ENDIF
256         ENDIF
257      ENDIF
258
259
260      ! Reference z-coordinate (depth - scale factor at T- and W-points)
261      ! ======================
262      IF( ppkth == 0._wp ) THEN            !  uniform vertical grid       
263         za1 = zhmax / FLOAT(jpk-1) 
264         DO jk = 1, jpk
265            zw = FLOAT( jk )
266            zt = FLOAT( jk ) + 0.5_wp
267            gdepw_1d(jk) = ( zw - 1 ) * za1
268            gdept_1d(jk) = ( zt - 1 ) * za1
269            e3w_1d  (jk) =  za1
270            e3t_1d  (jk) =  za1
271         END DO
272      ELSE                                ! Madec & Imbard 1996 function
273         IF( .NOT. ldbletanh ) THEN
274            DO jk = 1, jpk
275               zw = REAL( jk , wp )
276               zt = REAL( jk , wp ) + 0.5_wp
277               gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) )  )
278               gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) )  )
279               e3w_1d  (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   )
280               e3t_1d  (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   )
281            END DO
282         ELSE
283            DO jk = 1, jpk
284               zw = FLOAT( jk )
285               zt = FLOAT( jk ) + 0.5_wp
286               ! Double tanh function
287               gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr  ) )    &
288                  &                             + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) )  )
289               gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr  ) )    &
290                  &                             + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) )  )
291               e3w_1d  (jk) =          za0      + za1        * TANH(       (zw-zkth ) / zacr  )      &
292                  &                             + za2        * TANH(       (zw-zkth2) / zacr2 )
293               e3t_1d  (jk) =          za0      + za1        * TANH(       (zt-zkth ) / zacr  )      &
294                  &                             + za2        * TANH(       (zt-zkth2) / zacr2 )
295            END DO
296         ENDIF
297         gdepw_1d(1) = 0._wp                    ! force first w-level to be exactly at zero
298      ENDIF
299
300!      IF ( ln_isfcav ) THEN
301! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth)
302! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively
303         DO jk = 1, jpkm1
304            e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) 
305         END DO
306         e3t_1d(jpk) = e3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO
307
308         DO jk = 2, jpk
309            e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1) 
310         END DO
311         e3w_1d(1  ) = 2._wp * (gdept_1d(1) - gdepw_1d(1)) 
312!      END IF
313
314!!gm BUG in s-coordinate this does not work!
315      ! deepest/shallowest W level Above/Below ~10m
316      zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d )                   ! ref. depth with tolerance (10% of minimum layer thickness)
317      nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
318      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m
319!!gm end bug
320
321      IF(lwp) THEN                        ! control print
322         WRITE(numout,*)
323         WRITE(numout,*) '              Reference z-coordinate depth and scale factors:'
324         WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
325         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk )
326      ENDIF
327      DO jk = 1, jpk                      ! control positivity
328         IF( e3w_1d  (jk) <= 0._wp .OR. e3t_1d  (jk) <= 0._wp )   CALL ctl_stop( 'dom:zgr_z: e3w_1d or e3t_1d =< 0 '    )
329         IF( gdepw_1d(jk) <  0._wp .OR. gdept_1d(jk) <  0._wp )   CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' )
330      END DO
331      !
332      IF( nn_timing == 1 )  CALL timing_stop('zgr_z')
333      !
334   END SUBROUTINE zgr_z
335
336
337   SUBROUTINE zgr_bat
338      !!----------------------------------------------------------------------
339      !!                    ***  ROUTINE zgr_bat  ***
340      !!
341      !! ** Purpose :   set bathymetry both in levels and meters
342      !!
343      !! ** Method  :   read or define mbathy and bathy arrays
344      !!       * level bathymetry:
345      !!      The ocean basin geometry is given by a two-dimensional array,
346      !!      mbathy, which is defined as follow :
347      !!            mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
348      !!                              at t-point (ji,jj).
349      !!                            = 0  over the continental t-point.
350      !!      The array mbathy is checked to verified its consistency with
351      !!      model option. in particular:
352      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
353      !!                  along closed boundary.
354      !!            mbathy must be cyclic IF jperio=1.
355      !!            mbathy must be lower or equal to jpk-1.
356      !!            isolated ocean grid points are suppressed from mbathy
357      !!                  since they are only connected to remaining
358      !!                  ocean through vertical diffusion.
359      !!      ntopo=-1 :   rectangular channel or bassin with a bump
360      !!      ntopo= 0 :   flat rectangular channel or basin
361      !!      ntopo= 1 :   mbathy is read in 'bathy_level.nc' NetCDF file
362      !!                   bathy  is read in 'bathy_meter.nc' NetCDF file
363      !!
364      !! ** Action  : - mbathy: level bathymetry (in level index)
365      !!              - bathy : meter bathymetry (in meters)
366      !!----------------------------------------------------------------------
367      INTEGER  ::   ji, jj, jl, jk            ! dummy loop indices
368      INTEGER  ::   inum                      ! temporary logical unit
369      INTEGER  ::   ierror                    ! error flag
370      INTEGER  ::   ii_bump, ij_bump, ih      ! bump center position
371      INTEGER  ::   ii0, ii1, ij0, ij1, ik    ! local indices
372      REAL(wp) ::   r_bump , h_bump , h_oce   ! bump characteristics
373      REAL(wp) ::   zi, zj, zh, zhmin         ! local scalars
374      INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   idta   ! global domain integer data
375      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdta   ! global domain scalar data
376      !!----------------------------------------------------------------------
377      !
378      IF( nn_timing == 1 )  CALL timing_start('zgr_bat')
379      !
380      IF(lwp) WRITE(numout,*)
381      IF(lwp) WRITE(numout,*) '    zgr_bat : defines level and meter bathymetry'
382      IF(lwp) WRITE(numout,*) '    ~~~~~~~'
383      !                                               ! ================== !
384      IF( ntopo == 0 .OR. ntopo == -1 ) THEN          !   defined by hand  !
385         !                                            ! ================== !
386         !                                            ! global domain level and meter bathymetry (idta,zdta)
387         !
388         ALLOCATE( idta(jpidta,jpjdta), STAT=ierror )
389         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate idta array' )
390         ALLOCATE( zdta(jpidta,jpjdta), STAT=ierror )
391         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zdta array' )
392         !
393         IF( ntopo == 0 ) THEN                        ! flat basin
394            IF(lwp) WRITE(numout,*)
395            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin'
396            IF( rn_bathy > 0.01 ) THEN
397               IF(lwp) WRITE(numout,*) '         Depth = rn_bathy read in namelist'
398               zdta(:,:) = rn_bathy
399               IF( ln_sco ) THEN                                   ! s-coordinate (zsc       ): idta()=jpk
400                  idta(:,:) = jpkm1
401               ELSE                                                ! z-coordinate (zco or zps): step-like topography
402                  idta(:,:) = jpkm1
403                  DO jk = 1, jpkm1
404                     WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) )   idta(:,:) = jk
405                  END DO
406               ENDIF
407            ELSE
408               IF(lwp) WRITE(numout,*) '         Depth = depthw(jpkm1)'
409               idta(:,:) = jpkm1                            ! before last level
410               zdta(:,:) = gdepw_1d(jpk)                     ! last w-point depth
411               h_oce     = gdepw_1d(jpk)
412            ENDIF
413         ELSE                                         ! bump centered in the basin
414            IF(lwp) WRITE(numout,*)
415            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin with a bump'
416            ii_bump = jpidta / 2                           ! i-index of the bump center
417            ij_bump = jpjdta / 2                           ! j-index of the bump center
418            r_bump  = 50000._wp                            ! bump radius (meters)       
419            h_bump  =  2700._wp                            ! bump height (meters)
420            h_oce   = gdepw_1d(jpk)                        ! background ocean depth (meters)
421            IF(lwp) WRITE(numout,*) '            bump characteristics: '
422            IF(lwp) WRITE(numout,*) '               bump center (i,j)   = ', ii_bump, ii_bump
423            IF(lwp) WRITE(numout,*) '               bump height         = ', h_bump , ' meters'
424            IF(lwp) WRITE(numout,*) '               bump radius         = ', r_bump , ' index'
425            IF(lwp) WRITE(numout,*) '            background ocean depth = ', h_oce  , ' meters'
426            !                                       
427            DO jj = 1, jpjdta                              ! zdta :
428               DO ji = 1, jpidta
429                  zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
430                  zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
431                  zdta(ji,jj) = h_oce - h_bump * EXP( -( zi*zi + zj*zj ) )
432               END DO
433            END DO
434            !                                              ! idta :
435            IF( ln_sco ) THEN                                   ! s-coordinate (zsc       ): idta()=jpk
436               idta(:,:) = jpkm1
437            ELSE                                                ! z-coordinate (zco or zps): step-like topography
438               idta(:,:) = jpkm1
439               DO jk = 1, jpkm1
440                  WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) )   idta(:,:) = jk
441               END DO
442            ENDIF
443         ENDIF
444         !                                            ! set GLOBAL boundary conditions
445         !                                            ! Caution : idta on the global domain: use of jperio, not nperio
446         IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
447            idta( :    , 1    ) = -1                ;      zdta( :    , 1    ) = -1._wp
448            idta( :    ,jpjdta) =  0                ;      zdta( :    ,jpjdta) =  0._wp
449         ELSEIF( jperio == 2 ) THEN
450            idta( :    , 1    ) = idta( : ,  3  )   ;      zdta( :    , 1    ) = zdta( : ,  3  )
451            idta( :    ,jpjdta) = 0                 ;      zdta( :    ,jpjdta) =  0._wp
452            idta( 1    , :    ) = 0                 ;      zdta( 1    , :    ) =  0._wp
453            idta(jpidta, :    ) = 0                 ;      zdta(jpidta, :    ) =  0._wp
454         ELSE
455            ih = 0                                  ;      zh = 0._wp
456            IF( ln_sco )   ih = jpkm1               ;      IF( ln_sco )   zh = h_oce
457            idta( :    , 1    ) = ih                ;      zdta( :    , 1    ) =  zh
458            idta( :    ,jpjdta) = ih                ;      zdta( :    ,jpjdta) =  zh
459            idta( 1    , :    ) = ih                ;      zdta( 1    , :    ) =  zh
460            idta(jpidta, :    ) = ih                ;      zdta(jpidta, :    ) =  zh
461         ENDIF
462
463         !                                            ! local domain level and meter bathymetries (mbathy,bathy)
464         mbathy(:,:) = 0                                   ! set to zero extra halo points
465         bathy (:,:) = 0._wp                               ! (require for mpp case)
466         DO jj = 1, nlcj                                   ! interior values
467            DO ji = 1, nlci
468               mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
469               bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
470            END DO
471         END DO
472         risfdep(:,:)=0.e0
473         misfdep(:,:)=1
474         !
475         ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code
476         IF( cp_cfg == "isomip" .AND. ln_isfcav ) THEN
477            risfdep(:,:)=200.e0 
478            misfdep(:,:)=1 
479            ij0 = 1 ; ij1 = 40 
480            DO jj = mj0(ij0), mj1(ij1) 
481               risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 
482            END DO
483            WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp 
484         !
485         ELSEIF ( cp_cfg == "isomip2" .AND. ln_isfcav ) THEN
486         !
487            risfdep(:,:)=0.e0
488            misfdep(:,:)=1
489            ij0 = 1 ; ij1 = 40
490            DO jj = mj0(ij0), mj1(ij1)
491               risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp
492            END DO
493            WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp
494         END IF
495         !
496         DEALLOCATE( idta, zdta )
497         !
498         !                                            ! ================ !
499      ELSEIF( ntopo == 1 ) THEN                       !   read in file   ! (over the local domain)
500         !                                            ! ================ !
501         !
502         IF( ln_zco )   THEN                          ! zco : read level bathymetry
503            CALL iom_open ( 'bathy_level.nc', inum ) 
504            CALL iom_get  ( inum, jpdom_data, 'Bathy_level', bathy )
505            CALL iom_close( inum )
506            mbathy(:,:) = INT( bathy(:,:) )
507            !                                                ! =====================
508            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
509               !                                             ! =====================
510               IF( nn_cla == 0 ) THEN
511                  ii0 = 140   ;   ii1 = 140                  ! Gibraltar Strait open
512                  ij0 = 102   ;   ij1 = 102                  ! (Thomson, Ocean Modelling, 1995)
513                  DO ji = mi0(ii0), mi1(ii1)
514                     DO jj = mj0(ij0), mj1(ij1)
515                        mbathy(ji,jj) = 15
516                     END DO
517                  END DO
518                  IF(lwp) WRITE(numout,*)
519                  IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
520                  !
521                  ii0 = 160   ;   ii1 = 160                  ! Bab el mandeb Strait open
522                  ij0 = 88    ;   ij1 = 88                   ! (Thomson, Ocean Modelling, 1995)
523                  DO ji = mi0(ii0), mi1(ii1)
524                     DO jj = mj0(ij0), mj1(ij1)
525                        mbathy(ji,jj) = 12
526                     END DO
527                  END DO
528                  IF(lwp) WRITE(numout,*)
529                  IF(lwp) WRITE(numout,*) '      orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
530               ENDIF
531               !
532            ENDIF
533            !
534         ENDIF
535         IF( ln_zps .OR. ln_sco )   THEN              ! zps or sco : read meter bathymetry
536            CALL iom_open ( 'bathy_meter.nc', inum ) 
537            CALL iom_get  ( inum, jpdom_data, 'Bathymetry', bathy )
538            CALL iom_close( inum )
539            !                                               
540            risfdep(:,:)=0._wp         
541            misfdep(:,:)=1             
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           ! 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 jk = 2,jpk                         
1151            DO jj = 1, jpjm1 
1152               DO ji = 1, fs_jpim1   ! vector opt.
1153                  e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj  ,jk) ) &
1154                    &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj  ,jk-1) )
1155                  e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji  ,jj+1,jk) ) &
1156                    &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji  ,jj+1,jk-1) )
1157               END DO
1158            END DO
1159         END DO
1160      !END IF
1161      !End (ISF)
1162     
1163      CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions
1164      CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp )
1165      !
1166      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
1167         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk)
1168         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk)
1169         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk)
1170         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk)
1171      END DO
1172     
1173      ! Scale factor at F-point
1174      DO jk = 1, jpk                        ! initialisation to z-scale factors
1175         e3f_0(:,:,jk) = e3t_1d(jk)
1176      END DO
1177      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors
1178         DO jj = 1, jpjm1
1179            DO ji = 1, fs_jpim1   ! vector opt.
1180               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
1181            END DO
1182         END DO
1183      END DO
1184      CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions
1185      !
1186      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
1187         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk)
1188      END DO
1189!!gm  bug ? :  must be a do loop with mj0,mj1
1190      !
1191      e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2
1192      e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 
1193      e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 
1194      e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 
1195      e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 
1196
1197      ! Control of the sign
1198      IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' )
1199      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' )
1200      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' )
1201      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' )
1202     
1203      ! Compute gdep3w_0 (vertical sum of e3w)
1204      IF ( ln_isfcav ) THEN ! if cavity
1205         WHERE (misfdep == 0) misfdep = 1
1206         DO jj = 1,jpj
1207            DO ji = 1,jpi
1208               gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)
1209               DO jk = 2, misfdep(ji,jj)
1210                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
1211               END DO
1212               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))
1213               DO jk = misfdep(ji,jj) + 1, jpk
1214                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
1215               END DO
1216            END DO
1217         END DO
1218      ELSE ! no cavity
1219         gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
1220         DO jk = 2, jpk
1221            gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk)
1222         END DO
1223      END IF
1224      !                                               ! ================= !
1225      IF(lwp .AND. ll_print) THEN                     !   Control print   !
1226         !                                            ! ================= !
1227         DO jj = 1,jpj
1228            DO ji = 1, jpi
1229               ik = MAX( mbathy(ji,jj), 1 )
1230               zprt(ji,jj,1) = e3t_0   (ji,jj,ik)
1231               zprt(ji,jj,2) = e3w_0   (ji,jj,ik)
1232               zprt(ji,jj,3) = e3u_0   (ji,jj,ik)
1233               zprt(ji,jj,4) = e3v_0   (ji,jj,ik)
1234               zprt(ji,jj,5) = e3f_0   (ji,jj,ik)
1235               zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)
1236            END DO
1237         END DO
1238         WRITE(numout,*)
1239         WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1240         WRITE(numout,*)
1241         WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1242         WRITE(numout,*)
1243         WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1244         WRITE(numout,*)
1245         WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1246         WRITE(numout,*)
1247         WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1248         WRITE(numout,*)
1249         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1250      ENDIF 
1251      !
1252      CALL wrk_dealloc( jpi, jpj, jpk, zprt )
1253      !
1254      IF( nn_timing == 1 )  CALL timing_stop('zgr_zps')
1255      !
1256   END SUBROUTINE zgr_zps
1257
1258   SUBROUTINE zgr_isf
1259      !!----------------------------------------------------------------------
1260      !!                    ***  ROUTINE zgr_bat_ctl  ***
1261      !!   
1262      !! ** Purpose :   check the bathymetry in levels
1263      !!   
1264      !! ** Method  :   THe water column have to contained at least 2 cells
1265      !!                Bathymetry and isfdraft are modified (dig/close) to respect
1266      !!                this criterion.
1267      !!                 
1268      !!   
1269      !! ** Action  : - test compatibility between isfdraft and bathy
1270      !!              - bathy and isfdraft are modified
1271      !!----------------------------------------------------------------------
1272      !!   
1273      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
1274      INTEGER  ::   ik, it           ! temporary integers
1275      INTEGER  ::   id, jd, nprocd
1276      INTEGER  ::   icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1   ! (ISF)
1277      LOGICAL  ::   ll_print         ! Allow  control print for debugging
1278      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
1279      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
1280      REAL(wp) ::   zmax, zmin       ! Maximum and minimum depth
1281      REAL(wp) ::   zdiff            ! temporary scalar
1282      REAL(wp) ::   zrefdep          ! temporary scalar
1283      REAL(wp) ::   zbathydiff, zrisfdepdiff  ! isf temporary scalar
1284      REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 2D workspace (ISH)
1285      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep         ! 2D workspace (ISH)
1286      !!---------------------------------------------------------------------
1287      !
1288      IF( nn_timing == 1 )  CALL timing_start('zgr_isf')
1289      !
1290      CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep)
1291      CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy )
1292
1293
1294      ! (ISF) compute misfdep
1295      WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1 
1296      ELSEWHERE                      ;                          misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level
1297      END WHERE 
1298
1299      ! Compute misfdep for ocean points (i.e. first wet level)
1300      ! find the first ocean level such that the first level thickness
1301      ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where
1302      ! e3t_0 is the reference level thickness
1303      DO jk = 2, jpkm1 
1304         zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
1305         WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1 
1306      END DO
1307      WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp)
1308         risfdep(:,:) = 0. ; misfdep(:,:) = 1
1309      END WHERE
1310 
1311! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation
1312      icompt = 0 
1313! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together
1314      DO jl = 1, 10     
1315         WHERE (bathy(:,:) .EQ. risfdep(:,:) )
1316            misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp
1317            mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp
1318         END WHERE
1319         WHERE (mbathy(:,:) <= 0) 
1320            misfdep(:,:) = 0; risfdep(:,:) = 0._wp 
1321            mbathy (:,:) = 0; bathy  (:,:) = 0._wp
1322         ENDWHERE
1323         IF( lk_mpp ) THEN
1324            zbathy(:,:) = FLOAT( misfdep(:,:) )
1325            CALL lbc_lnk( zbathy, 'T', 1. )
1326            misfdep(:,:) = INT( zbathy(:,:) )
1327            CALL lbc_lnk( risfdep, 'T', 1. )
1328            CALL lbc_lnk( bathy, 'T', 1. )
1329            zbathy(:,:) = FLOAT( mbathy(:,:) )
1330            CALL lbc_lnk( zbathy, 'T', 1. )
1331            mbathy(:,:) = INT( zbathy(:,:) )
1332         ENDIF
1333         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
1334            misfdep( 1 ,:) = misfdep(jpim1,:)           ! local domain is cyclic east-west
1335            misfdep(jpi,:) = misfdep(  2  ,:) 
1336         ENDIF
1337
1338         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
1339            mbathy( 1 ,:) = mbathy(jpim1,:)             ! local domain is cyclic east-west
1340            mbathy(jpi,:) = mbathy(  2  ,:)
1341         ENDIF
1342
1343         ! split last cell if possible (only where water column is 2 cell or less)
1344         DO jk = jpkm1, 1, -1
1345            zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1346            WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)
1347               mbathy(:,:) = jk
1348               bathy(:,:)  = zmax
1349            END WHERE
1350         END DO
1351 
1352         ! split top cell if possible (only where water column is 2 cell or less)
1353         DO jk = 2, jpkm1
1354            zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1355            WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy)
1356               misfdep(:,:) = jk
1357               risfdep(:,:) = zmax
1358            END WHERE
1359         END DO
1360
1361 
1362 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition
1363         DO jj = 1, jpj
1364            DO ji = 1, jpi
1365               ! find the minimum change option:
1366               ! test bathy
1367               IF (risfdep(ji,jj) .GT. 1) THEN
1368               zbathydiff =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) &
1369                 &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
1370               zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) &
1371                 &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
1372 
1373                  IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN
1374                     IF (zbathydiff .LE. zrisfdepdiff) THEN
1375                        bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )
1376                        mbathy(ji,jj)= mbathy(ji,jj) + 1
1377                     ELSE
1378                        risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )
1379                        misfdep(ji,jj) = misfdep(ji,jj) - 1
1380                     END IF
1381                  END IF
1382               END IF
1383            END DO
1384         END DO
1385 
1386          ! At least 2 levels for water thickness at T, U, and V point.
1387         DO jj = 1, jpj
1388            DO ji = 1, jpi
1389               ! find the minimum change option:
1390               ! test bathy
1391               IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
1392                  zbathydiff  =ABS(bathy(ji,jj)  - (gdepw_1d(mbathy (ji,jj)+1)&
1393                    & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
1394                  zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  )  &
1395                    & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
1396                  IF (zbathydiff .LE. zrisfdepdiff) THEN
1397                     mbathy(ji,jj) = mbathy(ji,jj) + 1
1398                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
1399                  ELSE
1400                     misfdep(ji,jj)= misfdep(ji,jj) - 1
1401                     risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )
1402                  END IF
1403               ENDIF
1404            END DO
1405         END DO
1406 
1407 ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1)
1408         DO jj = 1, jpjm1
1409            DO ji = 1, jpim1
1410               IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
1411                  zbathydiff =ABS(bathy(ji,jj    ) - (gdepw_1d(mbathy (ji,jj)+1)   &
1412                    &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat )))
1413                  zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1))  &
1414                    &  - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))
1415                  IF (zbathydiff .LE. zrisfdepdiff) THEN
1416                     mbathy(ji,jj) = mbathy(ji,jj) + 1
1417                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) &
1418                   &    + MIN( e3zps_min, e3t_1d(mbathy(ji,jj   )+1)*e3zps_rat )
1419                  ELSE
1420                     misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1
1421                     risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) &
1422                   &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )
1423                  END IF
1424               ENDIF
1425            END DO
1426         END DO
1427 
1428         IF( lk_mpp ) THEN
1429            zbathy(:,:) = FLOAT( misfdep(:,:) )
1430            CALL lbc_lnk( zbathy, 'T', 1. )
1431            misfdep(:,:) = INT( zbathy(:,:) )
1432            CALL lbc_lnk( risfdep, 'T', 1. )
1433            CALL lbc_lnk( bathy, 'T', 1. )
1434            zbathy(:,:) = FLOAT( mbathy(:,:) )
1435            CALL lbc_lnk( zbathy, 'T', 1. )
1436            mbathy(:,:) = INT( zbathy(:,:) )
1437         ENDIF
1438 ! point V misdep(ji,jj) EQ mbathy(ji,jj+1)
1439         DO jj = 1, jpjm1
1440            DO ji = 1, jpim1
1441               IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN
1442                  zbathydiff =ABS(  bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) &
1443                   &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))
1444                  zrisfdepdiff=ABS(risfdep(ji,jj  ) - (gdepw_1d(misfdep(ji,jj  )  )  &
1445                   &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat )))
1446                  IF (zbathydiff .LE. zrisfdepdiff) THEN
1447                     mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1
1448                     bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )
1449                  ELSE
1450                     misfdep(ji,jj)   = misfdep(ji,jj) - 1
1451                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat )
1452                  END IF
1453               ENDIF
1454            END DO
1455         END DO
1456 
1457 
1458         IF( lk_mpp ) THEN
1459            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
1460            CALL lbc_lnk( zbathy, 'T', 1. ) 
1461            misfdep(:,:) = INT( zbathy(:,:) ) 
1462            CALL lbc_lnk( risfdep, 'T', 1. ) 
1463            CALL lbc_lnk( bathy, 'T', 1. )
1464            zbathy(:,:) = FLOAT( mbathy(:,:) )
1465            CALL lbc_lnk( zbathy, 'T', 1. )
1466            mbathy(:,:) = INT( zbathy(:,:) )
1467         ENDIF 
1468 
1469 ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1)
1470         DO jj = 1, jpjm1
1471            DO ji = 1, jpim1
1472               IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
1473                  zbathydiff =ABS(  bathy(ji  ,jj) - (gdepw_1d(mbathy (ji,jj)+1) &
1474                    &   + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat )))
1475                  zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) &
1476                    &  - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))
1477                  IF (zbathydiff .LE. zrisfdepdiff) THEN
1478                     mbathy(ji,jj) = mbathy(ji,jj) + 1
1479                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
1480                  ELSE
1481                     misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1
1482                     risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )
1483                  END IF
1484               ENDIF
1485            ENDDO
1486         ENDDO
1487 
1488         IF( lk_mpp ) THEN
1489            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
1490            CALL lbc_lnk( zbathy, 'T', 1. ) 
1491            misfdep(:,:) = INT( zbathy(:,:) ) 
1492            CALL lbc_lnk( risfdep, 'T', 1. ) 
1493            CALL lbc_lnk( bathy, 'T', 1. )
1494            zbathy(:,:) = FLOAT( mbathy(:,:) )
1495            CALL lbc_lnk( zbathy, 'T', 1. )
1496            mbathy(:,:) = INT( zbathy(:,:) )
1497         ENDIF 
1498 
1499 ! point U misfdep(ji,jj) EQ bathy(ji,jj+1)
1500         DO jj = 1, jpjm1
1501            DO ji = 1, jpim1
1502               IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
1503                  zbathydiff =ABS(  bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) &
1504                      &   + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))
1505                  zrisfdepdiff=ABS(risfdep(ji  ,jj) - (gdepw_1d(misfdep(ji  ,jj)  )  &
1506                      &  - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat )))
1507                  IF (zbathydiff .LE. zrisfdepdiff) THEN
1508                     mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1
1509                     bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  )  &
1510                      &   + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat )
1511                  ELSE
1512                     misfdep(ji,jj)   = misfdep(ji  ,jj) - 1
1513                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) &
1514                      &   - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat )
1515                  END IF
1516               ENDIF
1517            ENDDO
1518         ENDDO
1519 
1520         IF( lk_mpp ) THEN
1521            zbathy(:,:) = FLOAT( misfdep(:,:) )
1522            CALL lbc_lnk( zbathy, 'T', 1. )
1523            misfdep(:,:) = INT( zbathy(:,:) )
1524            CALL lbc_lnk( risfdep, 'T', 1. )
1525            CALL lbc_lnk( bathy, 'T', 1. )
1526            zbathy(:,:) = FLOAT( mbathy(:,:) )
1527            CALL lbc_lnk( zbathy, 'T', 1. )
1528            mbathy(:,:) = INT( zbathy(:,:) )
1529         ENDIF
1530      END DO
1531      ! end dig bathy/ice shelf to be compatible
1532      ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness
1533      DO jl = 1,20
1534 
1535 ! remove single point "bay" on isf coast line in the ice shelf draft'
1536         DO jk = 1, jpk
1537            WHERE (misfdep==0) misfdep=jpk
1538            zmask=0
1539            WHERE (misfdep .LE. jk) zmask=1
1540            DO jj = 2, jpjm1
1541               DO ji = 2, jpim1
1542                  IF (misfdep(ji,jj) .EQ. jk) THEN
1543                     ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
1544                     IF (ibtest .LE. 1) THEN
1545                        risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1
1546                        IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) misfdep(ji,jj) = jpk
1547                     END IF
1548                  END IF
1549               END DO
1550            END DO
1551         END DO
1552         WHERE (misfdep==jpk)
1553             misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.
1554         END WHERE
1555         IF( lk_mpp ) THEN
1556            zbathy(:,:) = FLOAT( misfdep(:,:) )
1557            CALL lbc_lnk( zbathy, 'T', 1. )
1558            misfdep(:,:) = INT( zbathy(:,:) )
1559            CALL lbc_lnk( risfdep, 'T', 1. )
1560            CALL lbc_lnk( bathy, 'T', 1. )
1561            zbathy(:,:) = FLOAT( mbathy(:,:) )
1562            CALL lbc_lnk( zbathy, 'T', 1. )
1563            mbathy(:,:) = INT( zbathy(:,:) )
1564         ENDIF
1565 
1566 ! remove single point "bay" on bathy coast line beneath an ice shelf'
1567         DO jk = jpk,1,-1
1568            zmask=0
1569            WHERE (mbathy .GE. jk ) zmask=1
1570            DO jj = 2, jpjm1
1571               DO ji = 2, jpim1
1572                  IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) THEN
1573                     ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
1574                     IF (ibtest .LE. 1) THEN
1575                        bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1
1576                        IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 0
1577                     END IF
1578                  END IF
1579               END DO
1580            END DO
1581         END DO
1582         WHERE (mbathy==0)
1583             misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.
1584         END WHERE
1585         IF( lk_mpp ) THEN
1586            zbathy(:,:) = FLOAT( misfdep(:,:) )
1587            CALL lbc_lnk( zbathy, 'T', 1. )
1588            misfdep(:,:) = INT( zbathy(:,:) )
1589            CALL lbc_lnk( risfdep, 'T', 1. )
1590            CALL lbc_lnk( bathy, 'T', 1. )
1591            zbathy(:,:) = FLOAT( mbathy(:,:) )
1592            CALL lbc_lnk( zbathy, 'T', 1. )
1593            mbathy(:,:) = INT( zbathy(:,:) )
1594         ENDIF
1595 
1596 ! fill hole in ice shelf
1597         zmisfdep = misfdep
1598         zrisfdep = risfdep
1599         WHERE (zmisfdep .LE. 1) zmisfdep=jpk
1600         DO jj = 2, jpjm1
1601            DO ji = 2, jpim1
1602               ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  )
1603               ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1)
1604               IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj  ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj  ) - 1)
1605               IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj  ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj  ) - 1)
1606               IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji  ,jj-1) - 1)
1607               IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji  ,jj+1) - 1)
1608               ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
1609               IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN
1610                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp
1611               END IF
1612               IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN
1613                  misfdep(ji,jj) = ibtest
1614                  risfdep(ji,jj) = gdepw_1d(ibtest)
1615               ENDIF
1616            ENDDO
1617         ENDDO
1618 
1619         IF( lk_mpp ) THEN
1620            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
1621            CALL lbc_lnk( zbathy, 'T', 1. ) 
1622            misfdep(:,:) = INT( zbathy(:,:) ) 
1623            CALL lbc_lnk( risfdep, 'T', 1. ) 
1624            CALL lbc_lnk( bathy, 'T', 1. )
1625            zbathy(:,:) = FLOAT( mbathy(:,:) )
1626            CALL lbc_lnk( zbathy, 'T', 1. )
1627            mbathy(:,:) = INT( zbathy(:,:) )
1628         ENDIF 
1629 !
1630 !! fill hole in bathymetry
1631         zmbathy (:,:)=mbathy (:,:)
1632         DO jj = 2, jpjm1
1633            DO ji = 2, jpim1
1634               ibtestim1 = zmbathy(ji-1,jj  ) ; ibtestip1 = zmbathy(ji+1,jj  )
1635               ibtestjm1 = zmbathy(ji  ,jj-1) ; ibtestjp1 = zmbathy(ji  ,jj+1)
1636               IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj  ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj  ) + 1)
1637               IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj  ) ) ibtestip1 = 0
1638               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj-1) ) ibtestjm1 = 0
1639               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj+1) ) ibtestjp1 = 0
1640               ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
1641               IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN
1642                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ;
1643               END IF
1644               IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN
1645                  mbathy(ji,jj) = ibtest
1646                  bathy(ji,jj)  = gdepw_1d(ibtest+1) 
1647               ENDIF
1648            END DO
1649         END DO
1650         IF( lk_mpp ) THEN
1651            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
1652            CALL lbc_lnk( zbathy, 'T', 1. ) 
1653            misfdep(:,:) = INT( zbathy(:,:) ) 
1654            CALL lbc_lnk( risfdep, 'T', 1. ) 
1655            CALL lbc_lnk( bathy, 'T', 1. )
1656            zbathy(:,:) = FLOAT( mbathy(:,:) )
1657            CALL lbc_lnk( zbathy, 'T', 1. )
1658            mbathy(:,:) = INT( zbathy(:,:) )
1659         ENDIF 
1660 ! if not compatible after all check (ie U point water column less than 2 cells), mask U
1661         DO jj = 1, jpjm1
1662            DO ji = 1, jpim1
1663               IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN
1664                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ;
1665               END IF
1666            END DO
1667         END DO
1668         IF( lk_mpp ) THEN
1669            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
1670            CALL lbc_lnk( zbathy, 'T', 1. ) 
1671            misfdep(:,:) = INT( zbathy(:,:) ) 
1672            CALL lbc_lnk( risfdep, 'T', 1. ) 
1673            CALL lbc_lnk( bathy, 'T', 1. )
1674            zbathy(:,:) = FLOAT( mbathy(:,:) )
1675            CALL lbc_lnk( zbathy, 'T', 1. )
1676            mbathy(:,:) = INT( zbathy(:,:) )
1677         ENDIF 
1678 ! if not compatible after all check (ie U point water column less than 2 cells), mask U
1679         DO jj = 1, jpjm1
1680            DO ji = 1, jpim1
1681               IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN
1682                  mbathy(ji+1,jj)  = mbathy(ji+1,jj) - 1;   bathy(ji+1,jj)   = gdepw_1d(mbathy(ji+1,jj)+1) ;
1683               END IF
1684            END DO
1685         END DO
1686         IF( lk_mpp ) THEN
1687            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
1688            CALL lbc_lnk( zbathy, 'T', 1. ) 
1689            misfdep(:,:) = INT( zbathy(:,:) ) 
1690            CALL lbc_lnk( risfdep, 'T', 1. ) 
1691            CALL lbc_lnk( bathy, 'T', 1. )
1692            zbathy(:,:) = FLOAT( mbathy(:,:) )
1693            CALL lbc_lnk( zbathy, 'T', 1. )
1694            mbathy(:,:) = INT( zbathy(:,:) )
1695         ENDIF 
1696 ! if not compatible after all check (ie V point water column less than 2 cells), mask V
1697         DO jj = 1, jpjm1
1698            DO ji = 1, jpi
1699               IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN
1700                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ;
1701               END IF
1702            END DO
1703         END DO
1704         IF( lk_mpp ) THEN
1705            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
1706            CALL lbc_lnk( zbathy, 'T', 1. ) 
1707            misfdep(:,:) = INT( zbathy(:,:) ) 
1708            CALL lbc_lnk( risfdep, 'T', 1. ) 
1709            CALL lbc_lnk( bathy, 'T', 1. )
1710            zbathy(:,:) = FLOAT( mbathy(:,:) )
1711            CALL lbc_lnk( zbathy, 'T', 1. )
1712            mbathy(:,:) = INT( zbathy(:,:) )
1713         ENDIF 
1714 ! if not compatible after all check (ie V point water column less than 2 cells), mask V
1715         DO jj = 1, jpjm1
1716            DO ji = 1, jpi
1717               IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN
1718                  mbathy(ji,jj+1)  = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1)   = gdepw_1d(mbathy(ji,jj+1)+1) ;
1719               END IF
1720            END DO
1721         END DO
1722         IF( lk_mpp ) THEN
1723            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
1724            CALL lbc_lnk( zbathy, 'T', 1. ) 
1725            misfdep(:,:) = INT( zbathy(:,:) ) 
1726            CALL lbc_lnk( risfdep, 'T', 1. ) 
1727            CALL lbc_lnk( bathy, 'T', 1. )
1728            zbathy(:,:) = FLOAT( mbathy(:,:) )
1729            CALL lbc_lnk( zbathy, 'T', 1. )
1730            mbathy(:,:) = INT( zbathy(:,:) )
1731         ENDIF 
1732 ! if not compatible after all check, mask T
1733         DO jj = 1, jpj
1734            DO ji = 1, jpi
1735               IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN
1736                  misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj)  = 0 ; bathy(ji,jj)   = 0._wp ;
1737               END IF
1738            END DO
1739         END DO
1740 
1741         WHERE (mbathy(:,:) == 1)
1742            mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp
1743         END WHERE
1744      END DO 
1745! end check compatibility ice shelf/bathy
1746      ! remove very shallow ice shelf (less than ~ 10m if 75L)
1747      WHERE (misfdep(:,:) <= 5)
1748         misfdep = 1; risfdep = 0.0_wp;
1749      END WHERE
1750
1751      IF( icompt == 0 ) THEN
1752         IF(lwp) WRITE(numout,*)'     no points with ice shelf too close to bathymetry' 
1753      ELSE
1754         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry' 
1755      ENDIF
1756
1757      CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep )
1758      CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy )
1759
1760      IF( nn_timing == 1 )  CALL timing_stop('zgr_isf')
1761     
1762   END SUBROUTINE
1763
1764   SUBROUTINE zgr_sco
1765      !!----------------------------------------------------------------------
1766      !!                  ***  ROUTINE zgr_sco  ***
1767      !!                     
1768      !! ** Purpose :   define the s-coordinate system
1769      !!
1770      !! ** Method  :   s-coordinate
1771      !!         The depth of model levels is defined as the product of an
1772      !!      analytical function by the local bathymetry, while the vertical
1773      !!      scale factors are defined as the product of the first derivative
1774      !!      of the analytical function by the bathymetry.
1775      !!      (this solution save memory as depth and scale factors are not
1776      !!      3d fields)
1777      !!          - Read bathymetry (in meters) at t-point and compute the
1778      !!         bathymetry at u-, v-, and f-points.
1779      !!            hbatu = mi( hbatt )
1780      !!            hbatv = mj( hbatt )
1781      !!            hbatf = mi( mj( hbatt ) )
1782      !!          - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
1783      !!         function and its derivative given as function.
1784      !!            z_gsigt(k) = fssig (k    )
1785      !!            z_gsigw(k) = fssig (k-0.5)
1786      !!            z_esigt(k) = fsdsig(k    )
1787      !!            z_esigw(k) = fsdsig(k-0.5)
1788      !!      Three options for stretching are give, and they can be modified
1789      !!      following the users requirements. Nevertheless, the output as
1790      !!      well as the way to compute the model levels and scale factors
1791      !!      must be respected in order to insure second order accuracy
1792      !!      schemes.
1793      !!
1794      !!      The three methods for stretching available are:
1795      !!
1796      !!           s_sh94 (Song and Haidvogel 1994)
1797      !!                a sinh/tanh function that allows sigma and stretched sigma
1798      !!
1799      !!           s_sf12 (Siddorn and Furner 2012?)
1800      !!                allows the maintenance of fixed surface and or
1801      !!                bottom cell resolutions (cf. geopotential coordinates)
1802      !!                within an analytically derived stretched S-coordinate framework.
1803      !!
1804      !!          s_tanh  (Madec et al 1996)
1805      !!                a cosh/tanh function that gives stretched coordinates       
1806      !!
1807      !!----------------------------------------------------------------------
1808      !
1809      INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
1810      INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
1811      INTEGER  ::   ios                      ! Local integer output status for namelist read
1812      REAL(wp) ::   zrmax, ztaper   ! temporary scalars
1813      REAL(wp) ::   zrfact
1814      !
1815      REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
1816      REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
1817
1818      NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
1819                           rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
1820     !!----------------------------------------------------------------------
1821      !
1822      IF( nn_timing == 1 )  CALL timing_start('zgr_sco')
1823      !
1824      CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
1825      !
1826      REWIND( numnam_ref )              ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
1827      READ  ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901)
1828901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist', lwp )
1829
1830      REWIND( numnam_cfg )              ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters
1831      READ  ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 )
1832902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist', lwp )
1833      IF(lwm) WRITE ( numond, namzgr_sco )
1834
1835      IF(lwp) THEN                           ! control print
1836         WRITE(numout,*)
1837         WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
1838         WRITE(numout,*) '~~~~~~~~~~~'
1839         WRITE(numout,*) '   Namelist namzgr_sco'
1840         WRITE(numout,*) '     stretching coeffs '
1841         WRITE(numout,*) '        maximum depth of s-bottom surface (>0)       rn_sbot_max   = ',rn_sbot_max
1842         WRITE(numout,*) '        minimum depth of s-bottom surface (>0)       rn_sbot_min   = ',rn_sbot_min
1843         WRITE(numout,*) '        Critical depth                               rn_hc         = ',rn_hc
1844         WRITE(numout,*) '        maximum cut-off r-value allowed              rn_rmax       = ',rn_rmax
1845         WRITE(numout,*) '     Song and Haidvogel 1994 stretching              ln_s_sh94     = ',ln_s_sh94
1846         WRITE(numout,*) '        Song and Haidvogel 1994 stretching coefficients'
1847         WRITE(numout,*) '        surface control parameter (0<=rn_theta<=20)  rn_theta      = ',rn_theta
1848         WRITE(numout,*) '        bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ',rn_thetb
1849         WRITE(numout,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb
1850         WRITE(numout,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12
1851         WRITE(numout,*) '        switching to sigma (T) or Z (F) at H<Hc      ln_sigcrit    = ',ln_sigcrit
1852         WRITE(numout,*) '        Siddorn and Furner 2012 stretching coefficients'
1853         WRITE(numout,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha
1854         WRITE(numout,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold
1855         WRITE(numout,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs
1856         WRITE(numout,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a
1857         WRITE(numout,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b
1858         WRITE(numout,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b'
1859      ENDIF
1860
1861      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate
1862      hifu(:,:) = rn_sbot_min
1863      hifv(:,:) = rn_sbot_min
1864      hiff(:,:) = rn_sbot_min
1865
1866      !                                        ! set maximum ocean depth
1867      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1868
1869      DO jj = 1, jpj
1870         DO ji = 1, jpi
1871           IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1872         END DO
1873      END DO
1874      !                                        ! =============================
1875      !                                        ! Define the envelop bathymetry   (hbatt)
1876      !                                        ! =============================
1877      ! use r-value to create hybrid coordinates
1878      zenv(:,:) = bathy(:,:)
1879      !
1880      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
1881      DO jj = 1, jpj
1882         DO ji = 1, jpi
1883           IF( bathy(ji,jj) == 0._wp ) THEN
1884             iip1 = MIN( ji+1, jpi )
1885             ijp1 = MIN( jj+1, jpj )
1886             iim1 = MAX( ji-1, 1 )
1887             ijm1 = MAX( jj-1, 1 )
1888             IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) +              &
1889        &         bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
1890               zenv(ji,jj) = rn_sbot_min
1891             ENDIF
1892           ENDIF
1893         END DO
1894      END DO
1895      ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1896      CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1897      !
1898      ! smooth the bathymetry (if required)
1899      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
1900      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
1901      !
1902      jl = 0
1903      zrmax = 1._wp
1904      !   
1905      !     
1906      ! set scaling factor used in reducing vertical gradients
1907      zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
1908      !
1909      ! initialise temporary evelope depth arrays
1910      ztmpi1(:,:) = zenv(:,:)
1911      ztmpi2(:,:) = zenv(:,:)
1912      ztmpj1(:,:) = zenv(:,:)
1913      ztmpj2(:,:) = zenv(:,:)
1914      !
1915      ! initialise temporary r-value arrays
1916      zri(:,:) = 1._wp
1917      zrj(:,:) = 1._wp
1918      !                                                            ! ================ !
1919      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) !  Iterative loop  !
1920         !                                                         ! ================ !
1921         jl = jl + 1
1922         zrmax = 0._wp
1923         ! we set zrmax from previous r-values (zri and zrj) first
1924         ! if set after current r-value calculation (as previously)
1925         ! we could exit DO WHILE prematurely before checking r-value
1926         ! of current zenv
1927         DO jj = 1, nlcj
1928            DO ji = 1, nlci
1929               zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
1930            END DO
1931         END DO
1932         zri(:,:) = 0._wp
1933         zrj(:,:) = 0._wp
1934         DO jj = 1, nlcj
1935            DO ji = 1, nlci
1936               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
1937               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
1938               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
1939                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
1940               END IF
1941               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
1942                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
1943               END IF
1944               IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact
1945               IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact
1946               IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact
1947               IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact
1948            END DO
1949         END DO
1950         IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
1951         !
1952         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax
1953         !
1954         DO jj = 1, nlcj
1955            DO ji = 1, nlci
1956               zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
1957            END DO
1958         END DO
1959         ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1960         CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1961         !                                                  ! ================ !
1962      END DO                                                !     End loop     !
1963      !                                                     ! ================ !
1964      DO jj = 1, jpj
1965         DO ji = 1, jpi
1966            zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
1967         END DO
1968      END DO
1969      !
1970      ! Envelope bathymetry saved in hbatt
1971      hbatt(:,:) = zenv(:,:) 
1972      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
1973         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1974         DO jj = 1, jpj
1975            DO ji = 1, jpi
1976               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
1977               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
1978            END DO
1979         END DO
1980      ENDIF
1981      !
1982      IF(lwp) THEN                             ! Control print
1983         WRITE(numout,*)
1984         WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
1985         WRITE(numout,*)
1986         CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )
1987         IF( nprint == 1 )   THEN       
1988            WRITE(numout,*) ' bathy  MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
1989            WRITE(numout,*) ' hbatt  MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
1990         ENDIF
1991      ENDIF
1992
1993      !                                        ! ==============================
1994      !                                        !   hbatu, hbatv, hbatf fields
1995      !                                        ! ==============================
1996      IF(lwp) THEN
1997         WRITE(numout,*)
1998         WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
1999      ENDIF
2000      hbatu(:,:) = rn_sbot_min
2001      hbatv(:,:) = rn_sbot_min
2002      hbatf(:,:) = rn_sbot_min
2003      DO jj = 1, jpjm1
2004        DO ji = 1, jpim1   ! NO vector opt.
2005           hbatu(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
2006           hbatv(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
2007           hbatf(ji,jj) = 0.25_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
2008              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
2009        END DO
2010      END DO
2011      !
2012      ! Apply lateral boundary condition
2013!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
2014      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk( hbatu, 'U', 1._wp )
2015      DO jj = 1, jpj
2016         DO ji = 1, jpi
2017            IF( hbatu(ji,jj) == 0._wp ) THEN
2018               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
2019               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
2020            ENDIF
2021         END DO
2022      END DO
2023      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk( hbatv, 'V', 1._wp )
2024      DO jj = 1, jpj
2025         DO ji = 1, jpi
2026            IF( hbatv(ji,jj) == 0._wp ) THEN
2027               IF( zhbat(ji,jj) == 0._wp )   hbatv(ji,jj) = rn_sbot_min
2028               IF( zhbat(ji,jj) /= 0._wp )   hbatv(ji,jj) = zhbat(ji,jj)
2029            ENDIF
2030         END DO
2031      END DO
2032      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk( hbatf, 'F', 1._wp )
2033      DO jj = 1, jpj
2034         DO ji = 1, jpi
2035            IF( hbatf(ji,jj) == 0._wp ) THEN
2036               IF( zhbat(ji,jj) == 0._wp )   hbatf(ji,jj) = rn_sbot_min
2037               IF( zhbat(ji,jj) /= 0._wp )   hbatf(ji,jj) = zhbat(ji,jj)
2038            ENDIF
2039         END DO
2040      END DO
2041
2042!!bug:  key_helsinki a verifer
2043      hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
2044      hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
2045      hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
2046      hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
2047
2048      IF( nprint == 1 .AND. lwp )   THEN
2049         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
2050            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
2051         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
2052            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
2053         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
2054            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
2055         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
2056            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
2057      ENDIF
2058!! helsinki
2059
2060      !                                            ! =======================
2061      !                                            !   s-ccordinate fields     (gdep., e3.)
2062      !                                            ! =======================
2063      !
2064      ! non-dimensional "sigma" for model level depth at w- and t-levels
2065
2066
2067!========================================================================
2068! Song and Haidvogel  1994 (ln_s_sh94=T)
2069! Siddorn and Furner 2012 (ln_sf12=T)
2070! or  tanh function       (both false)                   
2071!========================================================================
2072      IF      ( ln_s_sh94 ) THEN
2073                           CALL s_sh94()
2074      ELSE IF ( ln_s_sf12 ) THEN
2075                           CALL s_sf12()
2076      ELSE                 
2077                           CALL s_tanh()
2078      ENDIF
2079
2080      CALL lbc_lnk( e3t_0 , 'T', 1._wp )
2081      CALL lbc_lnk( e3u_0 , 'U', 1._wp )
2082      CALL lbc_lnk( e3v_0 , 'V', 1._wp )
2083      CALL lbc_lnk( e3f_0 , 'F', 1._wp )
2084      CALL lbc_lnk( e3w_0 , 'W', 1._wp )
2085      CALL lbc_lnk( e3uw_0, 'U', 1._wp )
2086      CALL lbc_lnk( e3vw_0, 'V', 1._wp )
2087
2088      fsdepw(:,:,:) = gdepw_0 (:,:,:)
2089      fsde3w(:,:,:) = gdep3w_0(:,:,:)
2090      !
2091      where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0
2092      where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0
2093      where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0
2094      where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0
2095      where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0
2096      where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0
2097      where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0
2098
2099#if defined key_agrif
2100      ! Ensure meaningful vertical scale factors in ghost lines/columns
2101      IF( .NOT. Agrif_Root() ) THEN
2102         
2103         IF((nbondi == -1).OR.(nbondi == 2)) THEN
2104            e3u_0(1,:,:) = e3u_0(2,:,:)
2105         ENDIF
2106         !
2107         IF((nbondi ==  1).OR.(nbondi == 2)) THEN
2108            e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:)
2109         ENDIF
2110         !
2111         IF((nbondj == -1).OR.(nbondj == 2)) THEN
2112            e3v_0(:,1,:) = e3v_0(:,2,:)
2113         ENDIF
2114         !
2115         IF((nbondj ==  1).OR.(nbondj == 2)) THEN
2116            e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:)
2117         ENDIF
2118         !
2119      ENDIF
2120#endif
2121
2122      fsdept(:,:,:) = gdept_0 (:,:,:)
2123      fsdepw(:,:,:) = gdepw_0 (:,:,:)
2124      fsde3w(:,:,:) = gdep3w_0(:,:,:)
2125      fse3t (:,:,:) = e3t_0   (:,:,:)
2126      fse3u (:,:,:) = e3u_0   (:,:,:)
2127      fse3v (:,:,:) = e3v_0   (:,:,:)
2128      fse3f (:,:,:) = e3f_0   (:,:,:)
2129      fse3w (:,:,:) = e3w_0   (:,:,:)
2130      fse3uw(:,:,:) = e3uw_0  (:,:,:)
2131      fse3vw(:,:,:) = e3vw_0  (:,:,:)
2132!!
2133      ! HYBRID :
2134      DO jj = 1, jpj
2135         DO ji = 1, jpi
2136            DO jk = 1, jpkm1
2137               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
2138            END DO
2139            IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0
2140         END DO
2141      END DO
2142      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
2143         &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
2144
2145      IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain
2146         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) )
2147         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
2148            &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w '  , MINVAL( gdep3w_0(:,:,:) )
2149         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0   (:,:,:) ),   &
2150            &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0   (:,:,:) ),   &
2151            &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0  (:,:,:) ),   &
2152            &                          ' w ', MINVAL( e3w_0  (:,:,:) )
2153
2154         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
2155            &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w '  , MAXVAL( gdep3w_0(:,:,:) )
2156         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0   (:,:,:) ),   &
2157            &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0   (:,:,:) ),   &
2158            &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0  (:,:,:) ),   &
2159            &                          ' w ', MAXVAL( e3w_0  (:,:,:) )
2160      ENDIF
2161      !  END DO
2162      IF(lwp) THEN                                  ! selected vertical profiles
2163         WRITE(numout,*)
2164         WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
2165         WRITE(numout,*) ' ~~~~~~  --------------------'
2166         WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2167         WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     &
2168            &                                 e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
2169         DO jj = mj0(20), mj1(20)
2170            DO ji = mi0(20), mi1(20)
2171               WRITE(numout,*)
2172               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
2173               WRITE(numout,*) ' ~~~~~~  --------------------'
2174               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2175               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
2176                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2177            END DO
2178         END DO
2179         DO jj = mj0(74), mj1(74)
2180            DO ji = mi0(100), mi1(100)
2181               WRITE(numout,*)
2182               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
2183               WRITE(numout,*) ' ~~~~~~  --------------------'
2184               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2185               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
2186                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2187            END DO
2188         END DO
2189      ENDIF
2190
2191!================================================================================
2192! check the coordinate makes sense
2193!================================================================================
2194      DO ji = 1, jpi
2195         DO jj = 1, jpj
2196
2197            IF( hbatt(ji,jj) > 0._wp) THEN
2198               DO jk = 1, mbathy(ji,jj)
2199                 ! check coordinate is monotonically increasing
2200                 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
2201                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
2202                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
2203                    WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
2204                    WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
2205                    CALL ctl_stop( ctmp1 )
2206                 ENDIF
2207                 ! and check it has never gone negative
2208                 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
2209                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
2210                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk
2211                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2212                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2213                    CALL ctl_stop( ctmp1 )
2214                 ENDIF
2215                 ! and check it never exceeds the total depth
2216                 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
2217                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
2218                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
2219                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2220                    CALL ctl_stop( ctmp1 )
2221                 ENDIF
2222               END DO
2223
2224               DO jk = 1, mbathy(ji,jj)-1
2225                 ! and check it never exceeds the total depth
2226                IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
2227                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
2228                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
2229                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2230                    CALL ctl_stop( ctmp1 )
2231                 ENDIF
2232               END DO
2233
2234            ENDIF
2235
2236         END DO
2237      END DO
2238      !
2239      CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
2240      !
2241      IF( nn_timing == 1 )  CALL timing_stop('zgr_sco')
2242      !
2243   END SUBROUTINE zgr_sco
2244
2245!!======================================================================
2246   SUBROUTINE s_sh94()
2247
2248      !!----------------------------------------------------------------------
2249      !!                  ***  ROUTINE s_sh94  ***
2250      !!                     
2251      !! ** Purpose :   stretch the s-coordinate system
2252      !!
2253      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994
2254      !!                mixed S/sigma coordinate
2255      !!
2256      !! Reference : Song and Haidvogel 1994.
2257      !!----------------------------------------------------------------------
2258      !
2259      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2260      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
2261      !
2262      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2263      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
2264
2265      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2266      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2267
2268      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
2269      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
2270      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
2271      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
2272
2273      DO ji = 1, jpi
2274         DO jj = 1, jpj
2275
2276            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
2277               DO jk = 1, jpk
2278                  z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
2279                  z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
2280               END DO
2281            ELSE ! shallow water, uniform sigma
2282               DO jk = 1, jpk
2283                  z_gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
2284                  z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
2285                  END DO
2286            ENDIF
2287            !
2288            DO jk = 1, jpkm1
2289               z_esigt3(ji,jj,jk  ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2290               z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2291            END DO
2292            z_esigw3(ji,jj,1  ) = 2._wp * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) )
2293            z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
2294            !
2295            ! Coefficients for vertical depth as the sum of e3w scale factors
2296            z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
2297            DO jk = 2, jpk
2298               z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2299            END DO
2300            !
2301            DO jk = 1, jpk
2302               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2303               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2304               gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
2305               gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
2306               gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
2307            END DO
2308           !
2309         END DO   ! for all jj's
2310      END DO    ! for all ji's
2311
2312      DO ji = 1, jpim1
2313         DO jj = 1, jpjm1
2314            DO jk = 1, jpk
2315               z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
2316                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2317               z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
2318                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2319               z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     &
2320                  &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   &
2321                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2322               z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
2323                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2324               z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
2325                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2326               !
2327               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2328               e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2329               e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2330               e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2331               !
2332               e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2333               e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2334               e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2335            END DO
2336        END DO
2337      END DO
2338
2339      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2340      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2341
2342   END SUBROUTINE s_sh94
2343
2344   SUBROUTINE s_sf12
2345
2346      !!----------------------------------------------------------------------
2347      !!                  ***  ROUTINE s_sf12 ***
2348      !!                     
2349      !! ** Purpose :   stretch the s-coordinate system
2350      !!
2351      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
2352      !!                mixed S/sigma/Z coordinate
2353      !!
2354      !!                This method allows the maintenance of fixed surface and or
2355      !!                bottom cell resolutions (cf. geopotential coordinates)
2356      !!                within an analytically derived stretched S-coordinate framework.
2357      !!
2358      !!
2359      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
2360      !!----------------------------------------------------------------------
2361      !
2362      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2363      REAL(wp) ::   zsmth               ! smoothing around critical depth
2364      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
2365      !
2366      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2367      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
2368
2369      !
2370      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2371      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2372
2373      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
2374      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
2375      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
2376      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
2377
2378      DO ji = 1, jpi
2379         DO jj = 1, jpj
2380
2381          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
2382             
2383              zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
2384                                                     ! could be changed by users but care must be taken to do so carefully
2385              zzb = 1.0_wp-(zzb/hbatt(ji,jj))
2386           
2387              zzs = rn_zs / hbatt(ji,jj) 
2388             
2389              IF (rn_efold /= 0.0_wp) THEN
2390                zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
2391              ELSE
2392                zsmth = 1.0_wp 
2393              ENDIF
2394               
2395              DO jk = 1, jpk
2396                z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
2397                z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
2398              ENDDO
2399              z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth  )
2400              z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth  )
2401 
2402          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
2403
2404            DO jk = 1, jpk
2405              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
2406              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
2407            END DO
2408
2409          ELSE  ! shallow water, z coordinates
2410
2411            DO jk = 1, jpk
2412              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2413              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2414            END DO
2415
2416          ENDIF
2417
2418          DO jk = 1, jpkm1
2419             z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2420             z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2421          END DO
2422          z_esigw3(ji,jj,1  ) = 2.0_wp * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ))
2423          z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
2424
2425          ! Coefficients for vertical depth as the sum of e3w scale factors
2426          z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
2427          DO jk = 2, jpk
2428             z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2429          END DO
2430
2431          DO jk = 1, jpk
2432             gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
2433             gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
2434             gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
2435          END DO
2436
2437        ENDDO   ! for all jj's
2438      ENDDO    ! for all ji's
2439
2440      DO ji=1,jpi-1
2441        DO jj=1,jpj-1
2442
2443          DO jk = 1, jpk
2444                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / &
2445                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2446                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / &
2447                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2448                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  &
2449                                      hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / &
2450                                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2451                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / &
2452                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2453                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / &
2454                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2455
2456             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
2457             e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
2458             e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
2459             e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
2460             !
2461             e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
2462             e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
2463             e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
2464          END DO
2465
2466        ENDDO
2467      ENDDO
2468      !
2469      CALL lbc_lnk(e3t_0 ,'T',1.) ; CALL lbc_lnk(e3u_0 ,'T',1.)
2470      CALL lbc_lnk(e3v_0 ,'T',1.) ; CALL lbc_lnk(e3f_0 ,'T',1.)
2471      CALL lbc_lnk(e3w_0 ,'T',1.)
2472      CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.)
2473      !
2474      !                                               ! =============
2475
2476      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2477      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2478
2479   END SUBROUTINE s_sf12
2480
2481   SUBROUTINE s_tanh()
2482
2483      !!----------------------------------------------------------------------
2484      !!                  ***  ROUTINE s_tanh***
2485      !!                     
2486      !! ** Purpose :   stretch the s-coordinate system
2487      !!
2488      !! ** Method  :   s-coordinate stretch
2489      !!
2490      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
2491      !!----------------------------------------------------------------------
2492
2493      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2494      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
2495
2496      REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
2497      REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw
2498
2499      CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      )
2500      CALL wrk_alloc( jpk, z_esigt, z_esigw                                               )
2501
2502      z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp   ;   z_gsi3w  = 0._wp
2503      z_esigt  = 0._wp   ;   z_esigw  = 0._wp 
2504
2505      DO jk = 1, jpk
2506        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
2507        z_gsigt(jk) = -fssig( REAL(jk,wp)        )
2508      END DO
2509      IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk)
2510      !
2511      ! Coefficients for vertical scale factors at w-, t- levels
2512!!gm bug :  define it from analytical function, not like juste bellow....
2513!!gm        or betteroffer the 2 possibilities....
2514      DO jk = 1, jpkm1
2515         z_esigt(jk  ) = z_gsigw(jk+1) - z_gsigw(jk)
2516         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
2517      END DO
2518      z_esigw( 1 ) = 2._wp * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
2519      z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
2520      !
2521      ! Coefficients for vertical depth as the sum of e3w scale factors
2522      z_gsi3w(1) = 0.5_wp * z_esigw(1)
2523      DO jk = 2, jpk
2524         z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
2525      END DO
2526!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
2527      DO jk = 1, jpk
2528         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2529         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2530         gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
2531         gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
2532         gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )
2533      END DO
2534!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
2535      DO jj = 1, jpj
2536         DO ji = 1, jpi
2537            DO jk = 1, jpk
2538              e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2539              e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2540              e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2541              e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
2542              !
2543              e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2544              e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2545              e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2546            END DO
2547         END DO
2548      END DO
2549
2550      CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      )
2551      CALL wrk_dealloc( jpk, z_esigt, z_esigw                                               )
2552
2553   END SUBROUTINE s_tanh
2554
2555   FUNCTION fssig( pk ) RESULT( pf )
2556      !!----------------------------------------------------------------------
2557      !!                 ***  ROUTINE fssig ***
2558      !!       
2559      !! ** Purpose :   provide the analytical function in s-coordinate
2560      !!         
2561      !! ** Method  :   the function provide the non-dimensional position of
2562      !!                T and W (i.e. between 0 and 1)
2563      !!                T-points at integer values (between 1 and jpk)
2564      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2565      !!----------------------------------------------------------------------
2566      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
2567      REAL(wp)             ::   pf   ! sigma value
2568      !!----------------------------------------------------------------------
2569      !
2570      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   &
2571         &     - TANH( rn_thetb * rn_theta                                )  )   &
2572         & * (   COSH( rn_theta                           )                      &
2573         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              &
2574         & / ( 2._wp * SINH( rn_theta ) )
2575      !
2576   END FUNCTION fssig
2577
2578
2579   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
2580      !!----------------------------------------------------------------------
2581      !!                 ***  ROUTINE fssig1 ***
2582      !!
2583      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
2584      !!
2585      !! ** Method  :   the function provides the non-dimensional position of
2586      !!                T and W (i.e. between 0 and 1)
2587      !!                T-points at integer values (between 1 and jpk)
2588      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2589      !!----------------------------------------------------------------------
2590      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
2591      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
2592      REAL(wp)             ::   pf1   ! sigma value
2593      !!----------------------------------------------------------------------
2594      !
2595      IF ( rn_theta == 0 ) then      ! uniform sigma
2596         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
2597      ELSE                        ! stretched sigma
2598         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
2599            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
2600            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
2601      ENDIF
2602      !
2603   END FUNCTION fssig1
2604
2605
2606   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
2607      !!----------------------------------------------------------------------
2608      !!                 ***  ROUTINE fgamma  ***
2609      !!
2610      !! ** Purpose :   provide analytical function for the s-coordinate
2611      !!
2612      !! ** Method  :   the function provides the non-dimensional position of
2613      !!                T and W (i.e. between 0 and 1)
2614      !!                T-points at integer values (between 1 and jpk)
2615      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2616      !!
2617      !!                This method allows the maintenance of fixed surface and or
2618      !!                bottom cell resolutions (cf. geopotential coordinates)
2619      !!                within an analytically derived stretched S-coordinate framework.
2620      !!
2621      !! Reference  :   Siddorn and Furner, in prep
2622      !!----------------------------------------------------------------------
2623      REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate
2624      REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate
2625      REAL(wp), INTENT(in   ) ::   pzb           ! Bottom box depth
2626      REAL(wp), INTENT(in   ) ::   pzs           ! surface box depth
2627      REAL(wp), INTENT(in   ) ::   psmth       ! Smoothing parameter
2628      REAL(wp)                ::   za1,za2,za3    ! local variables
2629      REAL(wp)                ::   zn1,zn2        ! local variables
2630      REAL(wp)                ::   za,zb,zx       ! local variables
2631      integer                 ::   jk
2632      !!----------------------------------------------------------------------
2633      !
2634
2635      zn1  =  1./(jpk-1.)
2636      zn2  =  1. -  zn1
2637
2638      za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 
2639      za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp)
2640      za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
2641     
2642      za = pzb - za3*(pzs-za1)-za2
2643      za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) )
2644      zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
2645      zx = 1.0_wp-za/2.0_wp-zb
2646 
2647      DO jk = 1, jpk
2648        p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp +  &
2649                    & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &
2650                    &      (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
2651        p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth)
2652      ENDDO 
2653
2654      !
2655   END FUNCTION fgamma
2656
2657   !!======================================================================
2658END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.