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/2014/dev_r4765_CNRS_agrif/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2014/dev_r4765_CNRS_agrif/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 4982

Last change on this file since 4982 was 4982, checked in by jchanut, 9 years ago

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