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

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

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