source: utils/tools_ticket2457/DOMAINcfg/src/domzgr.F90 @ 12871

Last change on this file since 12871 was 12871, checked in by mathiot, 9 months ago

ticket #2457: add fix suggested in ticket.

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      ! bathymetry in level (from bathy_meter)
898      ! ===================
899      zmax = gdepw_1d(jpk) + e3t_1d(jpk)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) )
900      bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat)
901      WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0
902      ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level
903      END WHERE
904
905      ! Compute mbathy for ocean points (i.e. the number of ocean levels)
906      ! find the number of ocean levels such that the last level thickness
907      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where
908      ! e3t_1d is the reference level thickness
909      DO jk = jpkm1, 1, -1
910         zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
911         WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1
912      END DO
913
914      ! isf draft in z and level
915      ! ========================
916      ! if under ice shelf cavities, compute misfdep for ocean point (i.e. the top level)
917      IF ( ln_isfcav ) THEN
918         ! compute position of the ice shelf grounding line
919         CALL zgr_isf_zspace
920         ! compute misfdep and adjust isf draft if needed
921         CALL zgr_isf_kspace
922      END IF
923
924      ! Scale factors and depth at T- and W-points
925      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate
926         gdept_0(:,:,jk) = gdept_1d(jk)
927         gdepw_0(:,:,jk) = gdepw_1d(jk)
928         e3t_0  (:,:,jk) = e3t_1d  (jk)
929         e3w_0  (:,:,jk) = e3w_1d  (jk)
930      END DO
931     
932      ! Scale factors and depth at T- and W-points
933      DO jj = 1, jpj
934         DO ji = 1, jpi
935            ik = mbathy(ji,jj)
936            IF( ik > 0 ) THEN               ! ocean point only
937               ! max ocean level case
938               IF( ik == jpkm1 ) THEN
939                  zdepwp = bathy(ji,jj)
940                  ze3tp  = bathy(ji,jj) - gdepw_1d(ik)
941                  ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
942                  e3t_0(ji,jj,ik  ) = ze3tp
943                  e3t_0(ji,jj,ik+1) = ze3tp
944                  e3w_0(ji,jj,ik  ) = ze3wp
945                  e3w_0(ji,jj,ik+1) = ze3tp
946                  gdepw_0(ji,jj,ik+1) = zdepwp
947                  gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp
948                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
949                  !
950               ELSE                         ! standard case
951                  IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
952                  ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
953                  ENDIF
954!gm Bug?  check the gdepw_1d
955                  !       ... on ik
956                  gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   &
957                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   &
958                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) ))
959                  e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   & 
960                     &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ) 
961                  e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   &
962                     &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
963                  !       ... on ik+1
964                  e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
965                  e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
966                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
967               ENDIF
968            ENDIF
969         END DO
970      END DO
971      !
972      it = 0
973      DO jj = 1, jpj
974         DO ji = 1, jpi
975            ik = mbathy(ji,jj)
976            IF( ik > 0 ) THEN               ! ocean point only
977               e3tp (ji,jj) = e3t_0(ji,jj,ik)
978               e3wp (ji,jj) = e3w_0(ji,jj,ik)
979               ! test
980               zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  )
981               IF( zdiff <= 0._wp .AND. lwp ) THEN
982                  it = it + 1
983                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj
984                  WRITE(numout,*) ' bathy = ', bathy(ji,jj)
985                  WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
986                  WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  )
987               ENDIF
988            ENDIF
989         END DO
990      END DO
991      !
992      ! compute top scale factor if ice shelf
993      IF (ln_isfcav) CALL zps_isf
994      !
995      ! Scale factors and depth at U-, V-, UW and VW-points
996      DO jk = 1, jpk                        ! initialisation to z-scale factors
997         e3u_0 (:,:,jk) = e3t_1d(jk)
998         e3v_0 (:,:,jk) = e3t_1d(jk)
999         e3uw_0(:,:,jk) = e3w_1d(jk)
1000         e3vw_0(:,:,jk) = e3w_1d(jk)
1001      END DO
1002
1003      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
1004         DO jj = 1, jpjm1
1005            DO ji = 1, jpim1   ! vector opt.
1006               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
1007               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
1008               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
1009               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
1010            END DO
1011         END DO
1012      END DO
1013
1014      ! update e3uw in case only 2 cells in the water column
1015      IF ( ln_isfcav ) CALL zps_isf_e3uv_w
1016      !
1017      CALL lbc_lnk('domzgr', e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk('domzgr', e3uw_0, 'U', 1._wp )   ! lateral boundary conditions
1018      CALL lbc_lnk('domzgr', e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp )
1019      !
1020      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
1021         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk)
1022         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk)
1023         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk)
1024         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk)
1025      END DO
1026     
1027      ! Scale factor at F-point
1028      DO jk = 1, jpk                        ! initialisation to z-scale factors
1029         e3f_0(:,:,jk) = e3t_1d(jk)
1030      END DO
1031      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors
1032         DO jj = 1, jpjm1
1033            DO ji = 1, jpim1   ! vector opt.
1034               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
1035            END DO
1036         END DO
1037      END DO
1038      CALL lbc_lnk('domzgr', e3f_0, 'F', 1._wp )       ! Lateral boundary conditions
1039      !
1040      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
1041         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk)
1042      END DO
1043!!gm  bug ? :  must be a do loop with mj0,mj1
1044      !
1045      e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2
1046      e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 
1047      e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 
1048      e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 
1049      e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 
1050
1051      ! Control of the sign
1052      IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' )
1053      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' )
1054      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' )
1055      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' )
1056      !
1057      ! if in the future gde3w_0 need to be compute, use the function defined in NEMO
1058      ! for now gde3w_0 computation is removed as not an output of domcfg
1059
1060      DEALLOCATE( zprt )
1061      !
1062   END SUBROUTINE zgr_zps
1063
1064   SUBROUTINE zgr_sco
1065      !!----------------------------------------------------------------------
1066      !!                  ***  ROUTINE zgr_sco  ***
1067      !!                     
1068      !! ** Purpose :   define the s-coordinate system
1069      !!
1070      !! ** Method  :   s-coordinate
1071      !!         The depth of model levels is defined as the product of an
1072      !!      analytical function by the local bathymetry, while the vertical
1073      !!      scale factors are defined as the product of the first derivative
1074      !!      of the analytical function by the bathymetry.
1075      !!      (this solution save memory as depth and scale factors are not
1076      !!      3d fields)
1077      !!          - Read bathymetry (in meters) at t-point and compute the
1078      !!         bathymetry at u-, v-, and f-points.
1079      !!            hbatu = mi( hbatt )
1080      !!            hbatv = mj( hbatt )
1081      !!            hbatf = mi( mj( hbatt ) )
1082      !!          - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
1083      !!         function and its derivative given as function.
1084      !!            z_gsigt(k) = fssig (k    )
1085      !!            z_gsigw(k) = fssig (k-0.5)
1086      !!            z_esigt(k) = fsdsig(k    )
1087      !!            z_esigw(k) = fsdsig(k-0.5)
1088      !!      Three options for stretching are give, and they can be modified
1089      !!      following the users requirements. Nevertheless, the output as
1090      !!      well as the way to compute the model levels and scale factors
1091      !!      must be respected in order to insure second order accuracy
1092      !!      schemes.
1093      !!
1094      !!      The three methods for stretching available are:
1095      !!
1096      !!           s_sh94 (Song and Haidvogel 1994)
1097      !!                a sinh/tanh function that allows sigma and stretched sigma
1098      !!
1099      !!           s_sf12 (Siddorn and Furner 2012?)
1100      !!                allows the maintenance of fixed surface and or
1101      !!                bottom cell resolutions (cf. geopotential coordinates)
1102      !!                within an analytically derived stretched S-coordinate framework.
1103      !!
1104      !!          s_tanh  (Madec et al 1996)
1105      !!                a cosh/tanh function that gives stretched coordinates       
1106      !!
1107      !!----------------------------------------------------------------------
1108      INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
1109      INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
1110      INTEGER  ::   ios                      ! Local integer output status for namelist read
1111      REAL(wp) ::   zrmax, ztaper   ! temporary scalars
1112      REAL(wp) ::   zrfact
1113      !
1114      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
1115      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
1116      !!
1117      NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
1118         &                rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
1119     !!----------------------------------------------------------------------
1120      !
1121      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) )
1122      !
1123      REWIND( numnam_ref )              ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
1124      READ  ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901)
1125901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist', lwp )
1126
1127      REWIND( numnam_cfg )              ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters
1128      READ  ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 )
1129902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist', lwp )
1130      IF(lwm) WRITE ( numond, namzgr_sco )
1131
1132      IF(lwp) THEN                           ! control print
1133         WRITE(numout,*)
1134         WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
1135         WRITE(numout,*) '~~~~~~~~~~~'
1136         WRITE(numout,*) '   Namelist namzgr_sco'
1137         WRITE(numout,*) '     stretching coeffs '
1138         WRITE(numout,*) '        maximum depth of s-bottom surface (>0)       rn_sbot_max   = ',rn_sbot_max
1139         WRITE(numout,*) '        minimum depth of s-bottom surface (>0)       rn_sbot_min   = ',rn_sbot_min
1140         WRITE(numout,*) '        Critical depth                               rn_hc         = ',rn_hc
1141         WRITE(numout,*) '        maximum cut-off r-value allowed              rn_rmax       = ',rn_rmax
1142         WRITE(numout,*) '     Song and Haidvogel 1994 stretching              ln_s_sh94     = ',ln_s_sh94
1143         WRITE(numout,*) '        Song and Haidvogel 1994 stretching coefficients'
1144         WRITE(numout,*) '        surface control parameter (0<=rn_theta<=20)  rn_theta      = ',rn_theta
1145         WRITE(numout,*) '        bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ',rn_thetb
1146         WRITE(numout,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb
1147         WRITE(numout,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12
1148         WRITE(numout,*) '        switching to sigma (T) or Z (F) at H<Hc      ln_sigcrit    = ',ln_sigcrit
1149         WRITE(numout,*) '        Siddorn and Furner 2012 stretching coefficients'
1150         WRITE(numout,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha
1151         WRITE(numout,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold
1152         WRITE(numout,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs
1153         WRITE(numout,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a
1154         WRITE(numout,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b
1155         WRITE(numout,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b'
1156      ENDIF
1157
1158      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate
1159      hifu(:,:) = rn_sbot_min
1160      hifv(:,:) = rn_sbot_min
1161      hiff(:,:) = rn_sbot_min
1162
1163      !                                        ! set maximum ocean depth
1164      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1165
1166         DO jj = 1, jpj
1167            DO ji = 1, jpi
1168              IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1169            END DO
1170         END DO
1171      !                                        ! =============================
1172      !                                        ! Define the envelop bathymetry   (hbatt)
1173      !                                        ! =============================
1174      ! use r-value to create hybrid coordinates
1175      zenv(:,:) = bathy(:,:)
1176      !
1177      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
1178         DO jj = 1, jpj
1179            DO ji = 1, jpi
1180               IF( bathy(ji,jj) == 0._wp ) THEN
1181                  iip1 = MIN( ji+1, jpi )
1182                  ijp1 = MIN( jj+1, jpj )
1183                  iim1 = MAX( ji-1, 1 )
1184                  ijm1 = MAX( jj-1, 1 )
1185!!gm BUG fix see ticket #1617
1186                  IF( ( + bathy(iim1,ijm1) + bathy(ji,ijp1) + bathy(iip1,ijp1)              &
1187                     &  + bathy(iim1,jj  )                  + bathy(iip1,jj  )              &
1188                     &  + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1)  ) > 0._wp ) &
1189                     &    zenv(ji,jj) = rn_sbot_min
1190!!gm
1191!!gm               IF( ( bathy(iip1,jj  ) + bathy(iim1,jj  ) + bathy(ji,ijp1  ) + bathy(ji,ijm1) +         &
1192!!gm                  &  bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
1193!!gm               zenv(ji,jj) = rn_sbot_min
1194!!gm             ENDIF
1195!!gm end
1196              ENDIF
1197            END DO
1198         END DO
1199
1200      ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1201      CALL lbc_lnk( 'domzgr',zenv, 'T', 1._wp, 'no0' )
1202      !
1203      ! smooth the bathymetry (if required)
1204      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
1205      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
1206      !
1207      jl = 0
1208      zrmax = 1._wp
1209      !   
1210      !     
1211      ! set scaling factor used in reducing vertical gradients
1212      zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
1213      !
1214      ! initialise temporary evelope depth arrays
1215      ztmpi1(:,:) = zenv(:,:)
1216      ztmpi2(:,:) = zenv(:,:)
1217      ztmpj1(:,:) = zenv(:,:)
1218      ztmpj2(:,:) = zenv(:,:)
1219      !
1220      ! initialise temporary r-value arrays
1221      zri(:,:) = 1._wp
1222      zrj(:,:) = 1._wp
1223      !                                                            ! ================ !
1224      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) !  Iterative loop  !
1225         !                                                         ! ================ !
1226         jl = jl + 1
1227         zrmax = 0._wp
1228         ! we set zrmax from previous r-values (zri and zrj) first
1229         ! if set after current r-value calculation (as previously)
1230         ! we could exit DO WHILE prematurely before checking r-value
1231         ! of current zenv
1232         DO jj = 1, nlcj
1233            DO ji = 1, nlci
1234               zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
1235            END DO
1236         END DO
1237         zri(:,:) = 0._wp
1238         zrj(:,:) = 0._wp
1239         DO jj = 1, nlcj
1240            DO ji = 1, nlci
1241               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
1242               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
1243               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
1244                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
1245               END IF
1246               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
1247                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
1248               END IF
1249               IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact
1250               IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact
1251               IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact
1252               IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact
1253            END DO
1254         END DO
1255  !       IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
1256         !
1257         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax
1258         !
1259         DO jj = 1, nlcj
1260            DO ji = 1, nlci
1261               zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
1262            END DO
1263         END DO
1264         ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1265         CALL lbc_lnk( 'domzgr',zenv, 'T', 1._wp, 'no0' )
1266         !                                                  ! ================ !
1267      END DO                                                !     End loop     !
1268      !                                                     ! ================ !
1269      DO jj = 1, jpj
1270         DO ji = 1, jpi
1271            zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
1272         END DO
1273      END DO
1274      !
1275      ! Envelope bathymetry saved in hbatt
1276      hbatt(:,:) = zenv(:,:) 
1277      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
1278         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1279         DO jj = 1, jpj
1280            DO ji = 1, jpi
1281               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
1282               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
1283            END DO
1284         END DO
1285      ENDIF
1286      !
1287      !                                        ! ==============================
1288      !                                        !   hbatu, hbatv, hbatf fields
1289      !                                        ! ==============================
1290      IF(lwp) THEN
1291         WRITE(numout,*)
1292           WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
1293      ENDIF
1294      hbatu(:,:) = rn_sbot_min
1295      hbatv(:,:) = rn_sbot_min
1296      hbatf(:,:) = rn_sbot_min
1297      DO jj = 1, jpjm1
1298        DO ji = 1, jpim1   ! NO vector opt.
1299           hbatu(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
1300           hbatv(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
1301           hbatf(ji,jj) = 0.25_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
1302              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
1303        END DO
1304      END DO
1305
1306      !
1307      ! Apply lateral boundary condition
1308!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
1309      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk('domzgr', hbatu, 'U', 1._wp )
1310      DO jj = 1, jpj
1311         DO ji = 1, jpi
1312            IF( hbatu(ji,jj) == 0._wp ) THEN
1313               !No worries about the following line when ln_wd == .true.
1314               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
1315               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
1316            ENDIF
1317         END DO
1318      END DO
1319      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk('domzgr', hbatv, 'V', 1._wp )
1320      DO jj = 1, jpj
1321         DO ji = 1, jpi
1322            IF( hbatv(ji,jj) == 0._wp ) THEN
1323               IF( zhbat(ji,jj) == 0._wp )   hbatv(ji,jj) = rn_sbot_min
1324               IF( zhbat(ji,jj) /= 0._wp )   hbatv(ji,jj) = zhbat(ji,jj)
1325            ENDIF
1326         END DO
1327      END DO
1328      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk('domzgr', hbatf, 'F', 1._wp )
1329      DO jj = 1, jpj
1330         DO ji = 1, jpi
1331            IF( hbatf(ji,jj) == 0._wp ) THEN
1332               IF( zhbat(ji,jj) == 0._wp )   hbatf(ji,jj) = rn_sbot_min
1333               IF( zhbat(ji,jj) /= 0._wp )   hbatf(ji,jj) = zhbat(ji,jj)
1334            ENDIF
1335         END DO
1336      END DO
1337
1338!!bug:  key_helsinki a verifer
1339        hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1340        hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1341        hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1342        hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1343
1344      IF( nprint == 1 .AND. lwp )   THEN
1345         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
1346            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
1347         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
1348            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
1349         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
1350            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1351         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
1352            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1353      ENDIF
1354!! helsinki
1355
1356      !                                            ! =======================
1357      !                                            !   s-ccordinate fields     (gdep., e3.)
1358      !                                            ! =======================
1359      !
1360      ! non-dimensional "sigma" for model level depth at w- and t-levels
1361
1362
1363!========================================================================
1364! Song and Haidvogel  1994 (ln_s_sh94=T)
1365! Siddorn and Furner 2012 (ln_sf12=T)
1366! or  tanh function       (both false)                   
1367!========================================================================
1368      IF      ( ln_s_sh94 ) THEN
1369                           CALL s_sh94()
1370      ELSE IF ( ln_s_sf12 ) THEN
1371                           CALL s_sf12()
1372      ELSE                 
1373                           CALL s_tanh()
1374      ENDIF
1375
1376      CALL lbc_lnk( 'domzgr',e3t_0 , 'T', 1._wp )
1377      CALL lbc_lnk( 'domzgr',e3u_0 , 'U', 1._wp )
1378      CALL lbc_lnk( 'domzgr',e3v_0 , 'V', 1._wp )
1379      CALL lbc_lnk( 'domzgr',e3f_0 , 'F', 1._wp )
1380      CALL lbc_lnk( 'domzgr',e3w_0 , 'W', 1._wp )
1381      CALL lbc_lnk( 'domzgr',e3uw_0, 'U', 1._wp )
1382      CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp )
1383      !
1384        WHERE( e3t_0 (:,:,:) == 0._wp )   e3t_0 (:,:,:) = 1._wp
1385        WHERE( e3u_0 (:,:,:) == 0._wp )   e3u_0 (:,:,:) = 1._wp
1386        WHERE( e3v_0 (:,:,:) == 0._wp )   e3v_0 (:,:,:) = 1._wp
1387        WHERE( e3f_0 (:,:,:) == 0._wp )   e3f_0 (:,:,:) = 1._wp
1388        WHERE( e3w_0 (:,:,:) == 0._wp )   e3w_0 (:,:,:) = 1._wp
1389        WHERE( e3uw_0(:,:,:) == 0._wp )   e3uw_0(:,:,:) = 1._wp
1390        WHERE( e3vw_0(:,:,:) == 0._wp )   e3vw_0(:,:,:) = 1._wp
1391!!
1392      ! HYBRID :
1393      DO jj = 1, jpj
1394         DO ji = 1, jpi
1395            DO jk = 1, jpkm1
1396               IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
1397            END DO
1398         END DO
1399      END DO
1400      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
1401         &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
1402
1403      IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain
1404         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) )
1405         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
1406            &                          ' w ', MINVAL( gdepw_0(:,:,:) )
1407         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0  (:,:,:) ),   &
1408            &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0  (:,:,:) ),   &
1409            &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0 (:,:,:) ),   &
1410            &                          ' w ', MINVAL( e3w_0  (:,:,:) )
1411
1412         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
1413            &                          ' w ', MAXVAL( gdepw_0(:,:,:) )
1414         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0  (:,:,:) ),   &
1415            &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0  (:,:,:) ),   &
1416            &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0 (:,:,:) ),   &
1417            &                          ' w ', MAXVAL( e3w_0  (:,:,:) )
1418      ENDIF
1419      !  END DO
1420      IF(lwp) THEN                                  ! selected vertical profiles
1421         WRITE(numout,*)
1422         WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
1423         WRITE(numout,*) ' ~~~~~~  --------------------'
1424         WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
1425         WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     &
1426            &                                 e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
1427         DO jj = mj0(20), mj1(20)
1428            DO ji = mi0(20), mi1(20)
1429               WRITE(numout,*)
1430               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
1431               WRITE(numout,*) ' ~~~~~~  --------------------'
1432               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
1433               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
1434                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
1435            END DO
1436         END DO
1437         DO jj = mj0(74), mj1(74)
1438            DO ji = mi0(100), mi1(100)
1439               WRITE(numout,*)
1440               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
1441               WRITE(numout,*) ' ~~~~~~  --------------------'
1442               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
1443               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
1444                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
1445            END DO
1446         END DO
1447      ENDIF
1448      !
1449      !================================================================================
1450      ! check the coordinate makes sense
1451      !================================================================================
1452      DO ji = 1, jpi
1453         DO jj = 1, jpj
1454            !
1455            IF( hbatt(ji,jj) > 0._wp) THEN
1456               DO jk = 1, mbathy(ji,jj)
1457                 ! check coordinate is monotonically increasing
1458                 IF (e3w_0(ji,jj,jk) <= 0._wp .OR. e3t_0(ji,jj,jk) <= 0._wp ) THEN
1459                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
1460                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
1461                    WRITE(numout,*) 'e3w',e3w_0(ji,jj,:)
1462                    WRITE(numout,*) 'e3t',e3t_0(ji,jj,:)
1463                    CALL ctl_stop( ctmp1 )
1464                 ENDIF
1465                 ! and check it has never gone negative
1466                 IF( gdepw_0(ji,jj,jk) < 0._wp .OR. gdept_0(ji,jj,jk) < 0._wp ) THEN
1467                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
1468                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk
1469                    WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:)
1470                    WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
1471                    CALL ctl_stop( ctmp1 )
1472                 ENDIF
1473                 ! and check it never exceeds the total depth
1474                 IF( gdepw_0(ji,jj,jk) > hbatt(ji,jj) ) THEN
1475                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
1476                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
1477                    WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:)
1478                    CALL ctl_stop( ctmp1 )
1479                 ENDIF
1480               END DO
1481               !
1482               DO jk = 1, mbathy(ji,jj)-1
1483                 ! and check it never exceeds the total depth
1484                IF( gdept_0(ji,jj,jk) > hbatt(ji,jj) ) THEN
1485                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
1486                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
1487                    WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
1488                    CALL ctl_stop( ctmp1 )
1489                 ENDIF
1490               END DO
1491            ENDIF
1492         END DO
1493      END DO
1494      !
1495      DEALLOCATE( zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
1496      !
1497   END SUBROUTINE zgr_sco
1498
1499
1500   SUBROUTINE s_sh94()
1501      !!----------------------------------------------------------------------
1502      !!                  ***  ROUTINE s_sh94  ***
1503      !!                     
1504      !! ** Purpose :   stretch the s-coordinate system
1505      !!
1506      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994
1507      !!                mixed S/sigma coordinate
1508      !!
1509      !! Reference : Song and Haidvogel 1994.
1510      !!----------------------------------------------------------------------
1511      INTEGER  ::   ji, jj, jk           ! dummy loop argument
1512      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
1513      REAL(wp) ::   ztmpu,  ztmpv,  ztmpf
1514      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1
1515      !
1516      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3
1517      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
1518      !!----------------------------------------------------------------------
1519
1520      ALLOCATE( z_gsigw3 (jpi,jpj,jpk), z_gsigt3 (jpi,jpj,jpk) )
1521      ALLOCATE( z_esigt3 (jpi,jpj,jpk), z_esigw3 (jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk), z_esigtv3(jpi,jpj,jpk) )
1522      ALLOCATE( z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk), z_esigwv3(jpi,jpj,jpk) )
1523
1524      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp
1525      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
1526      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
1527      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
1528      !
1529      DO ji = 1, jpi
1530         DO jj = 1, jpj
1531            !
1532            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
1533               DO jk = 1, jpk
1534                  z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
1535                  z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
1536               END DO
1537            ELSE ! shallow water, uniform sigma
1538               DO jk = 1, jpk
1539                  z_gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
1540                  z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
1541                  END DO
1542            ENDIF
1543            !
1544            DO jk = 1, jpkm1
1545               z_esigt3(ji,jj,jk  ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
1546               z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
1547            END DO
1548            z_esigw3(ji,jj,1  ) = 2._wp * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) )
1549            z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
1550            !
1551            DO jk = 1, jpk
1552               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1553               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1554               gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
1555               gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
1556            END DO
1557           !
1558         END DO   ! for all jj's
1559      END DO    ! for all ji's
1560
1561      DO ji = 1, jpim1
1562         DO jj = 1, jpjm1
1563            ! extended for Wetting/Drying case
1564            ztmpu  = hbatt(ji,jj)+hbatt(ji+1,jj)
1565            ztmpv  = hbatt(ji,jj)+hbatt(ji,jj+1)
1566            ztmpf  = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
1567            ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
1568            ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
1569            ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
1570                   & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
1571            DO jk = 1, jpk
1572                 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
1573                        &              / ztmpu
1574                 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
1575                        &              / ztmpu
1576
1577                 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
1578                        &              / ztmpv
1579                 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
1580                        &              / ztmpv
1581
1582                 z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                               &
1583                        &            +   hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                               &
1584                        &            +   hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                               &
1585                        &            +   hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
1586
1587               !
1588               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1589               e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1590               e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1591               e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1592               !
1593               e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1594               e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1595               e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1596            END DO
1597        END DO
1598      END DO
1599      !
1600      DEALLOCATE( z_gsigw3, z_gsigt3                                                        )
1601      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1602      !
1603   END SUBROUTINE s_sh94
1604
1605
1606   SUBROUTINE s_sf12
1607      !!----------------------------------------------------------------------
1608      !!                  ***  ROUTINE s_sf12 ***
1609      !!                     
1610      !! ** Purpose :   stretch the s-coordinate system
1611      !!
1612      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
1613      !!                mixed S/sigma/Z coordinate
1614      !!
1615      !!                This method allows the maintenance of fixed surface and or
1616      !!                bottom cell resolutions (cf. geopotential coordinates)
1617      !!                within an analytically derived stretched S-coordinate framework.
1618      !!
1619      !!
1620      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
1621      !!----------------------------------------------------------------------
1622      INTEGER  ::   ji, jj, jk           ! dummy loop argument
1623      REAL(wp) ::   zsmth               ! smoothing around critical depth
1624      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
1625      REAL(wp) ::   ztmpu, ztmpv, ztmpf
1626      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1
1627      !
1628      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3
1629      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
1630      !!----------------------------------------------------------------------
1631      !
1632      ALLOCATE( z_gsigw3 (jpi,jpj,jpk), z_gsigt3 (jpi,jpj,jpk) )
1633      ALLOCATE( z_esigt3 (jpi,jpj,jpk), z_esigw3 (jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk), z_esigtv3(jpi,jpj,jpk))
1634      ALLOCATE( z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk), z_esigwv3(jpi,jpj,jpk) )
1635
1636      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp
1637      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
1638      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
1639      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
1640
1641      DO ji = 1, jpi
1642         DO jj = 1, jpj
1643
1644          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
1645             
1646              zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
1647                                                     ! could be changed by users but care must be taken to do so carefully
1648              zzb = 1.0_wp-(zzb/hbatt(ji,jj))
1649           
1650              zzs = rn_zs / hbatt(ji,jj) 
1651             
1652              IF (rn_efold /= 0.0_wp) THEN
1653                zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
1654              ELSE
1655                zsmth = 1.0_wp 
1656              ENDIF
1657               
1658              DO jk = 1, jpk
1659                z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
1660                z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
1661              ENDDO
1662              z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth  )
1663              z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth  )
1664 
1665          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
1666
1667            DO jk = 1, jpk
1668              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
1669              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
1670            END DO
1671
1672          ELSE  ! shallow water, z coordinates
1673
1674            DO jk = 1, jpk
1675              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
1676              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
1677            END DO
1678
1679          ENDIF
1680
1681          DO jk = 1, jpkm1
1682             z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
1683             z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
1684          END DO
1685          z_esigw3(ji,jj,1  ) = 2.0_wp * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ))
1686          z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
1687
1688          DO jk = 1, jpk
1689             gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
1690             gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
1691          END DO
1692
1693        ENDDO   ! for all jj's
1694      ENDDO    ! for all ji's
1695
1696      DO ji=1,jpi-1
1697        DO jj=1,jpj-1
1698
1699           ! extend to suit for Wetting/Drying case
1700           ztmpu  = hbatt(ji,jj)+hbatt(ji+1,jj)
1701           ztmpv  = hbatt(ji,jj)+hbatt(ji,jj+1)
1702           ztmpf  = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
1703           ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
1704           ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
1705           ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
1706                  & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
1707           DO jk = 1, jpk
1708                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
1709                       &              / ztmpu
1710                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
1711                       &              / ztmpu
1712                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
1713                       &              / ztmpv
1714                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
1715                       &              / ztmpv
1716
1717                z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                               &
1718                       &              + hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                               &
1719                       &              + hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                               &
1720                       &              + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
1721
1722             ! Code prior to wetting and drying option (for reference)
1723             !z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
1724             !                     /( hbatt(ji,jj)+hbatt(ji+1,jj) )
1725             !
1726             !z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
1727             !                     /( hbatt(ji,jj)+hbatt(ji+1,jj) )
1728             !
1729             !z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
1730             !                     /( hbatt(ji,jj)+hbatt(ji,jj+1) )
1731             !
1732             !z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
1733             !                     /( hbatt(ji,jj)+hbatt(ji,jj+1) )
1734             !
1735             !z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                                 &
1736             !                    &  +hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                                 &
1737             !                       +hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                                 &
1738             !                    &  +hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )                               &
1739             !                     /( hbatt(ji  ,jj  )+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
1740
1741             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
1742             e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
1743             e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
1744             e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
1745             !
1746             e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
1747             e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
1748             e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
1749          END DO
1750
1751        ENDDO
1752      ENDDO
1753      !
1754      CALL lbc_lnk('domzgr',e3t_0 ,'T',1.) ; CALL lbc_lnk('domzgr',e3u_0 ,'T',1.)
1755      CALL lbc_lnk('domzgr',e3v_0 ,'T',1.) ; CALL lbc_lnk('domzgr',e3f_0 ,'T',1.)
1756      CALL lbc_lnk('domzgr',e3w_0 ,'T',1.)
1757      CALL lbc_lnk('domzgr',e3uw_0,'T',1.) ; CALL lbc_lnk('domzgr',e3vw_0,'T',1.)
1758      !
1759      DEALLOCATE( z_gsigw3, z_gsigt3                                                        )
1760      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1761      !
1762   END SUBROUTINE s_sf12
1763
1764
1765   SUBROUTINE s_tanh()
1766      !!----------------------------------------------------------------------
1767      !!                  ***  ROUTINE s_tanh***
1768      !!                     
1769      !! ** Purpose :   stretch the s-coordinate system
1770      !!
1771      !! ** Method  :   s-coordinate stretch
1772      !!
1773      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
1774      !!----------------------------------------------------------------------
1775      INTEGER  ::   ji, jj, jk       ! dummy loop argument
1776      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
1777      REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_gsigw, z_gsigt
1778      REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_esigt, z_esigw
1779      !!----------------------------------------------------------------------
1780
1781      ALLOCATE( z_gsigw(jpk), z_gsigt(jpk) )
1782      ALLOCATE( z_esigt(jpk), z_esigw(jpk) )
1783
1784      z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp
1785      z_esigt  = 0._wp   ;   z_esigw  = 0._wp 
1786
1787      DO jk = 1, jpk
1788        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
1789        z_gsigt(jk) = -fssig( REAL(jk,wp)        )
1790      END DO
1791      IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk)
1792      !
1793      ! Coefficients for vertical scale factors at w-, t- levels
1794!!gm bug :  define it from analytical function, not like juste bellow....
1795!!gm        or betteroffer the 2 possibilities....
1796      DO jk = 1, jpkm1
1797         z_esigt(jk  ) = z_gsigw(jk+1) - z_gsigw(jk)
1798         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
1799      END DO
1800      z_esigw( 1 ) = 2._wp * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
1801      z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
1802      !
1803      DO jk = 1, jpk
1804         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1805         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1806         gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
1807         gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
1808      END DO
1809
1810      DO jj = 1, jpj
1811         DO ji = 1, jpi
1812            DO jk = 1, jpk
1813              e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1814              e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1815              e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1816              e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
1817              !
1818              e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1819              e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1820              e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1821            END DO
1822         END DO
1823      END DO
1824      !
1825      DEALLOCATE( z_gsigw, z_gsigt )
1826      DEALLOCATE( z_esigt, z_esigw )
1827      !
1828   END SUBROUTINE s_tanh
1829
1830
1831   FUNCTION fssig( pk ) RESULT( pf )
1832      !!----------------------------------------------------------------------
1833      !!                 ***  ROUTINE fssig ***
1834      !!       
1835      !! ** Purpose :   provide the analytical function in s-coordinate
1836      !!         
1837      !! ** Method  :   the function provide the non-dimensional position of
1838      !!                T and W (i.e. between 0 and 1)
1839      !!                T-points at integer values (between 1 and jpk)
1840      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1841      !!----------------------------------------------------------------------
1842      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
1843      REAL(wp)             ::   pf   ! sigma value
1844      !!----------------------------------------------------------------------
1845      !
1846      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   &
1847         &     - TANH( rn_thetb * rn_theta                                )  )   &
1848         & * (   COSH( rn_theta                           )                      &
1849         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              &
1850         & / ( 2._wp * SINH( rn_theta ) )
1851      !
1852   END FUNCTION fssig
1853
1854
1855   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
1856      !!----------------------------------------------------------------------
1857      !!                 ***  ROUTINE fssig1 ***
1858      !!
1859      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
1860      !!
1861      !! ** Method  :   the function provides the non-dimensional position of
1862      !!                T and W (i.e. between 0 and 1)
1863      !!                T-points at integer values (between 1 and jpk)
1864      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1865      !!----------------------------------------------------------------------
1866      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
1867      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
1868      REAL(wp)             ::   pf1   ! sigma value
1869      !!----------------------------------------------------------------------
1870      !
1871      IF ( rn_theta == 0 ) then      ! uniform sigma
1872         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
1873      ELSE                        ! stretched sigma
1874         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
1875            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
1876            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
1877      ENDIF
1878      !
1879   END FUNCTION fssig1
1880
1881
1882   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
1883      !!----------------------------------------------------------------------
1884      !!                 ***  ROUTINE fgamma  ***
1885      !!
1886      !! ** Purpose :   provide analytical function for the s-coordinate
1887      !!
1888      !! ** Method  :   the function provides the non-dimensional position of
1889      !!                T and W (i.e. between 0 and 1)
1890      !!                T-points at integer values (between 1 and jpk)
1891      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1892      !!
1893      !!                This method allows the maintenance of fixed surface and or
1894      !!                bottom cell resolutions (cf. geopotential coordinates)
1895      !!                within an analytically derived stretched S-coordinate framework.
1896      !!
1897      !! Reference  :   Siddorn and Furner, in prep
1898      !!----------------------------------------------------------------------
1899      REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate
1900      REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate
1901      REAL(wp), INTENT(in   ) ::   pzb            ! Bottom box depth
1902      REAL(wp), INTENT(in   ) ::   pzs            ! surface box depth
1903      REAL(wp), INTENT(in   ) ::   psmth          ! Smoothing parameter
1904      !
1905      INTEGER  ::   jk             ! dummy loop index
1906      REAL(wp) ::   za1,za2,za3    ! local scalar
1907      REAL(wp) ::   zn1,zn2        !   -      -
1908      REAL(wp) ::   za,zb,zx       !   -      -
1909      !!----------------------------------------------------------------------
1910      !
1911      zn1  =  1._wp / REAL( jpkm1, wp )
1912      zn2  =  1._wp -  zn1
1913      !
1914      za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 
1915      za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp)
1916      za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
1917      !
1918      za = pzb - za3*(pzs-za1)-za2
1919      za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) )
1920      zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
1921      zx = 1.0_wp-za/2.0_wp-zb
1922      !
1923      DO jk = 1, jpk
1924         p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp +  &
1925            &          zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &
1926            &               (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
1927        p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth)
1928      END DO
1929      !
1930   END FUNCTION fgamma
1931
1932   !!======================================================================
1933END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.