source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 3998

Last change on this file since 3998 was 3998, checked in by jamesharle, 8 years ago

remove new evelope smoothing code and revert back to r3925 version of the code - see ticket 1119

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