source: branches/2013/dev_r3853_CNRS9_ConfSetting/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 3989

Last change on this file since 3989 was 3989, checked in by clevy, 8 years ago

Configuration setting/Step3 and doc, see ticket:#1074

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