New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domzgr.F90 in branches/2013/dev_MERCATOR_2013/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2013/dev_MERCATOR_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 4211

Last change on this file since 4211 was 4211, checked in by cbricaud, 10 years ago

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