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

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

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 4292

Last change on this file since 4292 was 4292, checked in by cetlod, 10 years ago

dev_MERGE_2013 : 1st step of the merge, see ticket #1185

  • Property svn:keywords set to Id
File size: 98.5 KB
Line 
1MODULE domzgr
2   !!==============================================================================
3   !!                       ***  MODULE domzgr   ***
4   !! Ocean initialization : domain initialization
5   !!==============================================================================
6   !! History :  OPA  ! 1995-12  (G. Madec)  Original code : s vertical coordinate
7   !!                 ! 1997-07  (G. Madec)  lbc_lnk call
8   !!                 ! 1997-04  (J.-O. Beismann)
9   !!            8.5  ! 2002-09  (A. Bozec, G. Madec)  F90: Free form and module
10   !!             -   ! 2002-09  (A. de Miranda)  rigid-lid + islands
11   !!  NEMO      1.0  ! 2003-08  (G. Madec)  F90: Free form and module
12   !!             -   ! 2005-10  (A. Beckmann)  modifications for hybrid s-ccordinates & new stretching function
13   !!            2.0  ! 2006-04  (R. Benshila, G. Madec)  add zgr_zco
14   !!            3.0  ! 2008-06  (G. Madec)  insertion of domzgr_zps.h90 & conding style
15   !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option
16   !!            3.3  ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level
17   !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn and Furner stretching function
18   !!            3.4  ! 2012-12  (R. Bourdalle-Badie and G. Reffray)  modify C1D case 
19   !!----------------------------------------------------------------------
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      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      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( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
1447               IF( scobot(ji,jj) == 0._wp            )   mbathy(ji,jj) = 0
1448            END DO
1449         END DO
1450      END DO
1451      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
1452         &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
1453
1454      IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain
1455         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) )
1456         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
1457            &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w '  , MINVAL( gdep3w_0(:,:,:) )
1458         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0   (:,:,:) ),   &
1459            &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0   (:,:,:) ),   &
1460            &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0  (:,:,:) ),   &
1461            &                          ' w ', MINVAL( e3w_0  (:,:,:) )
1462
1463         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
1464            &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w '  , MAXVAL( gdep3w_0(:,:,:) )
1465         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0   (:,:,:) ),   &
1466            &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0   (:,:,:) ),   &
1467            &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0  (:,:,:) ),   &
1468            &                          ' w ', MAXVAL( e3w_0  (:,:,:) )
1469      ENDIF
1470      !  END DO
1471      IF(lwp) THEN                                  ! selected vertical profiles
1472         WRITE(numout,*)
1473         WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
1474         WRITE(numout,*) ' ~~~~~~  --------------------'
1475         WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
1476         WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     &
1477            &                                 e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
1478         DO jj = mj0(20), mj1(20)
1479            DO ji = mi0(20), mi1(20)
1480               WRITE(numout,*)
1481               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
1482               WRITE(numout,*) ' ~~~~~~  --------------------'
1483               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
1484               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
1485                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
1486            END DO
1487         END DO
1488         DO jj = mj0(74), mj1(74)
1489            DO ji = mi0(100), mi1(100)
1490               WRITE(numout,*)
1491               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
1492               WRITE(numout,*) ' ~~~~~~  --------------------'
1493               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
1494               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
1495                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
1496            END DO
1497         END DO
1498      ENDIF
1499
1500!================================================================================
1501! check the coordinate makes sense
1502!================================================================================
1503      DO ji = 1, jpi
1504         DO jj = 1, jpj
1505
1506            IF( hbatt(ji,jj) > 0._wp) THEN
1507               DO jk = 1, mbathy(ji,jj)
1508                 ! check coordinate is monotonically increasing
1509                 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
1510                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
1511                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
1512                    WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
1513                    WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
1514                    CALL ctl_stop( ctmp1 )
1515                 ENDIF
1516                 ! and check it has never gone negative
1517                 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
1518                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
1519                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk
1520                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
1521                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
1522                    CALL ctl_stop( ctmp1 )
1523                 ENDIF
1524                 ! and check it never exceeds the total depth
1525                 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
1526                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
1527                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
1528                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
1529                    CALL ctl_stop( ctmp1 )
1530                 ENDIF
1531               END DO
1532
1533               DO jk = 1, mbathy(ji,jj)-1
1534                 ! and check it never exceeds the total depth
1535                IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
1536                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
1537                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
1538                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
1539                    CALL ctl_stop( ctmp1 )
1540                 ENDIF
1541               END DO
1542
1543            ENDIF
1544
1545         END DO
1546      END DO
1547      !
1548      CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
1549      !
1550      IF( nn_timing == 1 )  CALL timing_stop('zgr_sco')
1551      !
1552   END SUBROUTINE zgr_sco
1553
1554!!======================================================================
1555   SUBROUTINE s_sh94()
1556
1557      !!----------------------------------------------------------------------
1558      !!                  ***  ROUTINE s_sh94  ***
1559      !!                     
1560      !! ** Purpose :   stretch the s-coordinate system
1561      !!
1562      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994
1563      !!                mixed S/sigma coordinate
1564      !!
1565      !! Reference : Song and Haidvogel 1994.
1566      !!----------------------------------------------------------------------
1567      !
1568      INTEGER  ::   ji, jj, jk           ! dummy loop argument
1569      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
1570      !
1571      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
1572      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
1573
1574      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
1575      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1576
1577      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
1578      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
1579      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
1580      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
1581
1582      DO ji = 1, jpi
1583         DO jj = 1, jpj
1584
1585            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
1586               DO jk = 1, jpk
1587                  z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
1588                  z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
1589               END DO
1590            ELSE ! shallow water, uniform sigma
1591               DO jk = 1, jpk
1592                  z_gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
1593                  z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
1594                  END DO
1595            ENDIF
1596            !
1597            DO jk = 1, jpkm1
1598               z_esigt3(ji,jj,jk  ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
1599               z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
1600            END DO
1601            z_esigw3(ji,jj,1  ) = 2._wp * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) )
1602            z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
1603            !
1604            ! Coefficients for vertical depth as the sum of e3w scale factors
1605            z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
1606            DO jk = 2, jpk
1607               z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
1608            END DO
1609            !
1610            DO jk = 1, jpk
1611               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1612               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1613               gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
1614               gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
1615               gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
1616            END DO
1617           !
1618         END DO   ! for all jj's
1619      END DO    ! for all ji's
1620
1621      DO ji = 1, jpim1
1622         DO jj = 1, jpjm1
1623            DO jk = 1, jpk
1624               z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
1625                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1626               z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
1627                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1628               z_esigtf3(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+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   &
1630                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
1631               z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
1632                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1633               z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
1634                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1635               !
1636               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1637               e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1638               e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1639               e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1640               !
1641               e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1642               e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1643               e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1644            END DO
1645        END DO
1646      END DO
1647
1648      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
1649      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1650
1651   END SUBROUTINE s_sh94
1652
1653   SUBROUTINE s_sf12
1654
1655      !!----------------------------------------------------------------------
1656      !!                  ***  ROUTINE s_sf12 ***
1657      !!                     
1658      !! ** Purpose :   stretch the s-coordinate system
1659      !!
1660      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
1661      !!                mixed S/sigma/Z coordinate
1662      !!
1663      !!                This method allows the maintenance of fixed surface and or
1664      !!                bottom cell resolutions (cf. geopotential coordinates)
1665      !!                within an analytically derived stretched S-coordinate framework.
1666      !!
1667      !!
1668      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
1669      !!----------------------------------------------------------------------
1670      !
1671      INTEGER  ::   ji, jj, jk           ! dummy loop argument
1672      REAL(wp) ::   zsmth               ! smoothing around critical depth
1673      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
1674      !
1675      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
1676      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
1677
1678      !
1679      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
1680      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1681
1682      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
1683      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
1684      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
1685      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
1686
1687      DO ji = 1, jpi
1688         DO jj = 1, jpj
1689
1690          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
1691             
1692              zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
1693                                                     ! could be changed by users but care must be taken to do so carefully
1694              zzb = 1.0_wp-(zzb/hbatt(ji,jj))
1695           
1696              zzs = rn_zs / hbatt(ji,jj) 
1697             
1698              IF (rn_efold /= 0.0_wp) THEN
1699                zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
1700              ELSE
1701                zsmth = 1.0_wp 
1702              ENDIF
1703               
1704              DO jk = 1, jpk
1705                z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
1706                z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
1707              ENDDO
1708              z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth  )
1709              z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth  )
1710 
1711          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
1712
1713            DO jk = 1, jpk
1714              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
1715              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
1716            END DO
1717
1718          ELSE  ! shallow water, z coordinates
1719
1720            DO jk = 1, jpk
1721              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
1722              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
1723            END DO
1724
1725          ENDIF
1726
1727          DO jk = 1, jpkm1
1728             z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
1729             z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
1730          END DO
1731          z_esigw3(ji,jj,1  ) = 2.0_wp * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ))
1732          z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
1733
1734          ! Coefficients for vertical depth as the sum of e3w scale factors
1735          z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
1736          DO jk = 2, jpk
1737             z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
1738          END DO
1739
1740          DO jk = 1, jpk
1741             gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
1742             gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
1743             gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
1744          END DO
1745
1746        ENDDO   ! for all jj's
1747      ENDDO    ! for all ji's
1748
1749      DO ji=1,jpi-1
1750        DO jj=1,jpj-1
1751
1752          DO jk = 1, jpk
1753                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / &
1754                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1755                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / &
1756                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1757                z_esigtf3(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+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / &
1759                                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
1760                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / &
1761                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1762                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / &
1763                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1764
1765             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
1766             e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
1767             e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
1768             e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
1769             !
1770             e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
1771             e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
1772             e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
1773          END DO
1774
1775        ENDDO
1776      ENDDO
1777      !
1778      CALL lbc_lnk(e3t_0 ,'T',1.) ; CALL lbc_lnk(e3u_0 ,'T',1.)
1779      CALL lbc_lnk(e3v_0 ,'T',1.) ; CALL lbc_lnk(e3f_0 ,'T',1.)
1780      CALL lbc_lnk(e3w_0 ,'T',1.)
1781      CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.)
1782      !
1783      !                                               ! =============
1784
1785      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
1786      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1787
1788   END SUBROUTINE s_sf12
1789
1790   SUBROUTINE s_tanh()
1791
1792      !!----------------------------------------------------------------------
1793      !!                  ***  ROUTINE s_tanh***
1794      !!                     
1795      !! ** Purpose :   stretch the s-coordinate system
1796      !!
1797      !! ** Method  :   s-coordinate stretch
1798      !!
1799      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
1800      !!----------------------------------------------------------------------
1801
1802      INTEGER  ::   ji, jj, jk           ! dummy loop argument
1803      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
1804
1805      REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
1806      REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw
1807
1808      CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      )
1809      CALL wrk_alloc( jpk, z_esigt, z_esigw                                               )
1810
1811      z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp   ;   z_gsi3w  = 0._wp
1812      z_esigt  = 0._wp   ;   z_esigw  = 0._wp 
1813
1814      DO jk = 1, jpk
1815        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
1816        z_gsigt(jk) = -fssig( REAL(jk,wp)        )
1817      END DO
1818      IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk)
1819      !
1820      ! Coefficients for vertical scale factors at w-, t- levels
1821!!gm bug :  define it from analytical function, not like juste bellow....
1822!!gm        or betteroffer the 2 possibilities....
1823      DO jk = 1, jpkm1
1824         z_esigt(jk  ) = z_gsigw(jk+1) - z_gsigw(jk)
1825         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
1826      END DO
1827      z_esigw( 1 ) = 2._wp * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
1828      z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
1829      !
1830      ! Coefficients for vertical depth as the sum of e3w scale factors
1831      z_gsi3w(1) = 0.5_wp * z_esigw(1)
1832      DO jk = 2, jpk
1833         z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
1834      END DO
1835!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
1836      DO jk = 1, jpk
1837         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1838         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1839         gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
1840         gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
1841         gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )
1842      END DO
1843!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
1844      DO jj = 1, jpj
1845         DO ji = 1, jpi
1846            DO jk = 1, jpk
1847              e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1848              e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1849              e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1850              e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
1851              !
1852              e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1853              e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1854              e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1855            END DO
1856         END DO
1857      END DO
1858
1859      CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      )
1860      CALL wrk_dealloc( jpk, z_esigt, z_esigw                                               )
1861
1862   END SUBROUTINE s_tanh
1863
1864   FUNCTION fssig( pk ) RESULT( pf )
1865      !!----------------------------------------------------------------------
1866      !!                 ***  ROUTINE fssig ***
1867      !!       
1868      !! ** Purpose :   provide the analytical function in s-coordinate
1869      !!         
1870      !! ** Method  :   the function provide the non-dimensional position of
1871      !!                T and W (i.e. between 0 and 1)
1872      !!                T-points at integer values (between 1 and jpk)
1873      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1874      !!----------------------------------------------------------------------
1875      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
1876      REAL(wp)             ::   pf   ! sigma value
1877      !!----------------------------------------------------------------------
1878      !
1879      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   &
1880         &     - TANH( rn_thetb * rn_theta                                )  )   &
1881         & * (   COSH( rn_theta                           )                      &
1882         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              &
1883         & / ( 2._wp * SINH( rn_theta ) )
1884      !
1885   END FUNCTION fssig
1886
1887
1888   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
1889      !!----------------------------------------------------------------------
1890      !!                 ***  ROUTINE fssig1 ***
1891      !!
1892      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
1893      !!
1894      !! ** Method  :   the function provides the non-dimensional position of
1895      !!                T and W (i.e. between 0 and 1)
1896      !!                T-points at integer values (between 1 and jpk)
1897      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1898      !!----------------------------------------------------------------------
1899      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
1900      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
1901      REAL(wp)             ::   pf1   ! sigma value
1902      !!----------------------------------------------------------------------
1903      !
1904      IF ( rn_theta == 0 ) then      ! uniform sigma
1905         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
1906      ELSE                        ! stretched sigma
1907         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
1908            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
1909            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
1910      ENDIF
1911      !
1912   END FUNCTION fssig1
1913
1914
1915   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
1916      !!----------------------------------------------------------------------
1917      !!                 ***  ROUTINE fgamma  ***
1918      !!
1919      !! ** Purpose :   provide analytical function for the s-coordinate
1920      !!
1921      !! ** Method  :   the function provides the non-dimensional position of
1922      !!                T and W (i.e. between 0 and 1)
1923      !!                T-points at integer values (between 1 and jpk)
1924      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1925      !!
1926      !!                This method allows the maintenance of fixed surface and or
1927      !!                bottom cell resolutions (cf. geopotential coordinates)
1928      !!                within an analytically derived stretched S-coordinate framework.
1929      !!
1930      !! Reference  :   Siddorn and Furner, in prep
1931      !!----------------------------------------------------------------------
1932      REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate
1933      REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate
1934      REAL(wp), INTENT(in   ) ::   pzb           ! Bottom box depth
1935      REAL(wp), INTENT(in   ) ::   pzs           ! surface box depth
1936      REAL(wp), INTENT(in   ) ::   psmth       ! Smoothing parameter
1937      REAL(wp)                ::   za1,za2,za3    ! local variables
1938      REAL(wp)                ::   zn1,zn2        ! local variables
1939      REAL(wp)                ::   za,zb,zx       ! local variables
1940      integer                 ::   jk
1941      !!----------------------------------------------------------------------
1942      !
1943
1944      zn1  =  1./(jpk-1.)
1945      zn2  =  1. -  zn1
1946
1947      za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 
1948      za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp)
1949      za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
1950     
1951      za = pzb - za3*(pzs-za1)-za2
1952      za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) )
1953      zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
1954      zx = 1.0_wp-za/2.0_wp-zb
1955 
1956      DO jk = 1, jpk
1957        p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp +  &
1958                    & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &
1959                    &      (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
1960        p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth)
1961      ENDDO 
1962
1963      !
1964   END FUNCTION fgamma
1965
1966   !!======================================================================
1967END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.