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_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 4385

Last change on this file since 4385 was 4385, checked in by hliu, 10 years ago

1) Modified Wetting and Drying flux limiter part
2) I will apply all the modifications to v3.6beta as soon as everything runs smoothly at this version.

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