source: utils/tools/DOMAINcfg/src/domzgr.F90 @ 12414

Last change on this file since 12414 was 12414, checked in by smueller, 11 months ago

Reintegration of 2019 development branch /utils/tools_MERGE_2019 into the tools directory (/utils/tools)

File size: 98.1 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(:,:,:) )
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(:,:,:) )
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         DEALLOCATE( idta, zdta )
482         !
483         !                                            ! ================ !
484      ELSEIF( ntopo == 1 .OR. ntopo ==2 ) THEN                       !   read in file   ! (over the local domain)
485         !                                            ! ================ !
486         !
487         IF( ln_zco )   THEN                          ! zco : read level bathymetry
488            CALL iom_open ( 'bathy_level.nc', inum ) 
489            CALL iom_get  ( inum, jpdom_data, 'Bathy_level', bathy )
490            CALL iom_close( inum )
491            mbathy(:,:) = INT( bathy(:,:) )
492            ! initialisation isf variables
493            risfdep(:,:)=0._wp ; misfdep(:,:)=1             
494            !                                                ! =====================
495            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
496               !                                             ! =====================
497               !
498               ii0 = 140   ;   ii1 = 140                  ! Gibraltar Strait open
499               ij0 = 102   ;   ij1 = 102                  ! (Thomson, Ocean Modelling, 1995)
500               DO ji = mi0(ii0), mi1(ii1)
501                  DO jj = mj0(ij0), mj1(ij1)
502                     mbathy(ji,jj) = 15
503                  END DO
504               END DO
505               IF(lwp) WRITE(numout,*)
506               IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
507               !
508               ii0 = 160   ;   ii1 = 160                  ! Bab el mandeb Strait open
509               ij0 = 88    ;   ij1 = 88                   ! (Thomson, Ocean Modelling, 1995)
510               DO ji = mi0(ii0), mi1(ii1)
511                  DO jj = mj0(ij0), mj1(ij1)
512                     mbathy(ji,jj) = 12
513                  END DO
514               END DO
515               IF(lwp) WRITE(numout,*)
516               IF(lwp) WRITE(numout,*) '      orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
517               !
518            ENDIF
519            !
520         ENDIF
521         IF( ln_zps .OR. ln_sco )   THEN              ! zps or sco : read meter bathymetry
522#if defined key_agrif
523            IF (agrif_root()) THEN
524#endif
525            IF( ntopo == 1) THEN
526               CALL iom_open ( 'bathy_meter.nc', inum ) 
527               IF ( ln_isfcav ) THEN
528                  CALL iom_get  ( inum, jpdom_data, 'Bathymetry_isf', bathy, lrowattr=.false. )
529               ELSE
530                  CALL iom_get  ( inum, jpdom_data, 'Bathymetry'    , bathy, lrowattr=ln_use_jattr  )
531               END IF
532               CALL iom_close( inum )
533            ELSE
534               CALL dom_bat
535            ENDIF       
536#if defined key_agrif
537            ELSE
538               IF( ntopo == 1) THEN
539                  CALL agrif_create_bathy_meter()
540               ELSE
541                  CALL dom_bat 
542               ENDIF   
543            ENDIF
544#endif
545            !                                               
546            ! initialisation isf variables
547            risfdep(:,:)=0._wp ; misfdep(:,:)=1             
548            !
549            IF ( ln_isfcav ) THEN
550               CALL iom_open ( 'isf_draft_meter.nc', inum ) 
551               CALL iom_get  ( inum, jpdom_data, 'isf_draft', risfdep )
552               CALL iom_close( inum )
553            END IF
554            !       
555            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
556            !
557              ii0 = 140   ;   ii1 = 140                   ! Gibraltar Strait open
558              ij0 = 102   ;   ij1 = 102                   ! (Thomson, Ocean Modelling, 1995)
559              DO ji = mi0(ii0), mi1(ii1)
560                 DO jj = mj0(ij0), mj1(ij1)
561                    bathy(ji,jj) = 284._wp
562                 END DO
563               END DO
564              IF(lwp) WRITE(numout,*)     
565              IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
566              !
567              ii0 = 160   ;   ii1 = 160                   ! Bab el mandeb Strait open
568              ij0 = 88    ;   ij1 = 88                    ! (Thomson, Ocean Modelling, 1995)
569               DO ji = mi0(ii0), mi1(ii1)
570                 DO jj = mj0(ij0), mj1(ij1)
571                    bathy(ji,jj) = 137._wp
572                 END DO
573              END DO
574              IF(lwp) WRITE(numout,*)
575               IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
576              !
577           ENDIF
578            !
579        ENDIF
580         !                                            ! =============== !
581      ELSE                                            !      error      !
582         !                                            ! =============== !
583         WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
584         CALL ctl_stop( '    zgr_bat : '//trim(ctmp1) )
585      ENDIF
586      !
587      IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==!
588         IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level
589         ELSE                          ;   ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 )  ! from a depth
590         ENDIF
591         zhmin = gdepw_1d(ik+1)                                                        ! minimum depth = ik+1 w-levels
592         WHERE( bathy(:,:) <= 0._wp )   ;   bathy(:,:) = 0._wp                         ! min=0     over the lands
593         ELSE WHERE ( risfdep == 0._wp );   bathy(:,:) = MAX(  zhmin , bathy(:,:)  )   ! min=zhmin over the oceans
594         END WHERE
595         IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
596      ENDIF
597      !
598   END SUBROUTINE zgr_bat
599
600   SUBROUTINE zgr_bat_ctl
601      !!----------------------------------------------------------------------
602      !!                    ***  ROUTINE zgr_bat_ctl  ***
603      !!
604      !! ** Purpose :   check the bathymetry in levels
605      !!
606      !! ** Method  :   The array mbathy is checked to verified its consistency
607      !!      with the model options. in particular:
608      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
609      !!                  along closed boundary.
610      !!            mbathy must be cyclic IF jperio=1.
611      !!            mbathy must be lower or equal to jpk-1.
612      !!            isolated ocean grid points are suppressed from mbathy
613      !!                  since they are only connected to remaining
614      !!                  ocean through vertical diffusion.
615      !!      C A U T I O N : mbathy will be modified during the initializa-
616      !!      tion phase to become the number of non-zero w-levels of a water
617      !!      column, with a minimum value of 1.
618      !!
619      !! ** Action  : - update mbathy: level bathymetry (in level index)
620      !!              - update bathy : meter bathymetry (in meters)
621      !!----------------------------------------------------------------------
622      INTEGER ::   ji, jj, jl                    ! dummy loop indices
623      INTEGER ::   icompt, ibtest, ikmax         ! temporary integers
624      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  zbathy
625      !!----------------------------------------------------------------------
626      !
627      ALLOCATE(zbathy(jpi,jpj))
628      !
629      IF(lwp) WRITE(numout,*)
630      IF(lwp) WRITE(numout,*) '    zgr_bat_ctl : check the bathymetry'
631      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
632      !                                          ! Suppress isolated ocean grid points
633      IF(lwp) WRITE(numout,*)
634      IF(lwp) WRITE(numout,*)'                   suppress isolated ocean grid points'
635      IF(lwp) WRITE(numout,*)'                   -----------------------------------'
636      icompt = 0
637      DO jl = 1, 2
638         IF( l_Iperio ) THEN
639            mbathy( 1 ,:) = mbathy(jpim1,:)           ! local domain is cyclic east-west
640            mbathy(jpi,:) = mbathy(  2  ,:)
641         ENDIF
642         zbathy(:,:) = FLOAT( mbathy(:,:) )
643         CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp )
644         mbathy(:,:) = INT( zbathy(:,:) )
645         
646         DO jj = 2, jpjm1
647            DO ji = 2, jpim1
648               ibtest = MAX(  mbathy(ji-1,jj), mbathy(ji+1,jj),   &
649                  &           mbathy(ji,jj-1), mbathy(ji,jj+1)  )
650               IF( ibtest < mbathy(ji,jj) ) THEN
651                  IF(lwp) WRITE(numout,*) ' the number of ocean level at ',   &
652                     &   'grid-point (i,j) =  ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
653                  mbathy(ji,jj) = ibtest
654                  icompt = icompt + 1
655               ENDIF
656            END DO
657         END DO
658      END DO
659
660      IF( lk_mpp )   CALL mpp_sum( 'domzgr', icompt )
661      IF( icompt == 0 ) THEN
662         IF(lwp) WRITE(numout,*)'     no isolated ocean grid points'
663      ELSE
664         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed'
665      ENDIF
666
667      zbathy(:,:) = FLOAT( mbathy(:,:) )
668      CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp )
669      mbathy(:,:) = INT( zbathy(:,:) )
670
671      !                                          ! East-west cyclic boundary conditions
672      IF( jperio == 0 ) THEN
673         IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: jperio = ', jperio
674         IF( lk_mpp ) THEN
675            IF( nbondi == -1 .OR. nbondi == 2 ) THEN
676               IF( jperio /= 1 )   mbathy(1,:) = 0
677            ENDIF
678            IF( nbondi == 1 .OR. nbondi == 2 ) THEN
679               IF( jperio /= 1 )   mbathy(nlci,:) = 0
680            ENDIF
681         ELSE
682            IF( ln_zco .OR. ln_zps ) THEN
683               mbathy( 1 ,:) = 0
684               mbathy(jpi,:) = 0
685            ELSE
686               mbathy( 1 ,:) = jpkm1
687               mbathy(jpi,:) = jpkm1
688            ENDIF
689         ENDIF
690      ELSEIF( l_Iperio ) THEN
691         IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: jperio = ', jperio
692         mbathy( 1 ,:) = mbathy(jpim1,:)
693         mbathy(jpi,:) = mbathy(  2  ,:)
694      ELSEIF( jperio == 2 ) THEN
695         IF(lwp) WRITE(numout,*) '   equatorial boundary conditions on mbathy: jperio = ', jperio
696      ELSE
697         IF(lwp) WRITE(numout,*) '    e r r o r'
698         IF(lwp) WRITE(numout,*) '    parameter , jperio = ', jperio
699         !         STOP 'dom_mba'
700      ENDIF
701
702      !  Boundary condition on mbathy
703      IF( .NOT.lk_mpp ) THEN 
704!!gm     !!bug ???  think about it !
705         !   ... mono- or macro-tasking: T-point, >0, 2D array, no slab
706         zbathy(:,:) = FLOAT( mbathy(:,:) )
707         CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp )
708         mbathy(:,:) = INT( zbathy(:,:) )
709      ENDIF
710
711      ! Number of ocean level inferior or equal to jpkm1
712      zbathy(:,:) = FLOAT( mbathy(:,:) )
713      ikmax = glob_max( 'domzgr', zbathy(:,:) )
714
715      IF( ikmax > jpkm1 ) THEN
716         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' >  jpk-1'
717         IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
718      ELSE IF( ikmax < jpkm1 ) THEN
719         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1' 
720         IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
721      ENDIF
722      !
723      DEALLOCATE( zbathy )
724      !
725   END SUBROUTINE zgr_bat_ctl
726
727
728   SUBROUTINE zgr_bot_level
729      !!----------------------------------------------------------------------
730      !!                    ***  ROUTINE zgr_bot_level  ***
731      !!
732      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
733      !!
734      !! ** Method  :   computes from mbathy with a minimum value of 1 over land
735      !!
736      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest
737      !!                                     ocean level at t-, u- & v-points
738      !!                                     (min value = 1 over land)
739      !!----------------------------------------------------------------------
740      INTEGER ::   ji, jj   ! dummy loop indices
741      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  zmbk
742      !!----------------------------------------------------------------------
743      !
744      ALLOCATE( zmbk(jpi,jpj) )
745      !
746      IF(lwp) WRITE(numout,*)
747      IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
748      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
749      !
750      mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
751 
752      !                                     ! bottom k-index of W-level = mbkt+1
753      DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level
754         DO ji = 1, jpim1
755            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
756            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
757         END DO
758      END DO
759      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
760      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk('domzgr',zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
761      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk('domzgr',zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
762      !
763      DEALLOCATE( zmbk )
764      !
765   END SUBROUTINE zgr_bot_level
766
767
768   SUBROUTINE zgr_top_level
769      !!----------------------------------------------------------------------
770      !!                    ***  ROUTINE zgr_top_level  ***
771      !!
772      !! ** Purpose :   defines the vertical index of ocean top (mik. arrays)
773      !!
774      !! ** Method  :   computes from misfdep with a minimum value of 1
775      !!
776      !! ** Action  :   mikt, miku, mikv :   vertical indices of the shallowest
777      !!                                     ocean level at t-, u- & v-points
778      !!                                     (min value = 1)
779      !!----------------------------------------------------------------------
780      INTEGER ::   ji, jj   ! dummy loop indices
781      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  zmik
782      !!----------------------------------------------------------------------
783      !
784      ALLOCATE( zmik(jpi,jpj) )
785      !
786      IF(lwp) WRITE(numout,*)
787      IF(lwp) WRITE(numout,*) '    zgr_top_level : ocean top k-index of T-, U-, V- and W-levels '
788      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
789      !
790      mikt(:,:) = MAX( misfdep(:,:) , 1 )    ! top k-index of T-level (=1)
791      !                                      ! top k-index of W-level (=mikt)
792      DO jj = 1, jpjm1                       ! top k-index of U- (U-) level
793         DO ji = 1, jpim1
794            miku(ji,jj) = MAX(  mikt(ji+1,jj  ) , mikt(ji,jj)  )
795            mikv(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj)  )
796            mikf(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj), mikt(ji+1,jj  ), mikt(ji+1,jj+1)  )
797         END DO
798      END DO
799
800      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
801      zmik(:,:) = REAL( miku(:,:), wp )   ;   CALL lbc_lnk('domzgr',zmik,'U',1.)   ;   miku  (:,:) = MAX( INT( zmik(:,:) ), 1 )
802      zmik(:,:) = REAL( mikv(:,:), wp )   ;   CALL lbc_lnk('domzgr',zmik,'V',1.)   ;   mikv  (:,:) = MAX( INT( zmik(:,:) ), 1 )
803      zmik(:,:) = REAL( mikf(:,:), wp )   ;   CALL lbc_lnk('domzgr',zmik,'F',1.)   ;   mikf  (:,:) = MAX( INT( zmik(:,:) ), 1 )
804      !
805      DEALLOCATE( zmik )
806      !
807   END SUBROUTINE zgr_top_level
808
809
810   SUBROUTINE zgr_zco
811      !!----------------------------------------------------------------------
812      !!                  ***  ROUTINE zgr_zco  ***
813      !!
814      !! ** Purpose :   define the reference z-coordinate system
815      !!
816      !! ** Method  :   set 3D coord. arrays to reference 1D array
817      !!----------------------------------------------------------------------
818      INTEGER  ::   jk
819      !!----------------------------------------------------------------------
820      !
821      DO jk = 1, jpk
822         gdept_0(:,:,jk) = gdept_1d(jk)
823         gdepw_0(:,:,jk) = gdepw_1d(jk)
824         e3t_0  (:,:,jk) = e3t_1d  (jk)
825         e3u_0  (:,:,jk) = e3t_1d  (jk)
826         e3v_0  (:,:,jk) = e3t_1d  (jk)
827         e3f_0  (:,:,jk) = e3t_1d  (jk)
828         e3w_0  (:,:,jk) = e3w_1d  (jk)
829         e3uw_0 (:,:,jk) = e3w_1d  (jk)
830         e3vw_0 (:,:,jk) = e3w_1d  (jk)
831      END DO
832      !
833   END SUBROUTINE zgr_zco
834
835
836   SUBROUTINE zgr_zps
837      !!----------------------------------------------------------------------
838      !!                  ***  ROUTINE zgr_zps  ***
839      !!                     
840      !! ** Purpose :   the depth and vertical scale factor in partial step
841      !!              reference z-coordinate case
842      !!
843      !! ** Method  :   Partial steps : computes the 3D vertical scale factors
844      !!      of T-, U-, V-, W-, UW-, VW and F-points that are associated with
845      !!      a partial step representation of bottom topography.
846      !!
847      !!        The reference depth of model levels is defined from an analytical
848      !!      function the derivative of which gives the reference vertical
849      !!      scale factors.
850      !!        From  depth and scale factors reference, we compute there new value
851      !!      with partial steps  on 3d arrays ( i, j, k ).
852      !!
853      !!              w-level: gdepw_0(i,j,k)  = gdep(k)
854      !!                       e3w_0(i,j,k) = dk(gdep)(k)     = e3(i,j,k)
855      !!              t-level: gdept_0(i,j,k)  = gdep(k+0.5)
856      !!                       e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5)
857      !!
858      !!        With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
859      !!      we find the mbathy index of the depth at each grid point.
860      !!      This leads us to three cases:
861      !!
862      !!              - bathy = 0 => mbathy = 0
863      !!              - 1 < mbathy < jpkm1   
864      !!              - bathy > gdepw_0(jpk) => mbathy = jpkm1 
865      !!
866      !!        Then, for each case, we find the new depth at t- and w- levels
867      !!      and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
868      !!      and f-points.
869      !!
870      !!        This routine is given as an example, it must be modified
871      !!      following the user s desiderata. nevertheless, the output as
872      !!      well as the way to compute the model levels and scale factors
873      !!      must be respected in order to insure second order accuracy
874      !!      schemes.
875      !!
876      !!         c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives
877      !!         - - - - - - -   gdept_0, gdepw_0 and e3. are positives
878      !!     
879      !!  Reference :   Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
880      !!----------------------------------------------------------------------
881      INTEGER  ::   ji, jj, jk       ! dummy loop indices
882      INTEGER  ::   ik, it, ikb, ikt ! temporary integers
883      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
884      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
885      REAL(wp) ::   zdiff            ! temporary scalar
886      REAL(wp) ::   zmax             ! temporary scalar
887      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  zprt
888      !!---------------------------------------------------------------------
889      !
890      ALLOCATE( zprt(jpi,jpj,jpk) )
891      !
892      IF(lwp) WRITE(numout,*)
893      IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps'
894      IF(lwp) WRITE(numout,*) '    ~~~~~~~ '
895      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used'
896
897      ! compute position of the ice shelf grounding line
898      ! set bathy and isfdraft to 0 where grounded
899      IF ( ln_isfcav ) CALL zgr_isf_zspace
900
901      ! bathymetry in level (from bathy_meter)
902      ! ===================
903      zmax = gdepw_1d(jpk) + e3t_1d(jpk)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) )
904      bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat)
905      WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0
906      ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level
907      END WHERE
908
909      ! Compute mbathy for ocean points (i.e. the number of ocean levels)
910      ! find the number of ocean levels such that the last level thickness
911      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where
912      ! e3t_1d is the reference level thickness
913      DO jk = jpkm1, 1, -1
914         zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
915         WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1
916      END DO
917
918      ! Check compatibility between bathy and iceshelf draft
919      ! insure at least 2 wet level on the vertical under an ice shelf
920      ! compute misfdep and adjust isf draft if needed
921      IF ( ln_isfcav ) CALL zgr_isf_kspace
922
923      ! Scale factors and depth at T- and W-points
924      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate
925         gdept_0(:,:,jk) = gdept_1d(jk)
926         gdepw_0(:,:,jk) = gdepw_1d(jk)
927         e3t_0  (:,:,jk) = e3t_1d  (jk)
928         e3w_0  (:,:,jk) = e3w_1d  (jk)
929      END DO
930     
931      ! Scale factors and depth at T- and W-points
932      DO jj = 1, jpj
933         DO ji = 1, jpi
934            ik = mbathy(ji,jj)
935            IF( ik > 0 ) THEN               ! ocean point only
936               ! max ocean level case
937               IF( ik == jpkm1 ) THEN
938                  zdepwp = bathy(ji,jj)
939                  ze3tp  = bathy(ji,jj) - gdepw_1d(ik)
940                  ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
941                  e3t_0(ji,jj,ik  ) = ze3tp
942                  e3t_0(ji,jj,ik+1) = ze3tp
943                  e3w_0(ji,jj,ik  ) = ze3wp
944                  e3w_0(ji,jj,ik+1) = ze3tp
945                  gdepw_0(ji,jj,ik+1) = zdepwp
946                  gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp
947                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
948                  !
949               ELSE                         ! standard case
950                  IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
951                  ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
952                  ENDIF
953!gm Bug?  check the gdepw_1d
954                  !       ... on ik
955                  gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   &
956                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   &
957                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) ))
958                  e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   & 
959                     &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ) 
960                  e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   &
961                     &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
962                  !       ... on ik+1
963                  e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
964                  e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
965                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
966               ENDIF
967            ENDIF
968         END DO
969      END DO
970      !
971      it = 0
972      DO jj = 1, jpj
973         DO ji = 1, jpi
974            ik = mbathy(ji,jj)
975            IF( ik > 0 ) THEN               ! ocean point only
976               e3tp (ji,jj) = e3t_0(ji,jj,ik)
977               e3wp (ji,jj) = e3w_0(ji,jj,ik)
978               ! test
979               zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  )
980               IF( zdiff <= 0._wp .AND. lwp ) THEN
981                  it = it + 1
982                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj
983                  WRITE(numout,*) ' bathy = ', bathy(ji,jj)
984                  WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
985                  WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  )
986               ENDIF
987            ENDIF
988         END DO
989      END DO
990      !
991      ! compute top scale factor if ice shelf
992      IF (ln_isfcav) CALL zps_isf
993      !
994      ! Scale factors and depth at U-, V-, UW and VW-points
995      DO jk = 1, jpk                        ! initialisation to z-scale factors
996         e3u_0 (:,:,jk) = e3t_1d(jk)
997         e3v_0 (:,:,jk) = e3t_1d(jk)
998         e3uw_0(:,:,jk) = e3w_1d(jk)
999         e3vw_0(:,:,jk) = e3w_1d(jk)
1000      END DO
1001
1002      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
1003         DO jj = 1, jpjm1
1004            DO ji = 1, jpim1   ! vector opt.
1005               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
1006               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
1007               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
1008               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
1009            END DO
1010         END DO
1011      END DO
1012
1013      ! update e3uw in case only 2 cells in the water column
1014      IF ( ln_isfcav ) CALL zps_isf_e3uv_w
1015      !
1016      CALL lbc_lnk('domzgr', e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk('domzgr', e3uw_0, 'U', 1._wp )   ! lateral boundary conditions
1017      CALL lbc_lnk('domzgr', e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp )
1018      !
1019      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
1020         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk)
1021         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk)
1022         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk)
1023         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk)
1024      END DO
1025     
1026      ! Scale factor at F-point
1027      DO jk = 1, jpk                        ! initialisation to z-scale factors
1028         e3f_0(:,:,jk) = e3t_1d(jk)
1029      END DO
1030      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors
1031         DO jj = 1, jpjm1
1032            DO ji = 1, jpim1   ! vector opt.
1033               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
1034            END DO
1035         END DO
1036      END DO
1037      CALL lbc_lnk('domzgr', e3f_0, 'F', 1._wp )       ! Lateral boundary conditions
1038      !
1039      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
1040         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk)
1041      END DO
1042!!gm  bug ? :  must be a do loop with mj0,mj1
1043      !
1044      e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2
1045      e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 
1046      e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 
1047      e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 
1048      e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 
1049
1050      ! Control of the sign
1051      IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' )
1052      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' )
1053      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' )
1054      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' )
1055      !
1056      ! if in the future gde3w_0 need to be compute, use the function defined in NEMO
1057      ! for now gde3w_0 computation is removed as not an output of domcfg
1058
1059      DEALLOCATE( zprt )
1060      !
1061   END SUBROUTINE zgr_zps
1062
1063   SUBROUTINE zgr_sco
1064      !!----------------------------------------------------------------------
1065      !!                  ***  ROUTINE zgr_sco  ***
1066      !!                     
1067      !! ** Purpose :   define the s-coordinate system
1068      !!
1069      !! ** Method  :   s-coordinate
1070      !!         The depth of model levels is defined as the product of an
1071      !!      analytical function by the local bathymetry, while the vertical
1072      !!      scale factors are defined as the product of the first derivative
1073      !!      of the analytical function by the bathymetry.
1074      !!      (this solution save memory as depth and scale factors are not
1075      !!      3d fields)
1076      !!          - Read bathymetry (in meters) at t-point and compute the
1077      !!         bathymetry at u-, v-, and f-points.
1078      !!            hbatu = mi( hbatt )
1079      !!            hbatv = mj( hbatt )
1080      !!            hbatf = mi( mj( hbatt ) )
1081      !!          - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
1082      !!         function and its derivative given as function.
1083      !!            z_gsigt(k) = fssig (k    )
1084      !!            z_gsigw(k) = fssig (k-0.5)
1085      !!            z_esigt(k) = fsdsig(k    )
1086      !!            z_esigw(k) = fsdsig(k-0.5)
1087      !!      Three options for stretching are give, and they can be modified
1088      !!      following the users requirements. Nevertheless, the output as
1089      !!      well as the way to compute the model levels and scale factors
1090      !!      must be respected in order to insure second order accuracy
1091      !!      schemes.
1092      !!
1093      !!      The three methods for stretching available are:
1094      !!
1095      !!           s_sh94 (Song and Haidvogel 1994)
1096      !!                a sinh/tanh function that allows sigma and stretched sigma
1097      !!
1098      !!           s_sf12 (Siddorn and Furner 2012?)
1099      !!                allows the maintenance of fixed surface and or
1100      !!                bottom cell resolutions (cf. geopotential coordinates)
1101      !!                within an analytically derived stretched S-coordinate framework.
1102      !!
1103      !!          s_tanh  (Madec et al 1996)
1104      !!                a cosh/tanh function that gives stretched coordinates       
1105      !!
1106      !!----------------------------------------------------------------------
1107      INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
1108      INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
1109      INTEGER  ::   ios                      ! Local integer output status for namelist read
1110      REAL(wp) ::   zrmax, ztaper   ! temporary scalars
1111      REAL(wp) ::   zrfact
1112      !
1113      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
1114      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
1115      !!
1116      NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
1117         &                rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
1118     !!----------------------------------------------------------------------
1119      !
1120      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) )
1121      !
1122      REWIND( numnam_ref )              ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
1123      READ  ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901)
1124901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist', lwp )
1125
1126      REWIND( numnam_cfg )              ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters
1127      READ  ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 )
1128902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist', lwp )
1129      IF(lwm) WRITE ( numond, namzgr_sco )
1130
1131      IF(lwp) THEN                           ! control print
1132         WRITE(numout,*)
1133         WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
1134         WRITE(numout,*) '~~~~~~~~~~~'
1135         WRITE(numout,*) '   Namelist namzgr_sco'
1136         WRITE(numout,*) '     stretching coeffs '
1137         WRITE(numout,*) '        maximum depth of s-bottom surface (>0)       rn_sbot_max   = ',rn_sbot_max
1138         WRITE(numout,*) '        minimum depth of s-bottom surface (>0)       rn_sbot_min   = ',rn_sbot_min
1139         WRITE(numout,*) '        Critical depth                               rn_hc         = ',rn_hc
1140         WRITE(numout,*) '        maximum cut-off r-value allowed              rn_rmax       = ',rn_rmax
1141         WRITE(numout,*) '     Song and Haidvogel 1994 stretching              ln_s_sh94     = ',ln_s_sh94
1142         WRITE(numout,*) '        Song and Haidvogel 1994 stretching coefficients'
1143         WRITE(numout,*) '        surface control parameter (0<=rn_theta<=20)  rn_theta      = ',rn_theta
1144         WRITE(numout,*) '        bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ',rn_thetb
1145         WRITE(numout,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb
1146         WRITE(numout,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12
1147         WRITE(numout,*) '        switching to sigma (T) or Z (F) at H<Hc      ln_sigcrit    = ',ln_sigcrit
1148         WRITE(numout,*) '        Siddorn and Furner 2012 stretching coefficients'
1149         WRITE(numout,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha
1150         WRITE(numout,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold
1151         WRITE(numout,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs
1152         WRITE(numout,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a
1153         WRITE(numout,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b
1154         WRITE(numout,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b'
1155      ENDIF
1156
1157      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate
1158      hifu(:,:) = rn_sbot_min
1159      hifv(:,:) = rn_sbot_min
1160      hiff(:,:) = rn_sbot_min
1161
1162      !                                        ! set maximum ocean depth
1163      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1164
1165         DO jj = 1, jpj
1166            DO ji = 1, jpi
1167              IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1168            END DO
1169         END DO
1170      !                                        ! =============================
1171      !                                        ! Define the envelop bathymetry   (hbatt)
1172      !                                        ! =============================
1173      ! use r-value to create hybrid coordinates
1174      zenv(:,:) = bathy(:,:)
1175      !
1176      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
1177         DO jj = 1, jpj
1178            DO ji = 1, jpi
1179               IF( bathy(ji,jj) == 0._wp ) THEN
1180                  iip1 = MIN( ji+1, jpi )
1181                  ijp1 = MIN( jj+1, jpj )
1182                  iim1 = MAX( ji-1, 1 )
1183                  ijm1 = MAX( jj-1, 1 )
1184!!gm BUG fix see ticket #1617
1185                  IF( ( + bathy(iim1,ijm1) + bathy(ji,ijp1) + bathy(iip1,ijp1)              &
1186                     &  + bathy(iim1,jj  )                  + bathy(iip1,jj  )              &
1187                     &  + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1)  ) > 0._wp ) &
1188                     &    zenv(ji,jj) = rn_sbot_min
1189!!gm
1190!!gm               IF( ( bathy(iip1,jj  ) + bathy(iim1,jj  ) + bathy(ji,ijp1  ) + bathy(ji,ijm1) +         &
1191!!gm                  &  bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
1192!!gm               zenv(ji,jj) = rn_sbot_min
1193!!gm             ENDIF
1194!!gm end
1195              ENDIF
1196            END DO
1197         END DO
1198
1199      ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1200      CALL lbc_lnk( 'domzgr',zenv, 'T', 1._wp, 'no0' )
1201      !
1202      ! smooth the bathymetry (if required)
1203      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
1204      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
1205      !
1206      jl = 0
1207      zrmax = 1._wp
1208      !   
1209      !     
1210      ! set scaling factor used in reducing vertical gradients
1211      zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
1212      !
1213      ! initialise temporary evelope depth arrays
1214      ztmpi1(:,:) = zenv(:,:)
1215      ztmpi2(:,:) = zenv(:,:)
1216      ztmpj1(:,:) = zenv(:,:)
1217      ztmpj2(:,:) = zenv(:,:)
1218      !
1219      ! initialise temporary r-value arrays
1220      zri(:,:) = 1._wp
1221      zrj(:,:) = 1._wp
1222      !                                                            ! ================ !
1223      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) !  Iterative loop  !
1224         !                                                         ! ================ !
1225         jl = jl + 1
1226         zrmax = 0._wp
1227         ! we set zrmax from previous r-values (zri and zrj) first
1228         ! if set after current r-value calculation (as previously)
1229         ! we could exit DO WHILE prematurely before checking r-value
1230         ! of current zenv
1231         DO jj = 1, nlcj
1232            DO ji = 1, nlci
1233               zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
1234            END DO
1235         END DO
1236         zri(:,:) = 0._wp
1237         zrj(:,:) = 0._wp
1238         DO jj = 1, nlcj
1239            DO ji = 1, nlci
1240               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
1241               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
1242               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
1243                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
1244               END IF
1245               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
1246                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
1247               END IF
1248               IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact
1249               IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact
1250               IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact
1251               IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact
1252            END DO
1253         END DO
1254  !       IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
1255         !
1256         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax
1257         !
1258         DO jj = 1, nlcj
1259            DO ji = 1, nlci
1260               zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
1261            END DO
1262         END DO
1263         ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1264         CALL lbc_lnk( 'domzgr',zenv, 'T', 1._wp, 'no0' )
1265         !                                                  ! ================ !
1266      END DO                                                !     End loop     !
1267      !                                                     ! ================ !
1268      DO jj = 1, jpj
1269         DO ji = 1, jpi
1270            zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
1271         END DO
1272      END DO
1273      !
1274      ! Envelope bathymetry saved in hbatt
1275      hbatt(:,:) = zenv(:,:) 
1276      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
1277         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1278         DO jj = 1, jpj
1279            DO ji = 1, jpi
1280               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
1281               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
1282            END DO
1283         END DO
1284      ENDIF
1285      !
1286      !                                        ! ==============================
1287      !                                        !   hbatu, hbatv, hbatf fields
1288      !                                        ! ==============================
1289      IF(lwp) THEN
1290         WRITE(numout,*)
1291           WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
1292      ENDIF
1293      hbatu(:,:) = rn_sbot_min
1294      hbatv(:,:) = rn_sbot_min
1295      hbatf(:,:) = rn_sbot_min
1296      DO jj = 1, jpjm1
1297        DO ji = 1, jpim1   ! NO vector opt.
1298           hbatu(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
1299           hbatv(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
1300           hbatf(ji,jj) = 0.25_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
1301              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
1302        END DO
1303      END DO
1304
1305      !
1306      ! Apply lateral boundary condition
1307!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
1308      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk('domzgr', hbatu, 'U', 1._wp )
1309      DO jj = 1, jpj
1310         DO ji = 1, jpi
1311            IF( hbatu(ji,jj) == 0._wp ) THEN
1312               !No worries about the following line when ln_wd == .true.
1313               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
1314               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
1315            ENDIF
1316         END DO
1317      END DO
1318      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk('domzgr', hbatv, 'V', 1._wp )
1319      DO jj = 1, jpj
1320         DO ji = 1, jpi
1321            IF( hbatv(ji,jj) == 0._wp ) THEN
1322               IF( zhbat(ji,jj) == 0._wp )   hbatv(ji,jj) = rn_sbot_min
1323               IF( zhbat(ji,jj) /= 0._wp )   hbatv(ji,jj) = zhbat(ji,jj)
1324            ENDIF
1325         END DO
1326      END DO
1327      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk('domzgr', hbatf, 'F', 1._wp )
1328      DO jj = 1, jpj
1329         DO ji = 1, jpi
1330            IF( hbatf(ji,jj) == 0._wp ) THEN
1331               IF( zhbat(ji,jj) == 0._wp )   hbatf(ji,jj) = rn_sbot_min
1332               IF( zhbat(ji,jj) /= 0._wp )   hbatf(ji,jj) = zhbat(ji,jj)
1333            ENDIF
1334         END DO
1335      END DO
1336
1337!!bug:  key_helsinki a verifer
1338        hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1339        hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1340        hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1341        hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1342
1343      IF( nprint == 1 .AND. lwp )   THEN
1344         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
1345            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
1346         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
1347            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
1348         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
1349            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1350         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
1351            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1352      ENDIF
1353!! helsinki
1354
1355      !                                            ! =======================
1356      !                                            !   s-ccordinate fields     (gdep., e3.)
1357      !                                            ! =======================
1358      !
1359      ! non-dimensional "sigma" for model level depth at w- and t-levels
1360
1361
1362!========================================================================
1363! Song and Haidvogel  1994 (ln_s_sh94=T)
1364! Siddorn and Furner 2012 (ln_sf12=T)
1365! or  tanh function       (both false)                   
1366!========================================================================
1367      IF      ( ln_s_sh94 ) THEN
1368                           CALL s_sh94()
1369      ELSE IF ( ln_s_sf12 ) THEN
1370                           CALL s_sf12()
1371      ELSE                 
1372                           CALL s_tanh()
1373      ENDIF
1374
1375      CALL lbc_lnk( 'domzgr',e3t_0 , 'T', 1._wp )
1376      CALL lbc_lnk( 'domzgr',e3u_0 , 'U', 1._wp )
1377      CALL lbc_lnk( 'domzgr',e3v_0 , 'V', 1._wp )
1378      CALL lbc_lnk( 'domzgr',e3f_0 , 'F', 1._wp )
1379      CALL lbc_lnk( 'domzgr',e3w_0 , 'W', 1._wp )
1380      CALL lbc_lnk( 'domzgr',e3uw_0, 'U', 1._wp )
1381      CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp )
1382      !
1383        WHERE( e3t_0 (:,:,:) == 0._wp )   e3t_0 (:,:,:) = 1._wp
1384        WHERE( e3u_0 (:,:,:) == 0._wp )   e3u_0 (:,:,:) = 1._wp
1385        WHERE( e3v_0 (:,:,:) == 0._wp )   e3v_0 (:,:,:) = 1._wp
1386        WHERE( e3f_0 (:,:,:) == 0._wp )   e3f_0 (:,:,:) = 1._wp
1387        WHERE( e3w_0 (:,:,:) == 0._wp )   e3w_0 (:,:,:) = 1._wp
1388        WHERE( e3uw_0(:,:,:) == 0._wp )   e3uw_0(:,:,:) = 1._wp
1389        WHERE( e3vw_0(:,:,:) == 0._wp )   e3vw_0(:,:,:) = 1._wp
1390!!
1391      ! HYBRID :
1392      DO jj = 1, jpj
1393         DO ji = 1, jpi
1394            DO jk = 1, jpkm1
1395               IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
1396            END DO
1397         END DO
1398      END DO
1399      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
1400         &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
1401
1402      IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain
1403         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) )
1404         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
1405            &                          ' w ', MINVAL( gdepw_0(:,:,:) )
1406         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0  (:,:,:) ),   &
1407            &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0  (:,:,:) ),   &
1408            &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0 (:,:,:) ),   &
1409            &                          ' w ', MINVAL( e3w_0  (:,:,:) )
1410
1411         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
1412            &                          ' w ', MAXVAL( gdepw_0(:,:,:) )
1413         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0  (:,:,:) ),   &
1414            &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0  (:,:,:) ),   &
1415            &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0 (:,:,:) ),   &
1416            &                          ' w ', MAXVAL( e3w_0  (:,:,:) )
1417      ENDIF
1418      !  END DO
1419      IF(lwp) THEN                                  ! selected vertical profiles
1420         WRITE(numout,*)
1421         WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
1422         WRITE(numout,*) ' ~~~~~~  --------------------'
1423         WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
1424         WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     &
1425            &                                 e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
1426         DO jj = mj0(20), mj1(20)
1427            DO ji = mi0(20), mi1(20)
1428               WRITE(numout,*)
1429               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
1430               WRITE(numout,*) ' ~~~~~~  --------------------'
1431               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
1432               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
1433                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
1434            END DO
1435         END DO
1436         DO jj = mj0(74), mj1(74)
1437            DO ji = mi0(100), mi1(100)
1438               WRITE(numout,*)
1439               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
1440               WRITE(numout,*) ' ~~~~~~  --------------------'
1441               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
1442               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
1443                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
1444            END DO
1445         END DO
1446      ENDIF
1447      !
1448      !================================================================================
1449      ! check the coordinate makes sense
1450      !================================================================================
1451      DO ji = 1, jpi
1452         DO jj = 1, jpj
1453            !
1454            IF( hbatt(ji,jj) > 0._wp) THEN
1455               DO jk = 1, mbathy(ji,jj)
1456                 ! check coordinate is monotonically increasing
1457                 IF (e3w_0(ji,jj,jk) <= 0._wp .OR. e3t_0(ji,jj,jk) <= 0._wp ) THEN
1458                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
1459                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
1460                    WRITE(numout,*) 'e3w',e3w_0(ji,jj,:)
1461                    WRITE(numout,*) 'e3t',e3t_0(ji,jj,:)
1462                    CALL ctl_stop( ctmp1 )
1463                 ENDIF
1464                 ! and check it has never gone negative
1465                 IF( gdepw_0(ji,jj,jk) < 0._wp .OR. gdept_0(ji,jj,jk) < 0._wp ) THEN
1466                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
1467                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk
1468                    WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:)
1469                    WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
1470                    CALL ctl_stop( ctmp1 )
1471                 ENDIF
1472                 ! and check it never exceeds the total depth
1473                 IF( gdepw_0(ji,jj,jk) > hbatt(ji,jj) ) THEN
1474                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
1475                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
1476                    WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:)
1477                    CALL ctl_stop( ctmp1 )
1478                 ENDIF
1479               END DO
1480               !
1481               DO jk = 1, mbathy(ji,jj)-1
1482                 ! and check it never exceeds the total depth
1483                IF( gdept_0(ji,jj,jk) > hbatt(ji,jj) ) THEN
1484                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
1485                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
1486                    WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
1487                    CALL ctl_stop( ctmp1 )
1488                 ENDIF
1489               END DO
1490            ENDIF
1491         END DO
1492      END DO
1493      !
1494      DEALLOCATE( zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
1495      !
1496   END SUBROUTINE zgr_sco
1497
1498
1499   SUBROUTINE s_sh94()
1500      !!----------------------------------------------------------------------
1501      !!                  ***  ROUTINE s_sh94  ***
1502      !!                     
1503      !! ** Purpose :   stretch the s-coordinate system
1504      !!
1505      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994
1506      !!                mixed S/sigma coordinate
1507      !!
1508      !! Reference : Song and Haidvogel 1994.
1509      !!----------------------------------------------------------------------
1510      INTEGER  ::   ji, jj, jk           ! dummy loop argument
1511      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
1512      REAL(wp) ::   ztmpu,  ztmpv,  ztmpf
1513      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1
1514      !
1515      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3
1516      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
1517      !!----------------------------------------------------------------------
1518
1519      ALLOCATE( z_gsigw3 (jpi,jpj,jpk), z_gsigt3 (jpi,jpj,jpk) )
1520      ALLOCATE( z_esigt3 (jpi,jpj,jpk), z_esigw3 (jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk), z_esigtv3(jpi,jpj,jpk) )
1521      ALLOCATE( z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk), z_esigwv3(jpi,jpj,jpk) )
1522
1523      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp
1524      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
1525      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
1526      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
1527      !
1528      DO ji = 1, jpi
1529         DO jj = 1, jpj
1530            !
1531            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
1532               DO jk = 1, jpk
1533                  z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
1534                  z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
1535               END DO
1536            ELSE ! shallow water, uniform sigma
1537               DO jk = 1, jpk
1538                  z_gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
1539                  z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
1540                  END DO
1541            ENDIF
1542            !
1543            DO jk = 1, jpkm1
1544               z_esigt3(ji,jj,jk  ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
1545               z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
1546            END DO
1547            z_esigw3(ji,jj,1  ) = 2._wp * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) )
1548            z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
1549            !
1550            DO jk = 1, jpk
1551               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1552               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1553               gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
1554               gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
1555            END DO
1556           !
1557         END DO   ! for all jj's
1558      END DO    ! for all ji's
1559
1560      DO ji = 1, jpim1
1561         DO jj = 1, jpjm1
1562            ! extended for Wetting/Drying case
1563            ztmpu  = hbatt(ji,jj)+hbatt(ji+1,jj)
1564            ztmpv  = hbatt(ji,jj)+hbatt(ji,jj+1)
1565            ztmpf  = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
1566            ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
1567            ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
1568            ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
1569                   & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
1570            DO jk = 1, jpk
1571                 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
1572                        &              / ztmpu
1573                 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
1574                        &              / ztmpu
1575
1576                 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
1577                        &              / ztmpv
1578                 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
1579                        &              / ztmpv
1580
1581                 z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                               &
1582                        &            +   hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                               &
1583                        &            +   hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                               &
1584                        &            +   hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
1585
1586               !
1587               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1588               e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1589               e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1590               e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1591               !
1592               e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1593               e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1594               e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1595            END DO
1596        END DO
1597      END DO
1598      !
1599      DEALLOCATE( z_gsigw3, z_gsigt3                                                        )
1600      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1601      !
1602   END SUBROUTINE s_sh94
1603
1604
1605   SUBROUTINE s_sf12
1606      !!----------------------------------------------------------------------
1607      !!                  ***  ROUTINE s_sf12 ***
1608      !!                     
1609      !! ** Purpose :   stretch the s-coordinate system
1610      !!
1611      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
1612      !!                mixed S/sigma/Z coordinate
1613      !!
1614      !!                This method allows the maintenance of fixed surface and or
1615      !!                bottom cell resolutions (cf. geopotential coordinates)
1616      !!                within an analytically derived stretched S-coordinate framework.
1617      !!
1618      !!
1619      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
1620      !!----------------------------------------------------------------------
1621      INTEGER  ::   ji, jj, jk           ! dummy loop argument
1622      REAL(wp) ::   zsmth               ! smoothing around critical depth
1623      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
1624      REAL(wp) ::   ztmpu, ztmpv, ztmpf
1625      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1
1626      !
1627      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3
1628      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
1629      !!----------------------------------------------------------------------
1630      !
1631      ALLOCATE( z_gsigw3 (jpi,jpj,jpk), z_gsigt3 (jpi,jpj,jpk) )
1632      ALLOCATE( z_esigt3 (jpi,jpj,jpk), z_esigw3 (jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk), z_esigtv3(jpi,jpj,jpk))
1633      ALLOCATE( z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk), z_esigwv3(jpi,jpj,jpk) )
1634
1635      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp
1636      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
1637      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
1638      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
1639
1640      DO ji = 1, jpi
1641         DO jj = 1, jpj
1642
1643          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
1644             
1645              zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
1646                                                     ! could be changed by users but care must be taken to do so carefully
1647              zzb = 1.0_wp-(zzb/hbatt(ji,jj))
1648           
1649              zzs = rn_zs / hbatt(ji,jj) 
1650             
1651              IF (rn_efold /= 0.0_wp) THEN
1652                zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
1653              ELSE
1654                zsmth = 1.0_wp 
1655              ENDIF
1656               
1657              DO jk = 1, jpk
1658                z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
1659                z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
1660              ENDDO
1661              z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth  )
1662              z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth  )
1663 
1664          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
1665
1666            DO jk = 1, jpk
1667              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
1668              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
1669            END DO
1670
1671          ELSE  ! shallow water, z coordinates
1672
1673            DO jk = 1, jpk
1674              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
1675              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
1676            END DO
1677
1678          ENDIF
1679
1680          DO jk = 1, jpkm1
1681             z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
1682             z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
1683          END DO
1684          z_esigw3(ji,jj,1  ) = 2.0_wp * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ))
1685          z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
1686
1687          DO jk = 1, jpk
1688             gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
1689             gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
1690          END DO
1691
1692        ENDDO   ! for all jj's
1693      ENDDO    ! for all ji's
1694
1695      DO ji=1,jpi-1
1696        DO jj=1,jpj-1
1697
1698           ! extend to suit for Wetting/Drying case
1699           ztmpu  = hbatt(ji,jj)+hbatt(ji+1,jj)
1700           ztmpv  = hbatt(ji,jj)+hbatt(ji,jj+1)
1701           ztmpf  = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
1702           ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
1703           ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
1704           ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
1705                  & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
1706           DO jk = 1, jpk
1707                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
1708                       &              / ztmpu
1709                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
1710                       &              / ztmpu
1711                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
1712                       &              / ztmpv
1713                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
1714                       &              / ztmpv
1715
1716                z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                               &
1717                       &              + hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                               &
1718                       &              + hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                               &
1719                       &              + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
1720
1721             ! Code prior to wetting and drying option (for reference)
1722             !z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
1723             !                     /( hbatt(ji,jj)+hbatt(ji+1,jj) )
1724             !
1725             !z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
1726             !                     /( hbatt(ji,jj)+hbatt(ji+1,jj) )
1727             !
1728             !z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
1729             !                     /( hbatt(ji,jj)+hbatt(ji,jj+1) )
1730             !
1731             !z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
1732             !                     /( hbatt(ji,jj)+hbatt(ji,jj+1) )
1733             !
1734             !z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                                 &
1735             !                    &  +hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                                 &
1736             !                       +hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                                 &
1737             !                    &  +hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )                               &
1738             !                     /( hbatt(ji  ,jj  )+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
1739
1740             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
1741             e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
1742             e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
1743             e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
1744             !
1745             e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
1746             e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
1747             e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
1748          END DO
1749
1750        ENDDO
1751      ENDDO
1752      !
1753      CALL lbc_lnk('domzgr',e3t_0 ,'T',1.) ; CALL lbc_lnk('domzgr',e3u_0 ,'T',1.)
1754      CALL lbc_lnk('domzgr',e3v_0 ,'T',1.) ; CALL lbc_lnk('domzgr',e3f_0 ,'T',1.)
1755      CALL lbc_lnk('domzgr',e3w_0 ,'T',1.)
1756      CALL lbc_lnk('domzgr',e3uw_0,'T',1.) ; CALL lbc_lnk('domzgr',e3vw_0,'T',1.)
1757      !
1758      DEALLOCATE( z_gsigw3, z_gsigt3                                                        )
1759      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1760      !
1761   END SUBROUTINE s_sf12
1762
1763
1764   SUBROUTINE s_tanh()
1765      !!----------------------------------------------------------------------
1766      !!                  ***  ROUTINE s_tanh***
1767      !!                     
1768      !! ** Purpose :   stretch the s-coordinate system
1769      !!
1770      !! ** Method  :   s-coordinate stretch
1771      !!
1772      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
1773      !!----------------------------------------------------------------------
1774      INTEGER  ::   ji, jj, jk       ! dummy loop argument
1775      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
1776      REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_gsigw, z_gsigt
1777      REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_esigt, z_esigw
1778      !!----------------------------------------------------------------------
1779
1780      ALLOCATE( z_gsigw(jpk), z_gsigt(jpk) )
1781      ALLOCATE( z_esigt(jpk), z_esigw(jpk) )
1782
1783      z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp
1784      z_esigt  = 0._wp   ;   z_esigw  = 0._wp 
1785
1786      DO jk = 1, jpk
1787        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
1788        z_gsigt(jk) = -fssig( REAL(jk,wp)        )
1789      END DO
1790      IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk)
1791      !
1792      ! Coefficients for vertical scale factors at w-, t- levels
1793!!gm bug :  define it from analytical function, not like juste bellow....
1794!!gm        or betteroffer the 2 possibilities....
1795      DO jk = 1, jpkm1
1796         z_esigt(jk  ) = z_gsigw(jk+1) - z_gsigw(jk)
1797         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
1798      END DO
1799      z_esigw( 1 ) = 2._wp * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
1800      z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
1801      !
1802      DO jk = 1, jpk
1803         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1804         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1805         gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
1806         gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
1807      END DO
1808
1809      DO jj = 1, jpj
1810         DO ji = 1, jpi
1811            DO jk = 1, jpk
1812              e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1813              e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1814              e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1815              e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
1816              !
1817              e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1818              e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1819              e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1820            END DO
1821         END DO
1822      END DO
1823      !
1824      DEALLOCATE( z_gsigw, z_gsigt )
1825      DEALLOCATE( z_esigt, z_esigw )
1826      !
1827   END SUBROUTINE s_tanh
1828
1829
1830   FUNCTION fssig( pk ) RESULT( pf )
1831      !!----------------------------------------------------------------------
1832      !!                 ***  ROUTINE fssig ***
1833      !!       
1834      !! ** Purpose :   provide the analytical function in s-coordinate
1835      !!         
1836      !! ** Method  :   the function provide the non-dimensional position of
1837      !!                T and W (i.e. between 0 and 1)
1838      !!                T-points at integer values (between 1 and jpk)
1839      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1840      !!----------------------------------------------------------------------
1841      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
1842      REAL(wp)             ::   pf   ! sigma value
1843      !!----------------------------------------------------------------------
1844      !
1845      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   &
1846         &     - TANH( rn_thetb * rn_theta                                )  )   &
1847         & * (   COSH( rn_theta                           )                      &
1848         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              &
1849         & / ( 2._wp * SINH( rn_theta ) )
1850      !
1851   END FUNCTION fssig
1852
1853
1854   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
1855      !!----------------------------------------------------------------------
1856      !!                 ***  ROUTINE fssig1 ***
1857      !!
1858      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
1859      !!
1860      !! ** Method  :   the function provides the non-dimensional position of
1861      !!                T and W (i.e. between 0 and 1)
1862      !!                T-points at integer values (between 1 and jpk)
1863      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1864      !!----------------------------------------------------------------------
1865      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
1866      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
1867      REAL(wp)             ::   pf1   ! sigma value
1868      !!----------------------------------------------------------------------
1869      !
1870      IF ( rn_theta == 0 ) then      ! uniform sigma
1871         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
1872      ELSE                        ! stretched sigma
1873         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
1874            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
1875            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
1876      ENDIF
1877      !
1878   END FUNCTION fssig1
1879
1880
1881   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
1882      !!----------------------------------------------------------------------
1883      !!                 ***  ROUTINE fgamma  ***
1884      !!
1885      !! ** Purpose :   provide analytical function for the s-coordinate
1886      !!
1887      !! ** Method  :   the function provides the non-dimensional position of
1888      !!                T and W (i.e. between 0 and 1)
1889      !!                T-points at integer values (between 1 and jpk)
1890      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1891      !!
1892      !!                This method allows the maintenance of fixed surface and or
1893      !!                bottom cell resolutions (cf. geopotential coordinates)
1894      !!                within an analytically derived stretched S-coordinate framework.
1895      !!
1896      !! Reference  :   Siddorn and Furner, in prep
1897      !!----------------------------------------------------------------------
1898      REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate
1899      REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate
1900      REAL(wp), INTENT(in   ) ::   pzb            ! Bottom box depth
1901      REAL(wp), INTENT(in   ) ::   pzs            ! surface box depth
1902      REAL(wp), INTENT(in   ) ::   psmth          ! Smoothing parameter
1903      !
1904      INTEGER  ::   jk             ! dummy loop index
1905      REAL(wp) ::   za1,za2,za3    ! local scalar
1906      REAL(wp) ::   zn1,zn2        !   -      -
1907      REAL(wp) ::   za,zb,zx       !   -      -
1908      !!----------------------------------------------------------------------
1909      !
1910      zn1  =  1._wp / REAL( jpkm1, wp )
1911      zn2  =  1._wp -  zn1
1912      !
1913      za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 
1914      za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp)
1915      za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
1916      !
1917      za = pzb - za3*(pzs-za1)-za2
1918      za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) )
1919      zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
1920      zx = 1.0_wp-za/2.0_wp-zb
1921      !
1922      DO jk = 1, jpk
1923         p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp +  &
1924            &          zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &
1925            &               (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
1926        p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth)
1927      END DO
1928      !
1929   END FUNCTION fgamma
1930
1931   !!======================================================================
1932END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.