source: branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 3684

Last change on this file since 3684 was 3684, checked in by cbricaud, 9 years ago

split a line which too long in domzgr.F90

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