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

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

source: branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5066

Last change on this file since 5066 was 5066, checked in by rfurner, 9 years ago

added current state of wetting and drying code to test...note it does not work

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