source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 12 months ago

The Dr Hook changes from my perl code.

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