source: branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5074

Last change on this file since 5074 was 5074, checked in by hliu, 6 years ago

prevent nla10 from being 0

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