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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domzgr.F90 in branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5216

Last change on this file since 5216 was 5216, checked in by mathiot, 9 years ago

ISF branch: minor changes in zpshde, sbcisf + bug correction in domzgr (e3wu) only if ice shelf

  • Property svn:keywords set to Id
File size: 133.8 KB
Line 
1MODULE domzgr
2   !!==============================================================================
3   !!                       ***  MODULE domzgr   ***
4   !! Ocean initialization : domain initialization
5   !!==============================================================================
6   !! History :  OPA  ! 1995-12  (G. Madec)  Original code : s vertical coordinate
7   !!                 ! 1997-07  (G. Madec)  lbc_lnk call
8   !!                 ! 1997-04  (J.-O. Beismann)
9   !!            8.5  ! 2002-09  (A. Bozec, G. Madec)  F90: Free form and module
10   !!             -   ! 2002-09  (A. de Miranda)  rigid-lid + islands
11   !!  NEMO      1.0  ! 2003-08  (G. Madec)  F90: Free form and module
12   !!             -   ! 2005-10  (A. Beckmann)  modifications for hybrid s-ccordinates & new stretching function
13   !!            2.0  ! 2006-04  (R. Benshila, G. Madec)  add zgr_zco
14   !!            3.0  ! 2008-06  (G. Madec)  insertion of domzgr_zps.h90 & conding style
15   !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option
16   !!            3.3  ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level
17   !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn and Furner stretching function
18   !!            3.4  ! 2012-12  (R. Bourdalle-Badie and G. Reffray)  modify C1D case 
19   !!            3.6  ! 2014-11  (P. Mathiot and C. Harris) add ice shelf capabilitye 
20   !!----------------------------------------------------------------------
21
22   !!----------------------------------------------------------------------
23   !!   dom_zgr          : defined the ocean vertical coordinate system
24   !!       zgr_bat      : bathymetry fields (levels and meters)
25   !!       zgr_bat_zoom : modify the bathymetry field if zoom domain
26   !!       zgr_bat_ctl  : check the bathymetry files
27   !!       zgr_bot_level: deepest ocean level for t-, u, and v-points
28   !!       zgr_z        : reference z-coordinate
29   !!       zgr_zco      : z-coordinate
30   !!       zgr_zps      : z-coordinate with partial steps
31   !!       zgr_sco      : s-coordinate
32   !!       fssig        : tanh stretch function
33   !!       fssig1       : Song and Haidvogel 1994 stretch function
34   !!       fgamma       : Siddorn and Furner 2012 stretching function
35   !!---------------------------------------------------------------------
36   USE oce               ! ocean variables
37   USE dom_oce           ! ocean domain
38   USE closea            ! closed seas
39   USE c1d               ! 1D vertical configuration
40   USE in_out_manager    ! I/O manager
41   USE iom               ! I/O library
42   USE lbclnk            ! ocean lateral boundary conditions (or mpp link)
43   USE lib_mpp           ! distributed memory computing library
44   USE wrk_nemo          ! Memory allocation
45   USE timing            ! Timing
46
47   IMPLICIT NONE
48   PRIVATE
49
50   PUBLIC   dom_zgr        ! called by dom_init.F90
51
52   !                              !!* 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, 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            IF ( ln_isfcav ) THEN
538               CALL iom_get  ( inum, jpdom_data, 'Bathymetry_isf', bathy, lrowattr=.false. )
539            ELSE
540               CALL iom_get  ( inum, jpdom_data, 'Bathymetry'    , bathy, lrowattr=ln_use_jattr  )
541            END IF
542            CALL iom_close( inum )
543            !                                               
544            risfdep(:,:)=0._wp         
545            misfdep(:,:)=1             
546            IF ( ln_isfcav ) THEN
547               CALL iom_open ( 'isf_draft_meter.nc', inum ) 
548               CALL iom_get  ( inum, jpdom_data, 'isf_draft', risfdep )
549               CALL iom_close( inum )
550               WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp
551            END IF
552            !       
553            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
554               !
555              IF( nn_cla == 0 ) THEN
556                 ii0 = 140   ;   ii1 = 140                   ! Gibraltar Strait open
557                 ij0 = 102   ;   ij1 = 102                   ! (Thomson, Ocean Modelling, 1995)
558                 DO ji = mi0(ii0), mi1(ii1)
559                    DO jj = mj0(ij0), mj1(ij1)
560                       bathy(ji,jj) = 284._wp
561                    END DO
562                 END DO
563                 IF(lwp) WRITE(numout,*)     
564                 IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
565                 !
566                 ii0 = 160   ;   ii1 = 160                   ! Bab el mandeb Strait open
567                 ij0 = 88    ;   ij1 = 88                    ! (Thomson, Ocean Modelling, 1995)
568                 DO ji = mi0(ii0), mi1(ii1)
569                    DO jj = mj0(ij0), mj1(ij1)
570                       bathy(ji,jj) = 137._wp
571                    END DO
572                 END DO
573                 IF(lwp) WRITE(numout,*)
574                 IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
575              ENDIF
576              !
577           ENDIF
578            !
579        ENDIF
580         !                                            ! =============== !
581      ELSE                                            !      error      !
582         !                                            ! =============== !
583         WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
584         CALL ctl_stop( '    zgr_bat : '//trim(ctmp1) )
585      ENDIF
586      !
587      IF( nn_closea == 0 )   CALL clo_bat( bathy, mbathy )    !==  NO closed seas or lakes  ==!
588      !                       
589      IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==!
590         ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin
591         IF ( ln_isfcav ) THEN
592            WHERE (bathy == risfdep)
593               bathy   = 0.0_wp ; risfdep = 0.0_wp
594            END WHERE
595         END IF
596         ! end patch
597         IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level
598         ELSE                          ;   ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 )  ! from a depth
599         ENDIF
600         zhmin = gdepw_1d(ik+1)                                                         ! minimum depth = ik+1 w-levels
601         WHERE( bathy(:,:) <= 0._wp )   ;   bathy(:,:) = 0._wp                         ! min=0     over the lands
602         ELSE WHERE                     ;   bathy(:,:) = MAX(  zhmin , bathy(:,:)  )   ! min=zhmin over the oceans
603         END WHERE
604         IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
605      ENDIF
606      !
607      IF( nn_timing == 1 )  CALL timing_stop('zgr_bat')
608      !
609   END SUBROUTINE zgr_bat
610
611
612   SUBROUTINE zgr_bat_zoom
613      !!----------------------------------------------------------------------
614      !!                    ***  ROUTINE zgr_bat_zoom  ***
615      !!
616      !! ** Purpose : - Close zoom domain boundary if necessary
617      !!              - Suppress Med Sea from ORCA R2 and R05 arctic zoom
618      !!
619      !! ** Method  :
620      !!
621      !! ** Action  : - update mbathy: level bathymetry (in level index)
622      !!----------------------------------------------------------------------
623      INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integers
624      !!----------------------------------------------------------------------
625      !
626      IF(lwp) WRITE(numout,*)
627      IF(lwp) WRITE(numout,*) '    zgr_bat_zoom : modify the level bathymetry for zoom domain'
628      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~'
629      !
630      ! Zoom domain
631      ! ===========
632      !
633      ! Forced closed boundary if required
634      IF( lzoom_s )   mbathy(  : , mj0(jpjzoom):mj1(jpjzoom) )      = 0
635      IF( lzoom_w )   mbathy(      mi0(jpizoom):mi1(jpizoom) , :  ) = 0
636      IF( lzoom_e )   mbathy(      mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , :  ) = 0
637      IF( lzoom_n )   mbathy(  : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) )      = 0
638      !
639      ! Configuration specific domain modifications
640      ! (here, ORCA arctic configuration: suppress Med Sea)
641      IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN
642         SELECT CASE ( jp_cfg )
643         !                                        ! =======================
644         CASE ( 2 )                               !  ORCA_R2 configuration
645            !                                     ! =======================
646            IF(lwp) WRITE(numout,*) '                   ORCA R2 arctic zoom: suppress the Med Sea'
647            ii0 = 141   ;   ii1 = 162      ! Sea box i,j indices
648            ij0 =  98   ;   ij1 = 110
649            !                                     ! =======================
650         CASE ( 05 )                              !  ORCA_R05 configuration
651            !                                     ! =======================
652            IF(lwp) WRITE(numout,*) '                   ORCA R05 arctic zoom: suppress the Med Sea'
653            ii0 = 563   ;   ii1 = 642      ! zero over the Med Sea boxe
654            ij0 = 314   ;   ij1 = 370 
655         END SELECT
656         !
657         mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0   ! zero over the Med Sea boxe
658         !
659      ENDIF
660      !
661   END SUBROUTINE zgr_bat_zoom
662
663
664   SUBROUTINE zgr_bat_ctl
665      !!----------------------------------------------------------------------
666      !!                    ***  ROUTINE zgr_bat_ctl  ***
667      !!
668      !! ** Purpose :   check the bathymetry in levels
669      !!
670      !! ** Method  :   The array mbathy is checked to verified its consistency
671      !!      with the model options. in particular:
672      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
673      !!                  along closed boundary.
674      !!            mbathy must be cyclic IF jperio=1.
675      !!            mbathy must be lower or equal to jpk-1.
676      !!            isolated ocean grid points are suppressed from mbathy
677      !!                  since they are only connected to remaining
678      !!                  ocean through vertical diffusion.
679      !!      C A U T I O N : mbathy will be modified during the initializa-
680      !!      tion phase to become the number of non-zero w-levels of a water
681      !!      column, with a minimum value of 1.
682      !!
683      !! ** Action  : - update mbathy: level bathymetry (in level index)
684      !!              - update bathy : meter bathymetry (in meters)
685      !!----------------------------------------------------------------------
686      !!
687      INTEGER ::   ji, jj, jl                    ! dummy loop indices
688      INTEGER ::   icompt, ibtest, ikmax         ! temporary integers
689      REAL(wp), POINTER, DIMENSION(:,:) ::  zbathy
690
691      !!----------------------------------------------------------------------
692      !
693      IF( nn_timing == 1 )  CALL timing_start('zgr_bat_ctl')
694      !
695      CALL wrk_alloc( jpi, jpj, zbathy )
696      !
697      IF(lwp) WRITE(numout,*)
698      IF(lwp) WRITE(numout,*) '    zgr_bat_ctl : check the bathymetry'
699      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
700      !                                          ! Suppress isolated ocean grid points
701      IF(lwp) WRITE(numout,*)
702      IF(lwp) WRITE(numout,*)'                   suppress isolated ocean grid points'
703      IF(lwp) WRITE(numout,*)'                   -----------------------------------'
704      icompt = 0
705      DO jl = 1, 2
706         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
707            mbathy( 1 ,:) = mbathy(jpim1,:)           ! local domain is cyclic east-west
708            mbathy(jpi,:) = mbathy(  2  ,:)
709         ENDIF
710         DO jj = 2, jpjm1
711            DO ji = 2, jpim1
712               ibtest = MAX(  mbathy(ji-1,jj), mbathy(ji+1,jj),   &
713                  &           mbathy(ji,jj-1), mbathy(ji,jj+1)  )
714               IF( ibtest < mbathy(ji,jj) ) THEN
715                  IF(lwp) WRITE(numout,*) ' the number of ocean level at ',   &
716                     &   'grid-point (i,j) =  ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
717                  mbathy(ji,jj) = ibtest
718                  icompt = icompt + 1
719               ENDIF
720            END DO
721         END DO
722      END DO
723      IF( lk_mpp )   CALL mpp_sum( icompt )
724      IF( icompt == 0 ) THEN
725         IF(lwp) WRITE(numout,*)'     no isolated ocean grid points'
726      ELSE
727         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed'
728      ENDIF
729      IF( lk_mpp ) THEN
730         zbathy(:,:) = FLOAT( mbathy(:,:) )
731         CALL lbc_lnk( zbathy, 'T', 1._wp )
732         mbathy(:,:) = INT( zbathy(:,:) )
733      ENDIF
734      !                                          ! East-west cyclic boundary conditions
735      IF( nperio == 0 ) THEN
736         IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio
737         IF( lk_mpp ) THEN
738            IF( nbondi == -1 .OR. nbondi == 2 ) THEN
739               IF( jperio /= 1 )   mbathy(1,:) = 0
740            ENDIF
741            IF( nbondi == 1 .OR. nbondi == 2 ) THEN
742               IF( jperio /= 1 )   mbathy(nlci,:) = 0
743            ENDIF
744         ELSE
745            IF( ln_zco .OR. ln_zps ) THEN
746               mbathy( 1 ,:) = 0
747               mbathy(jpi,:) = 0
748            ELSE
749               mbathy( 1 ,:) = jpkm1
750               mbathy(jpi,:) = jpkm1
751            ENDIF
752         ENDIF
753      ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio ==  6 ) THEN
754         IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio
755         mbathy( 1 ,:) = mbathy(jpim1,:)
756         mbathy(jpi,:) = mbathy(  2  ,:)
757      ELSEIF( nperio == 2 ) THEN
758         IF(lwp) WRITE(numout,*) '   equatorial boundary conditions on mbathy: nperio = ', nperio
759      ELSE
760         IF(lwp) WRITE(numout,*) '    e r r o r'
761         IF(lwp) WRITE(numout,*) '    parameter , nperio = ', nperio
762         !         STOP 'dom_mba'
763      ENDIF
764      !  Boundary condition on mbathy
765      IF( .NOT.lk_mpp ) THEN 
766!!gm     !!bug ???  think about it !
767         !   ... mono- or macro-tasking: T-point, >0, 2D array, no slab
768         zbathy(:,:) = FLOAT( mbathy(:,:) )
769         CALL lbc_lnk( zbathy, 'T', 1._wp )
770         mbathy(:,:) = INT( zbathy(:,:) )
771      ENDIF
772      ! Number of ocean level inferior or equal to jpkm1
773      ikmax = 0
774      DO jj = 1, jpj
775         DO ji = 1, jpi
776            ikmax = MAX( ikmax, mbathy(ji,jj) )
777         END DO
778      END DO
779!!gm  !!! test to do:   ikmax = MAX( mbathy(:,:) )   ???
780      IF( ikmax > jpkm1 ) THEN
781         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' >  jpk-1'
782         IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
783      ELSE IF( ikmax < jpkm1 ) THEN
784         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1' 
785         IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
786      ENDIF
787
788      IF( lwp .AND. nprint == 1 ) THEN      ! control print
789         WRITE(numout,*)
790         WRITE(numout,*) ' bathymetric field :   number of non-zero T-levels '
791         WRITE(numout,*) ' ------------------'
792         CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )
793         WRITE(numout,*)
794      ENDIF
795      !
796      CALL wrk_dealloc( jpi, jpj, zbathy )
797      !
798      IF( nn_timing == 1 )  CALL timing_stop('zgr_bat_ctl')
799      !
800   END SUBROUTINE zgr_bat_ctl
801
802
803   SUBROUTINE zgr_bot_level
804      !!----------------------------------------------------------------------
805      !!                    ***  ROUTINE zgr_bot_level  ***
806      !!
807      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
808      !!
809      !! ** Method  :   computes from mbathy with a minimum value of 1 over land
810      !!
811      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest
812      !!                                     ocean level at t-, u- & v-points
813      !!                                     (min value = 1 over land)
814      !!----------------------------------------------------------------------
815      !!
816      INTEGER ::   ji, jj   ! dummy loop indices
817      REAL(wp), POINTER, DIMENSION(:,:) ::  zmbk
818      !!----------------------------------------------------------------------
819      !
820      IF( nn_timing == 1 )  CALL timing_start('zgr_bot_level')
821      !
822      CALL wrk_alloc( jpi, jpj, zmbk )
823      !
824      IF(lwp) WRITE(numout,*)
825      IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
826      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
827      !
828      mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
829 
830      !                                     ! bottom k-index of W-level = mbkt+1
831      DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level
832         DO ji = 1, jpim1
833            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
834            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
835         END DO
836      END DO
837      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
838      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
839      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
840      !
841      CALL wrk_dealloc( jpi, jpj, zmbk )
842      !
843      IF( nn_timing == 1 )  CALL timing_stop('zgr_bot_level')
844      !
845   END SUBROUTINE zgr_bot_level
846
847      SUBROUTINE zgr_top_level
848      !!----------------------------------------------------------------------
849      !!                    ***  ROUTINE zgr_bot_level  ***
850      !!
851      !! ** Purpose :   defines the vertical index of ocean top (mik. arrays)
852      !!
853      !! ** Method  :   computes from misfdep with a minimum value of 1
854      !!
855      !! ** Action  :   mikt, miku, mikv :   vertical indices of the shallowest
856      !!                                     ocean level at t-, u- & v-points
857      !!                                     (min value = 1)
858      !!----------------------------------------------------------------------
859      !!
860      INTEGER ::   ji, jj   ! dummy loop indices
861      REAL(wp), POINTER, DIMENSION(:,:) ::  zmik
862      !!----------------------------------------------------------------------
863      !
864      IF( nn_timing == 1 )  CALL timing_start('zgr_top_level')
865      !
866      CALL wrk_alloc( jpi, jpj, zmik )
867      !
868      IF(lwp) WRITE(numout,*)
869      IF(lwp) WRITE(numout,*) '    zgr_top_level : ocean top k-index of T-, U-, V- and W-levels '
870      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
871      !
872      mikt(:,:) = MAX( misfdep(:,:) , 1 )    ! top k-index of T-level (=1)
873      !                                      ! top k-index of W-level (=mikt)
874      DO jj = 1, jpjm1                       ! top k-index of U- (U-) level
875         DO ji = 1, jpim1
876            miku(ji,jj) = MAX(  mikt(ji+1,jj  ) , mikt(ji,jj)  )
877            mikv(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj)  )
878            mikf(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj), mikt(ji+1,jj  ), mikt(ji+1,jj+1)  )
879         END DO
880      END DO
881
882      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
883      zmik(:,:) = REAL( miku(:,:), wp )   ;   CALL lbc_lnk(zmik,'U',1.)   ;   miku  (:,:) = MAX( INT( zmik(:,:) ), 1 )
884      zmik(:,:) = REAL( mikv(:,:), wp )   ;   CALL lbc_lnk(zmik,'V',1.)   ;   mikv  (:,:) = MAX( INT( zmik(:,:) ), 1 )
885      zmik(:,:) = REAL( mikf(:,:), wp )   ;   CALL lbc_lnk(zmik,'F',1.)   ;   mikf  (:,:) = MAX( INT( zmik(:,:) ), 1 )
886      !
887      CALL wrk_dealloc( jpi, jpj, zmik )
888      !
889      IF( nn_timing == 1 )  CALL timing_stop('zgr_top_level')
890      !
891   END SUBROUTINE zgr_top_level
892
893   SUBROUTINE zgr_zco
894      !!----------------------------------------------------------------------
895      !!                  ***  ROUTINE zgr_zco  ***
896      !!
897      !! ** Purpose :   define the z-coordinate system
898      !!
899      !! ** Method  :   set 3D coord. arrays to reference 1D array
900      !!----------------------------------------------------------------------
901      INTEGER  ::   jk
902      !!----------------------------------------------------------------------
903      !
904      IF( nn_timing == 1 )  CALL timing_start('zgr_zco')
905      !
906      DO jk = 1, jpk
907         gdept_0 (:,:,jk) = gdept_1d(jk)
908         gdepw_0 (:,:,jk) = gdepw_1d(jk)
909         gdep3w_0(:,:,jk) = gdepw_1d(jk)
910         e3t_0   (:,:,jk) = e3t_1d  (jk)
911         e3u_0   (:,:,jk) = e3t_1d  (jk)
912         e3v_0   (:,:,jk) = e3t_1d  (jk)
913         e3f_0   (:,:,jk) = e3t_1d  (jk)
914         e3w_0   (:,:,jk) = e3w_1d  (jk)
915         e3uw_0  (:,:,jk) = e3w_1d  (jk)
916         e3vw_0  (:,:,jk) = e3w_1d  (jk)
917      END DO
918      !
919      IF( nn_timing == 1 )  CALL timing_stop('zgr_zco')
920      !
921   END SUBROUTINE zgr_zco
922
923
924   SUBROUTINE zgr_zps
925      !!----------------------------------------------------------------------
926      !!                  ***  ROUTINE zgr_zps  ***
927      !!                     
928      !! ** Purpose :   the depth and vertical scale factor in partial step
929      !!      z-coordinate case
930      !!
931      !! ** Method  :   Partial steps : computes the 3D vertical scale factors
932      !!      of T-, U-, V-, W-, UW-, VW and F-points that are associated with
933      !!      a partial step representation of bottom topography.
934      !!
935      !!        The reference depth of model levels is defined from an analytical
936      !!      function the derivative of which gives the reference vertical
937      !!      scale factors.
938      !!        From  depth and scale factors reference, we compute there new value
939      !!      with partial steps  on 3d arrays ( i, j, k ).
940      !!
941      !!              w-level: gdepw_0(i,j,k)  = gdep(k)
942      !!                       e3w_0(i,j,k) = dk(gdep)(k)     = e3(i,j,k)
943      !!              t-level: gdept_0(i,j,k)  = gdep(k+0.5)
944      !!                       e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5)
945      !!
946      !!        With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
947      !!      we find the mbathy index of the depth at each grid point.
948      !!      This leads us to three cases:
949      !!
950      !!              - bathy = 0 => mbathy = 0
951      !!              - 1 < mbathy < jpkm1   
952      !!              - bathy > gdepw_0(jpk) => mbathy = jpkm1 
953      !!
954      !!        Then, for each case, we find the new depth at t- and w- levels
955      !!      and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
956      !!      and f-points.
957      !!
958      !!        This routine is given as an example, it must be modified
959      !!      following the user s desiderata. nevertheless, the output as
960      !!      well as the way to compute the model levels and scale factors
961      !!      must be respected in order to insure second order accuracy
962      !!      schemes.
963      !!
964      !!         c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives
965      !!         - - - - - - -   gdept_0, gdepw_0 and e3. are positives
966      !!     
967      !!  Reference :   Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
968      !!----------------------------------------------------------------------
969      !!
970      INTEGER  ::   ji, jj, jk       ! dummy loop indices
971      INTEGER  ::   ik, it, ikb, ikt ! temporary integers
972      LOGICAL  ::   ll_print         ! Allow  control print for debugging
973      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
974      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
975      REAL(wp) ::   zdiff            ! temporary scalar
976      REAL(wp) ::   zmax             ! temporary scalar
977      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt
978      !!---------------------------------------------------------------------
979      !
980      IF( nn_timing == 1 )  CALL timing_start('zgr_zps')
981      !
982      CALL wrk_alloc( jpi, jpj, jpk, zprt )
983      !
984      IF(lwp) WRITE(numout,*)
985      IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps'
986      IF(lwp) WRITE(numout,*) '    ~~~~~~~ '
987      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used'
988
989      ll_print = .FALSE.                   ! Local variable for debugging
990     
991      IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth
992         WRITE(numout,*)
993         WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)'
994         CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
995      ENDIF
996
997
998      ! bathymetry in level (from bathy_meter)
999      ! ===================
1000      zmax = gdepw_1d(jpk) + e3t_1d(jpk)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) )
1001      bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat)
1002      WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0
1003      ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level
1004      END WHERE
1005
1006      ! Compute mbathy for ocean points (i.e. the number of ocean levels)
1007      ! find the number of ocean levels such that the last level thickness
1008      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where
1009      ! e3t_1d is the reference level thickness
1010      DO jk = jpkm1, 1, -1
1011         zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1012         WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1
1013      END DO
1014
1015      ! Scale factors and depth at T- and W-points
1016      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate
1017         gdept_0(:,:,jk) = gdept_1d(jk)
1018         gdepw_0(:,:,jk) = gdepw_1d(jk)
1019         e3t_0  (:,:,jk) = e3t_1d  (jk)
1020         e3w_0  (:,:,jk) = e3w_1d  (jk)
1021      END DO
1022     
1023      ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf
1024      IF ( ln_isfcav ) CALL zgr_isf
1025
1026      ! Scale factors and depth at T- and W-points
1027      IF ( .NOT. ln_isfcav ) THEN
1028         DO jj = 1, jpj
1029            DO ji = 1, jpi
1030               ik = mbathy(ji,jj)
1031               IF( ik > 0 ) THEN               ! ocean point only
1032                  ! max ocean level case
1033                  IF( ik == jpkm1 ) THEN
1034                     zdepwp = bathy(ji,jj)
1035                     ze3tp  = bathy(ji,jj) - gdepw_1d(ik)
1036                     ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
1037                     e3t_0(ji,jj,ik  ) = ze3tp
1038                     e3t_0(ji,jj,ik+1) = ze3tp
1039                     e3w_0(ji,jj,ik  ) = ze3wp
1040                     e3w_0(ji,jj,ik+1) = ze3tp
1041                     gdepw_0(ji,jj,ik+1) = zdepwp
1042                     gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp
1043                     gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
1044                     !
1045                  ELSE                         ! standard case
1046                     IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
1047                     ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
1048                     ENDIF
1049   !gm Bug?  check the gdepw_1d
1050                     !       ... on ik
1051                     gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   &
1052                        &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   &
1053                        &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) ))
1054                     e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   & 
1055                        &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ) 
1056                     e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   &
1057                        &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
1058                     !       ... on ik+1
1059                     e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
1060                     e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
1061                     gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
1062                  ENDIF
1063               ENDIF
1064            END DO
1065         END DO
1066         !
1067         it = 0
1068         DO jj = 1, jpj
1069            DO ji = 1, jpi
1070               ik = mbathy(ji,jj)
1071               IF( ik > 0 ) THEN               ! ocean point only
1072                  e3tp (ji,jj) = e3t_0(ji,jj,ik)
1073                  e3wp (ji,jj) = e3w_0(ji,jj,ik)
1074                  ! test
1075                  zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  )
1076                  IF( zdiff <= 0._wp .AND. lwp ) THEN
1077                     it = it + 1
1078                     WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj
1079                     WRITE(numout,*) ' bathy = ', bathy(ji,jj)
1080                     WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
1081                     WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  )
1082                  ENDIF
1083               ENDIF
1084            END DO
1085         END DO
1086      END IF
1087      !
1088      ! Scale factors and depth at U-, V-, UW and VW-points
1089      DO jk = 1, jpk                        ! initialisation to z-scale factors
1090         e3u_0 (:,:,jk) = e3t_1d(jk)
1091         e3v_0 (:,:,jk) = e3t_1d(jk)
1092         e3uw_0(:,:,jk) = e3w_1d(jk)
1093         e3vw_0(:,:,jk) = e3w_1d(jk)
1094      END DO
1095
1096      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
1097         DO jj = 1, jpjm1
1098            DO ji = 1, fs_jpim1   ! vector opt.
1099               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
1100               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
1101               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
1102               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
1103            END DO
1104         END DO
1105      END DO
1106      IF ( ln_isfcav ) THEN
1107      ! (ISF) define e3uw (adapted for 2 cells in the water column)
1108         DO jj = 2, jpjm1 
1109            DO ji = 2, fs_jpim1   ! vector opt.
1110               ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))
1111               ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))
1112               IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) &
1113                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) )
1114               ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))
1115               ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))
1116               IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) &
1117                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) )
1118            END DO
1119         END DO
1120      END IF
1121
1122      CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions
1123      CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp )
1124      !
1125
1126      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
1127         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk)
1128         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk)
1129         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk)
1130         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk)
1131      END DO
1132     
1133      ! Scale factor at F-point
1134      DO jk = 1, jpk                        ! initialisation to z-scale factors
1135         e3f_0(:,:,jk) = e3t_1d(jk)
1136      END DO
1137      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors
1138         DO jj = 1, jpjm1
1139            DO ji = 1, fs_jpim1   ! vector opt.
1140               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
1141            END DO
1142         END DO
1143      END DO
1144      CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions
1145      !
1146      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
1147         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk)
1148      END DO
1149!!gm  bug ? :  must be a do loop with mj0,mj1
1150      !
1151      e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2
1152      e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 
1153      e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 
1154      e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 
1155      e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 
1156
1157      ! Control of the sign
1158      IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' )
1159      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' )
1160      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' )
1161      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' )
1162     
1163      ! Compute gdep3w_0 (vertical sum of e3w)
1164      IF ( ln_isfcav ) THEN ! if cavity
1165         WHERE (misfdep == 0) misfdep = 1
1166         DO jj = 1,jpj
1167            DO ji = 1,jpi
1168               gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)
1169               DO jk = 2, misfdep(ji,jj)
1170                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
1171               END DO
1172               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))
1173               DO jk = misfdep(ji,jj) + 1, jpk
1174                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
1175               END DO
1176            END DO
1177         END DO
1178      ELSE ! no cavity
1179         gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
1180         DO jk = 2, jpk
1181            gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk)
1182         END DO
1183      END IF
1184      !                                               ! ================= !
1185      IF(lwp .AND. ll_print) THEN                     !   Control print   !
1186         !                                            ! ================= !
1187         DO jj = 1,jpj
1188            DO ji = 1, jpi
1189               ik = MAX( mbathy(ji,jj), 1 )
1190               zprt(ji,jj,1) = e3t_0   (ji,jj,ik)
1191               zprt(ji,jj,2) = e3w_0   (ji,jj,ik)
1192               zprt(ji,jj,3) = e3u_0   (ji,jj,ik)
1193               zprt(ji,jj,4) = e3v_0   (ji,jj,ik)
1194               zprt(ji,jj,5) = e3f_0   (ji,jj,ik)
1195               zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)
1196            END DO
1197         END DO
1198         WRITE(numout,*)
1199         WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1200         WRITE(numout,*)
1201         WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1202         WRITE(numout,*)
1203         WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1204         WRITE(numout,*)
1205         WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1206         WRITE(numout,*)
1207         WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1208         WRITE(numout,*)
1209         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1210      ENDIF 
1211      !
1212      CALL wrk_dealloc( jpi, jpj, jpk, zprt )
1213      !
1214      IF( nn_timing == 1 )  CALL timing_stop('zgr_zps')
1215      !
1216   END SUBROUTINE zgr_zps
1217
1218   SUBROUTINE zgr_isf
1219      !!----------------------------------------------------------------------
1220      !!                    ***  ROUTINE zgr_isf  ***
1221      !!   
1222      !! ** Purpose :   check the bathymetry in levels
1223      !!   
1224      !! ** Method  :   THe water column have to contained at least 2 cells
1225      !!                Bathymetry and isfdraft are modified (dig/close) to respect
1226      !!                this criterion.
1227      !!                 
1228      !!   
1229      !! ** Action  : - test compatibility between isfdraft and bathy
1230      !!              - bathy and isfdraft are modified
1231      !!----------------------------------------------------------------------
1232      !!   
1233      INTEGER  ::   ji, jj, jl, jk       ! dummy loop indices
1234      INTEGER  ::   ik, it               ! temporary integers
1235      INTEGER  ::   icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1   ! (ISF)
1236      REAL(wp) ::   zdepth           ! Ajusted ocean depth to avoid too small e3t
1237      REAL(wp) ::   zmax             ! Maximum and minimum depth
1238      REAL(wp) ::   zbathydiff, zrisfdepdiff  ! isf temporary scalar
1239      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
1240      REAL(wp) ::   zdepwp           ! Ajusted ocean depth to avoid too small e3t
1241      REAL(wp) ::   zdiff            ! temporary scalar
1242      REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 2D workspace (ISH)
1243      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep         ! 2D workspace (ISH)
1244      !!---------------------------------------------------------------------
1245      !
1246      IF( nn_timing == 1 )  CALL timing_start('zgr_isf')
1247      !
1248      CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep)
1249      CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy )
1250
1251
1252      ! (ISF) compute misfdep
1253      WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1 
1254      ELSEWHERE                      ;                          misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level
1255      END WHERE 
1256
1257      ! Compute misfdep for ocean points (i.e. first wet level)
1258      ! find the first ocean level such that the first level thickness
1259      ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where
1260      ! e3t_0 is the reference level thickness
1261      DO jk = 2, jpkm1 
1262         zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
1263         WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1 
1264      END DO
1265      WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp)
1266         risfdep(:,:) = 0. ; misfdep(:,:) = 1
1267      END WHERE
1268
1269      ! remove very shallow ice shelf (less than ~ 10m if 75L)
1270      IF ( cp_cfg .NE. "isomip" ) THEN
1271         WHERE (misfdep(:,:) <= 10 .AND. misfdep(:,:) .GT. 1)
1272            misfdep = 0; risfdep = 0.0_wp;
1273            mbathy  = 0; bathy   = 0.0_wp;
1274         END WHERE
1275         WHERE (bathy(:,:) <= 30.0 .AND. gphit < -60)
1276            misfdep = 0; risfdep = 0.0_wp;
1277            mbathy  = 0; bathy   = 0.0_wp;
1278         END WHERE
1279      END IF
1280 
1281! 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
1282      icompt = 0 
1283! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together
1284      DO jl = 1, 10     
1285         WHERE (bathy(:,:) .EQ. risfdep(:,:) )
1286            misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp
1287            mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp
1288         END WHERE
1289         WHERE (mbathy(:,:) <= 0) 
1290            misfdep(:,:) = 0; risfdep(:,:) = 0._wp 
1291            mbathy (:,:) = 0; bathy  (:,:) = 0._wp
1292         ENDWHERE
1293         IF( lk_mpp ) THEN
1294            zbathy(:,:) = FLOAT( misfdep(:,:) )
1295            CALL lbc_lnk( zbathy, 'T', 1. )
1296            misfdep(:,:) = INT( zbathy(:,:) )
1297            CALL lbc_lnk( risfdep, 'T', 1. )
1298            CALL lbc_lnk( bathy, 'T', 1. )
1299            zbathy(:,:) = FLOAT( mbathy(:,:) )
1300            CALL lbc_lnk( zbathy, 'T', 1. )
1301            mbathy(:,:) = INT( zbathy(:,:) )
1302         ENDIF
1303         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
1304            misfdep( 1 ,:) = misfdep(jpim1,:)           ! local domain is cyclic east-west
1305            misfdep(jpi,:) = misfdep(  2  ,:) 
1306         ENDIF
1307
1308         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
1309            mbathy( 1 ,:) = mbathy(jpim1,:)             ! local domain is cyclic east-west
1310            mbathy(jpi,:) = mbathy(  2  ,:)
1311         ENDIF
1312
1313         ! split last cell if possible (only where water column is 2 cell or less)
1314         DO jk = jpkm1, 1, -1
1315            zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1316            WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)
1317               mbathy(:,:) = jk
1318               bathy(:,:)  = zmax
1319            END WHERE
1320         END DO
1321 
1322         ! split top cell if possible (only where water column is 2 cell or less)
1323         DO jk = 2, jpkm1
1324            zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1325            WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy)
1326               misfdep(:,:) = jk
1327               risfdep(:,:) = zmax
1328            END WHERE
1329         END DO
1330
1331 
1332 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition
1333         DO jj = 1, jpj
1334            DO ji = 1, jpi
1335               ! find the minimum change option:
1336               ! test bathy
1337               IF (risfdep(ji,jj) .GT. 1) THEN
1338               zbathydiff =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) &
1339                 &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
1340               zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) &
1341                 &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
1342 
1343                  IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN
1344                     IF (zbathydiff .LE. zrisfdepdiff) THEN
1345                        bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )
1346                        mbathy(ji,jj)= mbathy(ji,jj) + 1
1347                     ELSE
1348                        risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )
1349                        misfdep(ji,jj) = misfdep(ji,jj) - 1
1350                     END IF
1351                  END IF
1352               END IF
1353            END DO
1354         END DO
1355 
1356          ! At least 2 levels for water thickness at T, U, and V point.
1357         DO jj = 1, jpj
1358            DO ji = 1, jpi
1359               ! find the minimum change option:
1360               ! test bathy
1361               IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
1362                  zbathydiff  =ABS(bathy(ji,jj)  - (gdepw_1d(mbathy (ji,jj)+1)&
1363                    & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
1364                  zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  )  &
1365                    & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
1366                  IF (zbathydiff .LE. zrisfdepdiff) THEN
1367                     mbathy(ji,jj) = mbathy(ji,jj) + 1
1368                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
1369                  ELSE
1370                     misfdep(ji,jj)= misfdep(ji,jj) - 1
1371                     risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )
1372                  END IF
1373               ENDIF
1374            END DO
1375         END DO
1376 
1377 ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1)
1378         DO jj = 1, jpjm1
1379            DO ji = 1, jpim1
1380               IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
1381                  zbathydiff =ABS(bathy(ji,jj    ) - (gdepw_1d(mbathy (ji,jj)+1)   &
1382                    &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat )))
1383                  zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1))  &
1384                    &  - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))
1385                  IF (zbathydiff .LE. zrisfdepdiff) THEN
1386                     mbathy(ji,jj) = mbathy(ji,jj) + 1
1387                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) &
1388                   &    + MIN( e3zps_min, e3t_1d(mbathy(ji,jj   )+1)*e3zps_rat )
1389                  ELSE
1390                     misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1
1391                     risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) &
1392                   &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )
1393                  END IF
1394               ENDIF
1395            END DO
1396         END DO
1397 
1398         IF( lk_mpp ) THEN
1399            zbathy(:,:) = FLOAT( misfdep(:,:) )
1400            CALL lbc_lnk( zbathy, 'T', 1. )
1401            misfdep(:,:) = INT( zbathy(:,:) )
1402            CALL lbc_lnk( risfdep, 'T', 1. )
1403            CALL lbc_lnk( bathy, 'T', 1. )
1404            zbathy(:,:) = FLOAT( mbathy(:,:) )
1405            CALL lbc_lnk( zbathy, 'T', 1. )
1406            mbathy(:,:) = INT( zbathy(:,:) )
1407         ENDIF
1408 ! point V misdep(ji,jj) EQ mbathy(ji,jj+1)
1409         DO jj = 1, jpjm1
1410            DO ji = 1, jpim1
1411               IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN
1412                  zbathydiff =ABS(  bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) &
1413                   &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))
1414                  zrisfdepdiff=ABS(risfdep(ji,jj  ) - (gdepw_1d(misfdep(ji,jj  )  )  &
1415                   &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat )))
1416                  IF (zbathydiff .LE. zrisfdepdiff) THEN
1417                     mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1
1418                     bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )
1419                  ELSE
1420                     misfdep(ji,jj)   = misfdep(ji,jj) - 1
1421                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat )
1422                  END IF
1423               ENDIF
1424            END DO
1425         END DO
1426 
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 
1439 ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1)
1440         DO jj = 1, jpjm1
1441            DO ji = 1, jpim1
1442               IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
1443                  zbathydiff =ABS(  bathy(ji  ,jj) - (gdepw_1d(mbathy (ji,jj)+1) &
1444                    &   + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat )))
1445                  zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) &
1446                    &  - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))
1447                  IF (zbathydiff .LE. zrisfdepdiff) THEN
1448                     mbathy(ji,jj) = mbathy(ji,jj) + 1
1449                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
1450                  ELSE
1451                     misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1
1452                     risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )
1453                  END IF
1454               ENDIF
1455            ENDDO
1456         ENDDO
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 misfdep(ji,jj) EQ bathy(ji,jj+1)
1470         DO jj = 1, jpjm1
1471            DO ji = 1, jpim1
1472               IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
1473                  zbathydiff =ABS(  bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) &
1474                      &   + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))
1475                  zrisfdepdiff=ABS(risfdep(ji  ,jj) - (gdepw_1d(misfdep(ji  ,jj)  )  &
1476                      &  - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat )))
1477                  IF (zbathydiff .LE. zrisfdepdiff) THEN
1478                     mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1
1479                     bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  )  &
1480                      &   + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat )
1481                  ELSE
1482                     misfdep(ji,jj)   = misfdep(ji  ,jj) - 1
1483                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) &
1484                      &   - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat )
1485                  END IF
1486               ENDIF
1487            ENDDO
1488         ENDDO
1489 
1490         IF( lk_mpp ) THEN
1491            zbathy(:,:) = FLOAT( misfdep(:,:) )
1492            CALL lbc_lnk( zbathy, 'T', 1. )
1493            misfdep(:,:) = INT( zbathy(:,:) )
1494            CALL lbc_lnk( risfdep, 'T', 1. )
1495            CALL lbc_lnk( bathy, 'T', 1. )
1496            zbathy(:,:) = FLOAT( mbathy(:,:) )
1497            CALL lbc_lnk( zbathy, 'T', 1. )
1498            mbathy(:,:) = INT( zbathy(:,:) )
1499         ENDIF
1500      END DO
1501      ! end dig bathy/ice shelf to be compatible
1502      ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness
1503      DO jl = 1,20
1504 
1505 ! remove single point "bay" on isf coast line in the ice shelf draft'
1506         DO jk = 2, jpk
1507            WHERE (misfdep==0) misfdep=jpk
1508            zmask=0
1509            WHERE (misfdep .LE. jk) zmask=1
1510            DO jj = 2, jpjm1
1511               DO ji = 2, jpim1
1512                  IF (misfdep(ji,jj) .EQ. jk) THEN
1513                     ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
1514                     IF (ibtest .LE. 1) THEN
1515                        risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1
1516                        IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) misfdep(ji,jj) = jpk
1517                     END IF
1518                  END IF
1519               END DO
1520            END DO
1521         END DO
1522         WHERE (misfdep==jpk)
1523             misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.
1524         END WHERE
1525         IF( lk_mpp ) THEN
1526            zbathy(:,:) = FLOAT( misfdep(:,:) )
1527            CALL lbc_lnk( zbathy, 'T', 1. )
1528            misfdep(:,:) = INT( zbathy(:,:) )
1529            CALL lbc_lnk( risfdep, 'T', 1. )
1530            CALL lbc_lnk( bathy, 'T', 1. )
1531            zbathy(:,:) = FLOAT( mbathy(:,:) )
1532            CALL lbc_lnk( zbathy, 'T', 1. )
1533            mbathy(:,:) = INT( zbathy(:,:) )
1534         ENDIF
1535 
1536 ! remove single point "bay" on bathy coast line beneath an ice shelf'
1537         DO jk = jpk,1,-1
1538            zmask=0
1539            WHERE (mbathy .GE. jk ) zmask=1
1540            DO jj = 2, jpjm1
1541               DO ji = 2, jpim1
1542                  IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) 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                        bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1
1546                        IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 0
1547                     END IF
1548                  END IF
1549               END DO
1550            END DO
1551         END DO
1552         WHERE (mbathy==0)
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 ! fill hole in ice shelf
1567         zmisfdep = misfdep
1568         zrisfdep = risfdep
1569         WHERE (zmisfdep .LE. 1) zmisfdep=jpk
1570         DO jj = 2, jpjm1
1571            DO ji = 2, jpim1
1572               ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  )
1573               ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1)
1574               IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj  ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj  ) - 1)
1575               IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj  ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj  ) - 1)
1576               IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji  ,jj-1) - 1)
1577               IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji  ,jj+1) - 1)
1578               ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
1579               IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN
1580                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp
1581               END IF
1582               IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN
1583                  misfdep(ji,jj) = ibtest
1584                  risfdep(ji,jj) = gdepw_1d(ibtest)
1585               ENDIF
1586            ENDDO
1587         ENDDO
1588 
1589         IF( lk_mpp ) THEN
1590            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
1591            CALL lbc_lnk( zbathy, 'T', 1. ) 
1592            misfdep(:,:) = INT( zbathy(:,:) ) 
1593            CALL lbc_lnk( risfdep, 'T', 1. ) 
1594            CALL lbc_lnk( bathy, 'T', 1. )
1595            zbathy(:,:) = FLOAT( mbathy(:,:) )
1596            CALL lbc_lnk( zbathy, 'T', 1. )
1597            mbathy(:,:) = INT( zbathy(:,:) )
1598         ENDIF 
1599 !
1600 !! fill hole in bathymetry
1601         zmbathy (:,:)=mbathy (:,:)
1602         DO jj = 2, jpjm1
1603            DO ji = 2, jpim1
1604               ibtestim1 = zmbathy(ji-1,jj  ) ; ibtestip1 = zmbathy(ji+1,jj  )
1605               ibtestjm1 = zmbathy(ji  ,jj-1) ; ibtestjp1 = zmbathy(ji  ,jj+1)
1606               IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj  ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj  ) + 1)
1607               IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj  ) ) ibtestip1 = 0
1608               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj-1) ) ibtestjm1 = 0
1609               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj+1) ) ibtestjp1 = 0
1610               ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
1611               IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN
1612                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ;
1613               END IF
1614               IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN
1615                  mbathy(ji,jj) = ibtest
1616                  bathy(ji,jj)  = gdepw_1d(ibtest+1) 
1617               ENDIF
1618            END DO
1619         END DO
1620         IF( lk_mpp ) THEN
1621            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
1622            CALL lbc_lnk( zbathy, 'T', 1. ) 
1623            misfdep(:,:) = INT( zbathy(:,:) ) 
1624            CALL lbc_lnk( risfdep, 'T', 1. ) 
1625            CALL lbc_lnk( bathy, 'T', 1. )
1626            zbathy(:,:) = FLOAT( mbathy(:,:) )
1627            CALL lbc_lnk( zbathy, 'T', 1. )
1628            mbathy(:,:) = INT( zbathy(:,:) )
1629         ENDIF 
1630 ! if not compatible after all check (ie U point water column less than 2 cells), mask U
1631         DO jj = 1, jpjm1
1632            DO ji = 1, jpim1
1633               IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN
1634                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ;
1635               END IF
1636            END DO
1637         END DO
1638         IF( lk_mpp ) THEN
1639            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
1640            CALL lbc_lnk( zbathy, 'T', 1. ) 
1641            misfdep(:,:) = INT( zbathy(:,:) ) 
1642            CALL lbc_lnk( risfdep, 'T', 1. ) 
1643            CALL lbc_lnk( bathy, 'T', 1. )
1644            zbathy(:,:) = FLOAT( mbathy(:,:) )
1645            CALL lbc_lnk( zbathy, 'T', 1. )
1646            mbathy(:,:) = INT( zbathy(:,:) )
1647         ENDIF 
1648 ! if not compatible after all check (ie U point water column less than 2 cells), mask U
1649         DO jj = 1, jpjm1
1650            DO ji = 1, jpim1
1651               IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN
1652                  mbathy(ji+1,jj)  = mbathy(ji+1,jj) - 1;   bathy(ji+1,jj)   = gdepw_1d(mbathy(ji+1,jj)+1) ;
1653               END IF
1654            END DO
1655         END DO
1656         IF( lk_mpp ) THEN
1657            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
1658            CALL lbc_lnk( zbathy, 'T', 1. ) 
1659            misfdep(:,:) = INT( zbathy(:,:) ) 
1660            CALL lbc_lnk( risfdep, 'T', 1. ) 
1661            CALL lbc_lnk( bathy, 'T', 1. )
1662            zbathy(:,:) = FLOAT( mbathy(:,:) )
1663            CALL lbc_lnk( zbathy, 'T', 1. )
1664            mbathy(:,:) = INT( zbathy(:,:) )
1665         ENDIF 
1666 ! if not compatible after all check (ie V point water column less than 2 cells), mask V
1667         DO jj = 1, jpjm1
1668            DO ji = 1, jpi
1669               IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN
1670                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ;
1671               END IF
1672            END DO
1673         END DO
1674         IF( lk_mpp ) THEN
1675            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
1676            CALL lbc_lnk( zbathy, 'T', 1. ) 
1677            misfdep(:,:) = INT( zbathy(:,:) ) 
1678            CALL lbc_lnk( risfdep, 'T', 1. ) 
1679            CALL lbc_lnk( bathy, 'T', 1. )
1680            zbathy(:,:) = FLOAT( mbathy(:,:) )
1681            CALL lbc_lnk( zbathy, 'T', 1. )
1682            mbathy(:,:) = INT( zbathy(:,:) )
1683         ENDIF 
1684 ! if not compatible after all check (ie V point water column less than 2 cells), mask V
1685         DO jj = 1, jpjm1
1686            DO ji = 1, jpi
1687               IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN
1688                  mbathy(ji,jj+1)  = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1)   = gdepw_1d(mbathy(ji,jj+1)+1) ;
1689               END IF
1690            END DO
1691         END DO
1692         IF( lk_mpp ) THEN
1693            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
1694            CALL lbc_lnk( zbathy, 'T', 1. ) 
1695            misfdep(:,:) = INT( zbathy(:,:) ) 
1696            CALL lbc_lnk( risfdep, 'T', 1. ) 
1697            CALL lbc_lnk( bathy, 'T', 1. )
1698            zbathy(:,:) = FLOAT( mbathy(:,:) )
1699            CALL lbc_lnk( zbathy, 'T', 1. )
1700            mbathy(:,:) = INT( zbathy(:,:) )
1701         ENDIF 
1702 ! if not compatible after all check, mask T
1703         DO jj = 1, jpj
1704            DO ji = 1, jpi
1705               IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN
1706                  misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj)  = 0 ; bathy(ji,jj)   = 0._wp ;
1707               END IF
1708            END DO
1709         END DO
1710 
1711         WHERE (mbathy(:,:) == 1)
1712            mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp
1713         END WHERE
1714      END DO 
1715! end check compatibility ice shelf/bathy
1716      ! remove very shallow ice shelf (less than ~ 10m if 75L)
1717      WHERE (misfdep(:,:) <= 5)
1718         misfdep = 1; risfdep = 0.0_wp;
1719      END WHERE
1720
1721      IF( icompt == 0 ) THEN
1722         IF(lwp) WRITE(numout,*)'     no points with ice shelf too close to bathymetry' 
1723      ELSE
1724         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry' 
1725      ENDIF 
1726
1727      ! compute scale factor and depth at T- and W- points
1728      DO jj = 1, jpj
1729         DO ji = 1, jpi
1730            ik = mbathy(ji,jj)
1731            IF( ik > 0 ) THEN               ! ocean point only
1732               ! max ocean level case
1733               IF( ik == jpkm1 ) THEN
1734                  zdepwp = bathy(ji,jj)
1735                  ze3tp  = bathy(ji,jj) - gdepw_1d(ik)
1736                  ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
1737                  e3t_0(ji,jj,ik  ) = ze3tp
1738                  e3t_0(ji,jj,ik+1) = ze3tp
1739                  e3w_0(ji,jj,ik  ) = ze3wp
1740                  e3w_0(ji,jj,ik+1) = ze3tp
1741                  gdepw_0(ji,jj,ik+1) = zdepwp
1742                  gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp
1743                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
1744                  !
1745               ELSE                         ! standard case
1746                  IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
1747                  ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
1748                  ENDIF
1749      !            gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
1750!gm Bug?  check the gdepw_1d
1751                  !       ... on ik
1752                  gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   &
1753                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   &
1754                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) ))
1755                  e3t_0  (ji,jj,ik  ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik  )
1756                  e3w_0  (ji,jj,ik  ) = gdept_0(ji,jj,ik  ) - gdept_1d(ik-1)
1757                  !       ... on ik+1
1758                  e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
1759                  e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
1760               ENDIF
1761            ENDIF
1762         END DO
1763      END DO
1764      !
1765      it = 0
1766      DO jj = 1, jpj
1767         DO ji = 1, jpi
1768            ik = mbathy(ji,jj)
1769            IF( ik > 0 ) THEN               ! ocean point only
1770               e3tp (ji,jj) = e3t_0(ji,jj,ik)
1771               e3wp (ji,jj) = e3w_0(ji,jj,ik)
1772               ! test
1773               zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  )
1774               IF( zdiff <= 0._wp .AND. lwp ) THEN
1775                  it = it + 1
1776                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj
1777                  WRITE(numout,*) ' bathy = ', bathy(ji,jj)
1778                  WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
1779                  WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  )
1780               ENDIF
1781            ENDIF
1782         END DO
1783      END DO
1784      !
1785      ! (ISF) Definition of e3t, u, v, w for ISF case
1786      DO jj = 1, jpj 
1787         DO ji = 1, jpi 
1788            ik = misfdep(ji,jj) 
1789            IF( ik > 1 ) THEN               ! ice shelf point only
1790               IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik) 
1791               gdepw_0(ji,jj,ik) = risfdep(ji,jj) 
1792!gm Bug?  check the gdepw_0
1793            !       ... on ik
1794               gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   & 
1795                  &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   & 
1796                  &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      ) 
1797               e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 
1798               e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik)
1799
1800               IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)
1801                  e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 
1802               ENDIF 
1803            !       ... on ik / ik-1
1804               e3w_0  (ji,jj,ik  ) = e3t_0  (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))
1805               e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
1806! The next line isn't required and doesn't affect results - included for consistency with bathymetry code
1807               gdept_0(ji,jj,ik-1) = gdept_1d(ik-1)
1808            ENDIF
1809         END DO
1810      END DO
1811   
1812      it = 0 
1813      DO jj = 1, jpj 
1814         DO ji = 1, jpi 
1815            ik = misfdep(ji,jj) 
1816            IF( ik > 1 ) THEN               ! ice shelf point only
1817               e3tp (ji,jj) = e3t_0(ji,jj,ik  ) 
1818               e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 
1819            ! test
1820               zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  ) 
1821               IF( zdiff <= 0. .AND. lwp ) THEN 
1822                  it = it + 1 
1823                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
1824                  WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 
1825                  WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
1826                  WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj) 
1827               ENDIF
1828            ENDIF
1829         END DO
1830      END DO
1831
1832      CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep )
1833      CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy )
1834
1835      IF( nn_timing == 1 )  CALL timing_stop('zgr_isf')
1836     
1837   END SUBROUTINE
1838
1839   SUBROUTINE zgr_sco
1840      !!----------------------------------------------------------------------
1841      !!                  ***  ROUTINE zgr_sco  ***
1842      !!                     
1843      !! ** Purpose :   define the s-coordinate system
1844      !!
1845      !! ** Method  :   s-coordinate
1846      !!         The depth of model levels is defined as the product of an
1847      !!      analytical function by the local bathymetry, while the vertical
1848      !!      scale factors are defined as the product of the first derivative
1849      !!      of the analytical function by the bathymetry.
1850      !!      (this solution save memory as depth and scale factors are not
1851      !!      3d fields)
1852      !!          - Read bathymetry (in meters) at t-point and compute the
1853      !!         bathymetry at u-, v-, and f-points.
1854      !!            hbatu = mi( hbatt )
1855      !!            hbatv = mj( hbatt )
1856      !!            hbatf = mi( mj( hbatt ) )
1857      !!          - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
1858      !!         function and its derivative given as function.
1859      !!            z_gsigt(k) = fssig (k    )
1860      !!            z_gsigw(k) = fssig (k-0.5)
1861      !!            z_esigt(k) = fsdsig(k    )
1862      !!            z_esigw(k) = fsdsig(k-0.5)
1863      !!      Three options for stretching are give, and they can be modified
1864      !!      following the users requirements. Nevertheless, the output as
1865      !!      well as the way to compute the model levels and scale factors
1866      !!      must be respected in order to insure second order accuracy
1867      !!      schemes.
1868      !!
1869      !!      The three methods for stretching available are:
1870      !!
1871      !!           s_sh94 (Song and Haidvogel 1994)
1872      !!                a sinh/tanh function that allows sigma and stretched sigma
1873      !!
1874      !!           s_sf12 (Siddorn and Furner 2012?)
1875      !!                allows the maintenance of fixed surface and or
1876      !!                bottom cell resolutions (cf. geopotential coordinates)
1877      !!                within an analytically derived stretched S-coordinate framework.
1878      !!
1879      !!          s_tanh  (Madec et al 1996)
1880      !!                a cosh/tanh function that gives stretched coordinates       
1881      !!
1882      !!----------------------------------------------------------------------
1883      !
1884      INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
1885      INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
1886      INTEGER  ::   ios                      ! Local integer output status for namelist read
1887      REAL(wp) ::   zrmax, ztaper   ! temporary scalars
1888      REAL(wp) ::   zrfact
1889      !
1890      REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
1891      REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
1892
1893      NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
1894                           rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
1895     !!----------------------------------------------------------------------
1896      !
1897      IF( nn_timing == 1 )  CALL timing_start('zgr_sco')
1898      !
1899      CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
1900      !
1901      REWIND( numnam_ref )              ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
1902      READ  ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901)
1903901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist', lwp )
1904
1905      REWIND( numnam_cfg )              ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters
1906      READ  ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 )
1907902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist', lwp )
1908      IF(lwm) WRITE ( numond, namzgr_sco )
1909
1910      IF(lwp) THEN                           ! control print
1911         WRITE(numout,*)
1912         WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
1913         WRITE(numout,*) '~~~~~~~~~~~'
1914         WRITE(numout,*) '   Namelist namzgr_sco'
1915         WRITE(numout,*) '     stretching coeffs '
1916         WRITE(numout,*) '        maximum depth of s-bottom surface (>0)       rn_sbot_max   = ',rn_sbot_max
1917         WRITE(numout,*) '        minimum depth of s-bottom surface (>0)       rn_sbot_min   = ',rn_sbot_min
1918         WRITE(numout,*) '        Critical depth                               rn_hc         = ',rn_hc
1919         WRITE(numout,*) '        maximum cut-off r-value allowed              rn_rmax       = ',rn_rmax
1920         WRITE(numout,*) '     Song and Haidvogel 1994 stretching              ln_s_sh94     = ',ln_s_sh94
1921         WRITE(numout,*) '        Song and Haidvogel 1994 stretching coefficients'
1922         WRITE(numout,*) '        surface control parameter (0<=rn_theta<=20)  rn_theta      = ',rn_theta
1923         WRITE(numout,*) '        bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ',rn_thetb
1924         WRITE(numout,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb
1925         WRITE(numout,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12
1926         WRITE(numout,*) '        switching to sigma (T) or Z (F) at H<Hc      ln_sigcrit    = ',ln_sigcrit
1927         WRITE(numout,*) '        Siddorn and Furner 2012 stretching coefficients'
1928         WRITE(numout,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha
1929         WRITE(numout,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold
1930         WRITE(numout,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs
1931         WRITE(numout,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a
1932         WRITE(numout,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b
1933         WRITE(numout,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b'
1934      ENDIF
1935
1936      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate
1937      hifu(:,:) = rn_sbot_min
1938      hifv(:,:) = rn_sbot_min
1939      hiff(:,:) = rn_sbot_min
1940
1941      !                                        ! set maximum ocean depth
1942      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1943
1944      DO jj = 1, jpj
1945         DO ji = 1, jpi
1946           IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1947         END DO
1948      END DO
1949      !                                        ! =============================
1950      !                                        ! Define the envelop bathymetry   (hbatt)
1951      !                                        ! =============================
1952      ! use r-value to create hybrid coordinates
1953      zenv(:,:) = bathy(:,:)
1954      !
1955      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
1956      DO jj = 1, jpj
1957         DO ji = 1, jpi
1958           IF( bathy(ji,jj) == 0._wp ) THEN
1959             iip1 = MIN( ji+1, jpi )
1960             ijp1 = MIN( jj+1, jpj )
1961             iim1 = MAX( ji-1, 1 )
1962             ijm1 = MAX( jj-1, 1 )
1963             IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) +              &
1964        &         bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
1965               zenv(ji,jj) = rn_sbot_min
1966             ENDIF
1967           ENDIF
1968         END DO
1969      END DO
1970      ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1971      CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1972      !
1973      ! smooth the bathymetry (if required)
1974      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
1975      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
1976      !
1977      jl = 0
1978      zrmax = 1._wp
1979      !   
1980      !     
1981      ! set scaling factor used in reducing vertical gradients
1982      zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
1983      !
1984      ! initialise temporary evelope depth arrays
1985      ztmpi1(:,:) = zenv(:,:)
1986      ztmpi2(:,:) = zenv(:,:)
1987      ztmpj1(:,:) = zenv(:,:)
1988      ztmpj2(:,:) = zenv(:,:)
1989      !
1990      ! initialise temporary r-value arrays
1991      zri(:,:) = 1._wp
1992      zrj(:,:) = 1._wp
1993      !                                                            ! ================ !
1994      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) !  Iterative loop  !
1995         !                                                         ! ================ !
1996         jl = jl + 1
1997         zrmax = 0._wp
1998         ! we set zrmax from previous r-values (zri and zrj) first
1999         ! if set after current r-value calculation (as previously)
2000         ! we could exit DO WHILE prematurely before checking r-value
2001         ! of current zenv
2002         DO jj = 1, nlcj
2003            DO ji = 1, nlci
2004               zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
2005            END DO
2006         END DO
2007         zri(:,:) = 0._wp
2008         zrj(:,:) = 0._wp
2009         DO jj = 1, nlcj
2010            DO ji = 1, nlci
2011               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
2012               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
2013               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
2014                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
2015               END IF
2016               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
2017                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
2018               END IF
2019               IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact
2020               IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact
2021               IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact
2022               IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact
2023            END DO
2024         END DO
2025         IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
2026         !
2027         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax
2028         !
2029         DO jj = 1, nlcj
2030            DO ji = 1, nlci
2031               zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
2032            END DO
2033         END DO
2034         ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
2035         CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
2036         !                                                  ! ================ !
2037      END DO                                                !     End loop     !
2038      !                                                     ! ================ !
2039      DO jj = 1, jpj
2040         DO ji = 1, jpi
2041            zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
2042         END DO
2043      END DO
2044      !
2045      ! Envelope bathymetry saved in hbatt
2046      hbatt(:,:) = zenv(:,:) 
2047      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
2048         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
2049         DO jj = 1, jpj
2050            DO ji = 1, jpi
2051               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
2052               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
2053            END DO
2054         END DO
2055      ENDIF
2056      !
2057      IF(lwp) THEN                             ! Control print
2058         WRITE(numout,*)
2059         WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
2060         WRITE(numout,*)
2061         CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )
2062         IF( nprint == 1 )   THEN       
2063            WRITE(numout,*) ' bathy  MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
2064            WRITE(numout,*) ' hbatt  MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
2065         ENDIF
2066      ENDIF
2067
2068      !                                        ! ==============================
2069      !                                        !   hbatu, hbatv, hbatf fields
2070      !                                        ! ==============================
2071      IF(lwp) THEN
2072         WRITE(numout,*)
2073         WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
2074      ENDIF
2075      hbatu(:,:) = rn_sbot_min
2076      hbatv(:,:) = rn_sbot_min
2077      hbatf(:,:) = rn_sbot_min
2078      DO jj = 1, jpjm1
2079        DO ji = 1, jpim1   ! NO vector opt.
2080           hbatu(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
2081           hbatv(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
2082           hbatf(ji,jj) = 0.25_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
2083              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
2084        END DO
2085      END DO
2086      !
2087      ! Apply lateral boundary condition
2088!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
2089      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk( hbatu, 'U', 1._wp )
2090      DO jj = 1, jpj
2091         DO ji = 1, jpi
2092            IF( hbatu(ji,jj) == 0._wp ) THEN
2093               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
2094               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
2095            ENDIF
2096         END DO
2097      END DO
2098      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk( hbatv, 'V', 1._wp )
2099      DO jj = 1, jpj
2100         DO ji = 1, jpi
2101            IF( hbatv(ji,jj) == 0._wp ) THEN
2102               IF( zhbat(ji,jj) == 0._wp )   hbatv(ji,jj) = rn_sbot_min
2103               IF( zhbat(ji,jj) /= 0._wp )   hbatv(ji,jj) = zhbat(ji,jj)
2104            ENDIF
2105         END DO
2106      END DO
2107      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk( hbatf, 'F', 1._wp )
2108      DO jj = 1, jpj
2109         DO ji = 1, jpi
2110            IF( hbatf(ji,jj) == 0._wp ) THEN
2111               IF( zhbat(ji,jj) == 0._wp )   hbatf(ji,jj) = rn_sbot_min
2112               IF( zhbat(ji,jj) /= 0._wp )   hbatf(ji,jj) = zhbat(ji,jj)
2113            ENDIF
2114         END DO
2115      END DO
2116
2117!!bug:  key_helsinki a verifer
2118      hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
2119      hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
2120      hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
2121      hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
2122
2123      IF( nprint == 1 .AND. lwp )   THEN
2124         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
2125            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
2126         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
2127            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
2128         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
2129            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
2130         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
2131            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
2132      ENDIF
2133!! helsinki
2134
2135      !                                            ! =======================
2136      !                                            !   s-ccordinate fields     (gdep., e3.)
2137      !                                            ! =======================
2138      !
2139      ! non-dimensional "sigma" for model level depth at w- and t-levels
2140
2141
2142!========================================================================
2143! Song and Haidvogel  1994 (ln_s_sh94=T)
2144! Siddorn and Furner 2012 (ln_sf12=T)
2145! or  tanh function       (both false)                   
2146!========================================================================
2147      IF      ( ln_s_sh94 ) THEN
2148                           CALL s_sh94()
2149      ELSE IF ( ln_s_sf12 ) THEN
2150                           CALL s_sf12()
2151      ELSE                 
2152                           CALL s_tanh()
2153      ENDIF
2154
2155      CALL lbc_lnk( e3t_0 , 'T', 1._wp )
2156      CALL lbc_lnk( e3u_0 , 'U', 1._wp )
2157      CALL lbc_lnk( e3v_0 , 'V', 1._wp )
2158      CALL lbc_lnk( e3f_0 , 'F', 1._wp )
2159      CALL lbc_lnk( e3w_0 , 'W', 1._wp )
2160      CALL lbc_lnk( e3uw_0, 'U', 1._wp )
2161      CALL lbc_lnk( e3vw_0, 'V', 1._wp )
2162
2163      fsdepw(:,:,:) = gdepw_0 (:,:,:)
2164      fsde3w(:,:,:) = gdep3w_0(:,:,:)
2165      !
2166      where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0
2167      where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0
2168      where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0
2169      where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0
2170      where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0
2171      where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0
2172      where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0
2173
2174#if defined key_agrif
2175      ! Ensure meaningful vertical scale factors in ghost lines/columns
2176      IF( .NOT. Agrif_Root() ) THEN
2177         
2178         IF((nbondi == -1).OR.(nbondi == 2)) THEN
2179            e3u_0(1,:,:) = e3u_0(2,:,:)
2180         ENDIF
2181         !
2182         IF((nbondi ==  1).OR.(nbondi == 2)) THEN
2183            e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:)
2184         ENDIF
2185         !
2186         IF((nbondj == -1).OR.(nbondj == 2)) THEN
2187            e3v_0(:,1,:) = e3v_0(:,2,:)
2188         ENDIF
2189         !
2190         IF((nbondj ==  1).OR.(nbondj == 2)) THEN
2191            e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:)
2192         ENDIF
2193         !
2194      ENDIF
2195#endif
2196
2197      fsdept(:,:,:) = gdept_0 (:,:,:)
2198      fsdepw(:,:,:) = gdepw_0 (:,:,:)
2199      fsde3w(:,:,:) = gdep3w_0(:,:,:)
2200      fse3t (:,:,:) = e3t_0   (:,:,:)
2201      fse3u (:,:,:) = e3u_0   (:,:,:)
2202      fse3v (:,:,:) = e3v_0   (:,:,:)
2203      fse3f (:,:,:) = e3f_0   (:,:,:)
2204      fse3w (:,:,:) = e3w_0   (:,:,:)
2205      fse3uw(:,:,:) = e3uw_0  (:,:,:)
2206      fse3vw(:,:,:) = e3vw_0  (:,:,:)
2207!!
2208      ! HYBRID :
2209      DO jj = 1, jpj
2210         DO ji = 1, jpi
2211            DO jk = 1, jpkm1
2212               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
2213            END DO
2214            IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0
2215         END DO
2216      END DO
2217      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
2218         &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
2219
2220      IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain
2221         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) )
2222         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
2223            &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w '  , MINVAL( gdep3w_0(:,:,:) )
2224         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0   (:,:,:) ),   &
2225            &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0   (:,:,:) ),   &
2226            &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0  (:,:,:) ),   &
2227            &                          ' w ', MINVAL( e3w_0  (:,:,:) )
2228
2229         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
2230            &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w '  , MAXVAL( gdep3w_0(:,:,:) )
2231         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0   (:,:,:) ),   &
2232            &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0   (:,:,:) ),   &
2233            &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0  (:,:,:) ),   &
2234            &                          ' w ', MAXVAL( e3w_0  (:,:,:) )
2235      ENDIF
2236      !  END DO
2237      IF(lwp) THEN                                  ! selected vertical profiles
2238         WRITE(numout,*)
2239         WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
2240         WRITE(numout,*) ' ~~~~~~  --------------------'
2241         WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2242         WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     &
2243            &                                 e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
2244         DO jj = mj0(20), mj1(20)
2245            DO ji = mi0(20), mi1(20)
2246               WRITE(numout,*)
2247               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
2248               WRITE(numout,*) ' ~~~~~~  --------------------'
2249               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2250               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
2251                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2252            END DO
2253         END DO
2254         DO jj = mj0(74), mj1(74)
2255            DO ji = mi0(100), mi1(100)
2256               WRITE(numout,*)
2257               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
2258               WRITE(numout,*) ' ~~~~~~  --------------------'
2259               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2260               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
2261                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2262            END DO
2263         END DO
2264      ENDIF
2265
2266!================================================================================
2267! check the coordinate makes sense
2268!================================================================================
2269      DO ji = 1, jpi
2270         DO jj = 1, jpj
2271
2272            IF( hbatt(ji,jj) > 0._wp) THEN
2273               DO jk = 1, mbathy(ji,jj)
2274                 ! check coordinate is monotonically increasing
2275                 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
2276                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
2277                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
2278                    WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
2279                    WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
2280                    CALL ctl_stop( ctmp1 )
2281                 ENDIF
2282                 ! and check it has never gone negative
2283                 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
2284                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
2285                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk
2286                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2287                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2288                    CALL ctl_stop( ctmp1 )
2289                 ENDIF
2290                 ! and check it never exceeds the total depth
2291                 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
2292                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
2293                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
2294                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2295                    CALL ctl_stop( ctmp1 )
2296                 ENDIF
2297               END DO
2298
2299               DO jk = 1, mbathy(ji,jj)-1
2300                 ! and check it never exceeds the total depth
2301                IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
2302                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
2303                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
2304                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2305                    CALL ctl_stop( ctmp1 )
2306                 ENDIF
2307               END DO
2308
2309            ENDIF
2310
2311         END DO
2312      END DO
2313      !
2314      CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
2315      !
2316      IF( nn_timing == 1 )  CALL timing_stop('zgr_sco')
2317      !
2318   END SUBROUTINE zgr_sco
2319
2320!!======================================================================
2321   SUBROUTINE s_sh94()
2322
2323      !!----------------------------------------------------------------------
2324      !!                  ***  ROUTINE s_sh94  ***
2325      !!                     
2326      !! ** Purpose :   stretch the s-coordinate system
2327      !!
2328      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994
2329      !!                mixed S/sigma coordinate
2330      !!
2331      !! Reference : Song and Haidvogel 1994.
2332      !!----------------------------------------------------------------------
2333      !
2334      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2335      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
2336      !
2337      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2338      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
2339
2340      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2341      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2342
2343      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
2344      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
2345      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
2346      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
2347
2348      DO ji = 1, jpi
2349         DO jj = 1, jpj
2350
2351            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
2352               DO jk = 1, jpk
2353                  z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
2354                  z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
2355               END DO
2356            ELSE ! shallow water, uniform sigma
2357               DO jk = 1, jpk
2358                  z_gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
2359                  z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
2360                  END DO
2361            ENDIF
2362            !
2363            DO jk = 1, jpkm1
2364               z_esigt3(ji,jj,jk  ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2365               z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2366            END DO
2367            z_esigw3(ji,jj,1  ) = 2._wp * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) )
2368            z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
2369            !
2370            ! Coefficients for vertical depth as the sum of e3w scale factors
2371            z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
2372            DO jk = 2, jpk
2373               z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2374            END DO
2375            !
2376            DO jk = 1, jpk
2377               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2378               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2379               gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
2380               gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
2381               gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
2382            END DO
2383           !
2384         END DO   ! for all jj's
2385      END DO    ! for all ji's
2386
2387      DO ji = 1, jpim1
2388         DO jj = 1, jpjm1
2389            DO jk = 1, jpk
2390               z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
2391                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2392               z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
2393                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2394               z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     &
2395                  &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   &
2396                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2397               z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
2398                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2399               z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
2400                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2401               !
2402               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2403               e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2404               e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2405               e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2406               !
2407               e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2408               e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2409               e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2410            END DO
2411        END DO
2412      END DO
2413
2414      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2415      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2416
2417   END SUBROUTINE s_sh94
2418
2419   SUBROUTINE s_sf12
2420
2421      !!----------------------------------------------------------------------
2422      !!                  ***  ROUTINE s_sf12 ***
2423      !!                     
2424      !! ** Purpose :   stretch the s-coordinate system
2425      !!
2426      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
2427      !!                mixed S/sigma/Z coordinate
2428      !!
2429      !!                This method allows the maintenance of fixed surface and or
2430      !!                bottom cell resolutions (cf. geopotential coordinates)
2431      !!                within an analytically derived stretched S-coordinate framework.
2432      !!
2433      !!
2434      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
2435      !!----------------------------------------------------------------------
2436      !
2437      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2438      REAL(wp) ::   zsmth               ! smoothing around critical depth
2439      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
2440      !
2441      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2442      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
2443
2444      !
2445      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2446      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2447
2448      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
2449      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
2450      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
2451      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
2452
2453      DO ji = 1, jpi
2454         DO jj = 1, jpj
2455
2456          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
2457             
2458              zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
2459                                                     ! could be changed by users but care must be taken to do so carefully
2460              zzb = 1.0_wp-(zzb/hbatt(ji,jj))
2461           
2462              zzs = rn_zs / hbatt(ji,jj) 
2463             
2464              IF (rn_efold /= 0.0_wp) THEN
2465                zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
2466              ELSE
2467                zsmth = 1.0_wp 
2468              ENDIF
2469               
2470              DO jk = 1, jpk
2471                z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
2472                z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
2473              ENDDO
2474              z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth  )
2475              z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth  )
2476 
2477          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
2478
2479            DO jk = 1, jpk
2480              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
2481              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
2482            END DO
2483
2484          ELSE  ! shallow water, z coordinates
2485
2486            DO jk = 1, jpk
2487              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2488              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2489            END DO
2490
2491          ENDIF
2492
2493          DO jk = 1, jpkm1
2494             z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2495             z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2496          END DO
2497          z_esigw3(ji,jj,1  ) = 2.0_wp * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ))
2498          z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
2499
2500          ! Coefficients for vertical depth as the sum of e3w scale factors
2501          z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
2502          DO jk = 2, jpk
2503             z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2504          END DO
2505
2506          DO jk = 1, jpk
2507             gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
2508             gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
2509             gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
2510          END DO
2511
2512        ENDDO   ! for all jj's
2513      ENDDO    ! for all ji's
2514
2515      DO ji=1,jpi-1
2516        DO jj=1,jpj-1
2517
2518          DO jk = 1, jpk
2519                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / &
2520                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2521                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / &
2522                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2523                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  &
2524                                      hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / &
2525                                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2526                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / &
2527                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2528                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / &
2529                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2530
2531             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
2532             e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
2533             e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
2534             e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
2535             !
2536             e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
2537             e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
2538             e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
2539          END DO
2540
2541        ENDDO
2542      ENDDO
2543      !
2544      CALL lbc_lnk(e3t_0 ,'T',1.) ; CALL lbc_lnk(e3u_0 ,'T',1.)
2545      CALL lbc_lnk(e3v_0 ,'T',1.) ; CALL lbc_lnk(e3f_0 ,'T',1.)
2546      CALL lbc_lnk(e3w_0 ,'T',1.)
2547      CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.)
2548      !
2549      !                                               ! =============
2550
2551      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2552      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2553
2554   END SUBROUTINE s_sf12
2555
2556   SUBROUTINE s_tanh()
2557
2558      !!----------------------------------------------------------------------
2559      !!                  ***  ROUTINE s_tanh***
2560      !!                     
2561      !! ** Purpose :   stretch the s-coordinate system
2562      !!
2563      !! ** Method  :   s-coordinate stretch
2564      !!
2565      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
2566      !!----------------------------------------------------------------------
2567
2568      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2569      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
2570
2571      REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
2572      REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw
2573
2574      CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      )
2575      CALL wrk_alloc( jpk, z_esigt, z_esigw                                               )
2576
2577      z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp   ;   z_gsi3w  = 0._wp
2578      z_esigt  = 0._wp   ;   z_esigw  = 0._wp 
2579
2580      DO jk = 1, jpk
2581        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
2582        z_gsigt(jk) = -fssig( REAL(jk,wp)        )
2583      END DO
2584      IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk)
2585      !
2586      ! Coefficients for vertical scale factors at w-, t- levels
2587!!gm bug :  define it from analytical function, not like juste bellow....
2588!!gm        or betteroffer the 2 possibilities....
2589      DO jk = 1, jpkm1
2590         z_esigt(jk  ) = z_gsigw(jk+1) - z_gsigw(jk)
2591         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
2592      END DO
2593      z_esigw( 1 ) = 2._wp * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
2594      z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
2595      !
2596      ! Coefficients for vertical depth as the sum of e3w scale factors
2597      z_gsi3w(1) = 0.5_wp * z_esigw(1)
2598      DO jk = 2, jpk
2599         z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
2600      END DO
2601!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
2602      DO jk = 1, jpk
2603         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2604         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2605         gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
2606         gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
2607         gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )
2608      END DO
2609!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
2610      DO jj = 1, jpj
2611         DO ji = 1, jpi
2612            DO jk = 1, jpk
2613              e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2614              e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2615              e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2616              e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
2617              !
2618              e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2619              e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2620              e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2621            END DO
2622         END DO
2623      END DO
2624
2625      CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      )
2626      CALL wrk_dealloc( jpk, z_esigt, z_esigw                                               )
2627
2628   END SUBROUTINE s_tanh
2629
2630   FUNCTION fssig( pk ) RESULT( pf )
2631      !!----------------------------------------------------------------------
2632      !!                 ***  ROUTINE fssig ***
2633      !!       
2634      !! ** Purpose :   provide the analytical function in s-coordinate
2635      !!         
2636      !! ** Method  :   the function provide the non-dimensional position of
2637      !!                T and W (i.e. between 0 and 1)
2638      !!                T-points at integer values (between 1 and jpk)
2639      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2640      !!----------------------------------------------------------------------
2641      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
2642      REAL(wp)             ::   pf   ! sigma value
2643      !!----------------------------------------------------------------------
2644      !
2645      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   &
2646         &     - TANH( rn_thetb * rn_theta                                )  )   &
2647         & * (   COSH( rn_theta                           )                      &
2648         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              &
2649         & / ( 2._wp * SINH( rn_theta ) )
2650      !
2651   END FUNCTION fssig
2652
2653
2654   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
2655      !!----------------------------------------------------------------------
2656      !!                 ***  ROUTINE fssig1 ***
2657      !!
2658      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
2659      !!
2660      !! ** Method  :   the function provides the non-dimensional position of
2661      !!                T and W (i.e. between 0 and 1)
2662      !!                T-points at integer values (between 1 and jpk)
2663      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2664      !!----------------------------------------------------------------------
2665      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
2666      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
2667      REAL(wp)             ::   pf1   ! sigma value
2668      !!----------------------------------------------------------------------
2669      !
2670      IF ( rn_theta == 0 ) then      ! uniform sigma
2671         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
2672      ELSE                        ! stretched sigma
2673         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
2674            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
2675            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
2676      ENDIF
2677      !
2678   END FUNCTION fssig1
2679
2680
2681   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
2682      !!----------------------------------------------------------------------
2683      !!                 ***  ROUTINE fgamma  ***
2684      !!
2685      !! ** Purpose :   provide analytical function for the s-coordinate
2686      !!
2687      !! ** Method  :   the function provides the non-dimensional position of
2688      !!                T and W (i.e. between 0 and 1)
2689      !!                T-points at integer values (between 1 and jpk)
2690      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2691      !!
2692      !!                This method allows the maintenance of fixed surface and or
2693      !!                bottom cell resolutions (cf. geopotential coordinates)
2694      !!                within an analytically derived stretched S-coordinate framework.
2695      !!
2696      !! Reference  :   Siddorn and Furner, in prep
2697      !!----------------------------------------------------------------------
2698      REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate
2699      REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate
2700      REAL(wp), INTENT(in   ) ::   pzb           ! Bottom box depth
2701      REAL(wp), INTENT(in   ) ::   pzs           ! surface box depth
2702      REAL(wp), INTENT(in   ) ::   psmth       ! Smoothing parameter
2703      REAL(wp)                ::   za1,za2,za3    ! local variables
2704      REAL(wp)                ::   zn1,zn2        ! local variables
2705      REAL(wp)                ::   za,zb,zx       ! local variables
2706      integer                 ::   jk
2707      !!----------------------------------------------------------------------
2708      !
2709
2710      zn1  =  1./(jpk-1.)
2711      zn2  =  1. -  zn1
2712
2713      za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 
2714      za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp)
2715      za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
2716     
2717      za = pzb - za3*(pzs-za1)-za2
2718      za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) )
2719      zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
2720      zx = 1.0_wp-za/2.0_wp-zb
2721 
2722      DO jk = 1, jpk
2723        p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp +  &
2724                    & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &
2725                    &      (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
2726        p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth)
2727      ENDDO 
2728
2729      !
2730   END FUNCTION fgamma
2731
2732   !!======================================================================
2733END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.