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

Last change on this file since 4754 was 4754, checked in by rfurner, 7 years ago

change to give a one level model in surge case

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