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 NEMO/branches/2019/ENHANCE-02_ISF_domcfg/src – NEMO

source: NEMO/branches/2019/ENHANCE-02_ISF_domcfg/src/domzgr.F90 @ 11568

Last change on this file since 11568 was 11568, checked in by mathiot, 5 years ago

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