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

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

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

Last change on this file since 5014 was 5014, checked in by hliu, 9 years ago

upload the modifications for W/D based on r:4826

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