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

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

source: branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 3453

Last change on this file since 3453 was 3453, checked in by rfurner, 12 years ago

changes to allow new strecthed coordinate, including namelist for AMM12 std config; ticket#985
  • Property svn:keywords set to Id
File size: 94.3 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 adn Furner stretching functio
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        : sigma coordinate non-dimensional function
31   !!       fgamma       : Siddorn and Furner stretching function
32   !!       dfssig       : derivative of the sigma coordinate function    !!gm  (currently missing!)
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   LOGICAL  ::   ln_sigcrit  = .false.      ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
54   !
55   REAL(wp) ::   rn_sbot_min =  300._wp     ! minimum depth of s-bottom surface (>0) (m)
56   REAL(wp) ::   rn_sbot_max = 5250._wp     ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
57   REAL(wp) ::   rn_rmax     =    0.15_wp   ! maximum cut-off r-value allowed (0<rn_rmax<1)
58   REAL(wp) ::   rn_hc       =  150._wp     ! Critical depth for transition from sigma to stretched coordinates
59   ! Song and Haidvogel 1994 stretching parameters
60   REAL(wp) ::   rn_theta    =    6.00_wp   ! surface control parameter (0<=rn_theta<=20)
61   REAL(wp) ::   rn_thetb    =    0.75_wp   ! bottom control parameter  (0<=rn_thetb<= 1)
62   REAL(wp) ::   rn_bb       =    0.80_wp   ! stretching parameter
63   !                                        ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
64   ! Siddorn and Furner stretching parameters
65   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 gsigt, gsigw, esigt, esigw from an analytical
1068      !!         function and its derivative given as function.
1069      !!            gsigt(k) = fssig (k    )
1070      !!            gsigw(k) = fssig (k-0.5)
1071      !!            esigt(k) = fsdsig(k    )
1072      !!            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,*) '        Siddorn and Furner 2012 stretching coefficients'
1128         WRITE(numout,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha
1129         WRITE(numout,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold
1130         WRITE(numout,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs
1131         WRITE(numout,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a
1132         WRITE(numout,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b
1133         WRITE(numout,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b'
1134      ENDIF
1135
1136      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate
1137      hifu(:,:) = rn_sbot_min
1138      hifv(:,:) = rn_sbot_min
1139      hiff(:,:) = rn_sbot_min
1140
1141      !                                        ! set maximum ocean depth
1142      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1143
1144      DO jj = 1, jpj
1145         DO ji = 1, jpi
1146           IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1147         END DO
1148      END DO
1149      !                                        ! =============================
1150      !                                        ! Define the envelop bathymetry   (hbatt)
1151      !                                        ! =============================
1152      ! use r-value to create hybrid coordinates
1153      DO jj = 1, jpj
1154         DO ji = 1, jpi
1155            zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min )
1156         END DO
1157      END DO
1158      !
1159      ! Smooth the bathymetry (if required)
1160      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
1161      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
1162      !
1163      jl = 0
1164      zrmax = 1._wp
1165      !                                                     ! ================ !
1166      DO WHILE( jl <= 10000 .AND. zrmax > rn_rmax )         !  Iterative loop  !
1167         !                                                  ! ================ !
1168         jl = jl + 1
1169         zrmax = 0._wp
1170         zmsk(:,:) = 0._wp
1171         DO jj = 1, nlcj
1172            DO ji = 1, nlci
1173               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
1174               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
1175               zri(ji,jj) = ABS( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
1176               zrj(ji,jj) = ABS( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
1177               zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) )
1178               IF( zri(ji,jj) > rn_rmax )   zmsk(ji  ,jj  ) = 1._wp
1179               IF( zri(ji,jj) > rn_rmax )   zmsk(iip1,jj  ) = 1._wp
1180               IF( zrj(ji,jj) > rn_rmax )   zmsk(ji  ,jj  ) = 1._wp
1181               IF( zrj(ji,jj) > rn_rmax )   zmsk(ji  ,ijp1) = 1._wp
1182            END DO
1183         END DO
1184         IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
1185         ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)
1186         ztmp(:,:) = zmsk(:,:)   ;   CALL lbc_lnk( zmsk, 'T', 1._wp )
1187         DO jj = 1, nlcj
1188            DO ji = 1, nlci
1189                zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )
1190            END DO
1191         END DO
1192         !
1193         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) )
1194         !
1195         DO jj = 1, nlcj
1196            DO ji = 1, nlci
1197               iip1 = MIN( ji+1, nlci )     ! last  line (ji=nlci)
1198               ijp1 = MIN( jj+1, nlcj )     ! last  raw  (jj=nlcj)
1199               iim1 = MAX( ji-1,  1  )      ! first line (ji=nlci)
1200               ijm1 = MAX( jj-1,  1  )      ! first raw  (jj=nlcj)
1201               IF( zmsk(ji,jj) == 1._wp ) THEN
1202                  ztmp(ji,jj) =   (                                                                                   &
1203             &      zenv(iim1,ijp1)*zmsk(iim1,ijp1) + zenv(ji,ijp1)*zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1)   &
1204             &    + zenv(iim1,jj  )*zmsk(iim1,jj  ) + zenv(ji,jj  )*    2._wp     + zenv(iip1,jj  )*zmsk(iip1,jj  )   &
1205             &    + zenv(iim1,ijm1)*zmsk(iim1,ijm1) + zenv(ji,ijm1)*zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1)   &
1206             &                    ) / (                                                                               &
1207             &                      zmsk(iim1,ijp1) +               zmsk(ji,ijp1) +                 zmsk(iip1,ijp1)   &
1208             &    +                 zmsk(iim1,jj  ) +                   2._wp     +                 zmsk(iip1,jj  )   &
1209             &    +                 zmsk(iim1,ijm1) +               zmsk(ji,ijm1) +                 zmsk(iip1,ijm1)   &
1210             &                        )
1211               ENDIF
1212            END DO
1213         END DO
1214         !
1215         DO jj = 1, nlcj
1216            DO ji = 1, nlci
1217               IF( zmsk(ji,jj) == 1._wp )   zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )
1218            END DO
1219         END DO
1220         !
1221         ! Apply lateral boundary condition   CAUTION: kept the value when the lbc field is zero
1222         ztmp(:,:) = zenv(:,:)   ;   CALL lbc_lnk( zenv, 'T', 1._wp )
1223         DO jj = 1, nlcj
1224            DO ji = 1, nlci
1225               IF( zenv(ji,jj) == 0._wp )   zenv(ji,jj) = ztmp(ji,jj)
1226            END DO
1227         END DO
1228         !                                                  ! ================ !
1229      END DO                                                !     End loop     !
1230      !                                                     ! ================ !
1231      !
1232      !                                        ! envelop bathymetry saved in hbatt
1233      hbatt(:,:) = zenv(:,:) 
1234      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
1235         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1236         DO jj = 1, jpj
1237            DO ji = 1, jpi
1238               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2 )
1239               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
1240            END DO
1241         END DO
1242      ENDIF
1243      !
1244      IF(lwp) THEN                             ! Control print
1245         WRITE(numout,*)
1246         WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
1247         WRITE(numout,*)
1248         CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )
1249         IF( nprint == 1 )   THEN       
1250            WRITE(numout,*) ' bathy  MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
1251            WRITE(numout,*) ' hbatt  MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
1252         ENDIF
1253      ENDIF
1254
1255      !                                        ! ==============================
1256      !                                        !   hbatu, hbatv, hbatf fields
1257      !                                        ! ==============================
1258      IF(lwp) THEN
1259         WRITE(numout,*)
1260         WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
1261      ENDIF
1262      hbatu(:,:) = rn_sbot_min
1263      hbatv(:,:) = rn_sbot_min
1264      hbatf(:,:) = rn_sbot_min
1265      DO jj = 1, jpjm1
1266        DO ji = 1, jpim1   ! NO vector opt.
1267           hbatu(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
1268           hbatv(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
1269           hbatf(ji,jj) = 0.25_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
1270              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
1271        END DO
1272      END DO
1273      !
1274      ! Apply lateral boundary condition
1275!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
1276      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk( hbatu, 'U', 1._wp )
1277      DO jj = 1, jpj
1278         DO ji = 1, jpi
1279            IF( hbatu(ji,jj) == 0._wp ) THEN
1280               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
1281               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
1282            ENDIF
1283         END DO
1284      END DO
1285      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk( hbatv, 'V', 1._wp )
1286      DO jj = 1, jpj
1287         DO ji = 1, jpi
1288            IF( hbatv(ji,jj) == 0._wp ) THEN
1289               IF( zhbat(ji,jj) == 0._wp )   hbatv(ji,jj) = rn_sbot_min
1290               IF( zhbat(ji,jj) /= 0._wp )   hbatv(ji,jj) = zhbat(ji,jj)
1291            ENDIF
1292         END DO
1293      END DO
1294      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk( hbatf, 'F', 1._wp )
1295      DO jj = 1, jpj
1296         DO ji = 1, jpi
1297            IF( hbatf(ji,jj) == 0._wp ) THEN
1298               IF( zhbat(ji,jj) == 0._wp )   hbatf(ji,jj) = rn_sbot_min
1299               IF( zhbat(ji,jj) /= 0._wp )   hbatf(ji,jj) = zhbat(ji,jj)
1300            ENDIF
1301         END DO
1302      END DO
1303
1304!!bug:  key_helsinki a verifer
1305      hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1306      hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1307      hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1308      hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1309
1310      IF( nprint == 1 .AND. lwp )   THEN
1311         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
1312            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
1313         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
1314            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
1315         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
1316            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1317         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
1318            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1319      ENDIF
1320!! helsinki
1321
1322      !                                            ! =======================
1323      !                                            !   s-ccordinate fields     (gdep., e3.)
1324      !                                            ! =======================
1325      !
1326      ! non-dimensional "sigma" for model level depth at w- and t-levels
1327
1328
1329!========================================================================
1330! Song and Haidvogel  1994 (ln_s_sh94=T)
1331! Siddorn and Furner 2012 (ln_sf12=T)
1332! or  tanh function       (both false)                   
1333!========================================================================
1334      IF      ( ln_s_sh94 ) THEN
1335                           CALL s_sh94()
1336      ELSE IF ( ln_s_sf12 ) THEN
1337                           CALL s_sf12()
1338      ELSE                 
1339                           CALL s_tanh()
1340      ENDIF
1341
1342      CALL lbc_lnk( e3t , 'T', 1._wp )
1343      CALL lbc_lnk( e3u , 'U', 1._wp )
1344      CALL lbc_lnk( e3v , 'V', 1._wp )
1345      CALL lbc_lnk( e3f , 'F', 1._wp )
1346      CALL lbc_lnk( e3w , 'W', 1._wp )
1347      CALL lbc_lnk( e3uw, 'U', 1._wp )
1348      CALL lbc_lnk( e3vw, 'V', 1._wp )
1349
1350      fsdepw(:,:,:) = gdepw (:,:,:)
1351      fsde3w(:,:,:) = gdep3w(:,:,:)
1352      !
1353      where (e3t   (:,:,:).eq.0.0)  e3t(:,:,:) = 1.0
1354      where (e3u   (:,:,:).eq.0.0)  e3u(:,:,:) = 1.0
1355      where (e3v   (:,:,:).eq.0.0)  e3v(:,:,:) = 1.0
1356      where (e3f   (:,:,:).eq.0.0)  e3f(:,:,:) = 1.0
1357      where (e3w   (:,:,:).eq.0.0)  e3w(:,:,:) = 1.0
1358      where (e3uw  (:,:,:).eq.0.0)  e3uw(:,:,:) = 1.0
1359      where (e3vw  (:,:,:).eq.0.0)  e3vw(:,:,:) = 1.0
1360
1361
1362      fsdept(:,:,:) = gdept (:,:,:)
1363      fsdepw(:,:,:) = gdepw (:,:,:)
1364      fsde3w(:,:,:) = gdep3w(:,:,:)
1365      fse3t (:,:,:) = e3t   (:,:,:)
1366      fse3u (:,:,:) = e3u   (:,:,:)
1367      fse3v (:,:,:) = e3v   (:,:,:)
1368      fse3f (:,:,:) = e3f   (:,:,:)
1369      fse3w (:,:,:) = e3w   (:,:,:)
1370      fse3uw(:,:,:) = e3uw  (:,:,:)
1371      fse3vw(:,:,:) = e3vw  (:,:,:)
1372!!
1373      ! HYBRID :
1374      DO jj = 1, jpj
1375         DO ji = 1, jpi
1376            DO jk = 1, jpkm1
1377               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
1378               IF( scobot(ji,jj) == 0._wp            )   mbathy(ji,jj) = 0
1379            END DO
1380         END DO
1381      END DO
1382      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
1383         &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
1384
1385      !                                               ! =============
1386      IF(lwp) THEN                                    ! Control print
1387         !                                            ! =============
1388         WRITE(numout,*) 
1389         WRITE(numout,*) ' domzgr: vertical coefficients for model level'
1390         WRITE(numout, "(9x,'  level    gsigt      gsigw      esigt      esigw      gsi3w')" )
1391         WRITE(numout, "(10x,i4,5f11.4)" ) ( jk, gsigt(jk), gsigw(jk), esigt(jk), esigw(jk), gsi3w(jk), jk=1,jpk )
1392      ENDIF
1393      IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain
1394         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)   ), ' MAX ', MAXVAL( mbathy(:,:) )
1395         WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ),   &
1396            &                          ' w ', MINVAL( fsdepw(:,:,:) ), '3w '  , MINVAL( fsde3w(:,:,:) )
1397         WRITE(numout,*) ' MIN val e3    t ', MINVAL( fse3t (:,:,:) ), ' f '  , MINVAL( fse3f (:,:,:) ),   &
1398            &                          ' u ', MINVAL( fse3u (:,:,:) ), ' u '  , MINVAL( fse3v (:,:,:) ),   &
1399            &                          ' uw', MINVAL( fse3uw(:,:,:) ), ' vw'  , MINVAL( fse3vw(:,:,:) ),   &
1400            &                          ' w ', MINVAL( fse3w (:,:,:) )
1401
1402         WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ),   &
1403            &                          ' w ', MAXVAL( fsdepw(:,:,:) ), '3w '  , MAXVAL( fsde3w(:,:,:) )
1404         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( fse3t (:,:,:) ), ' f '  , MAXVAL( fse3f (:,:,:) ),   &
1405            &                          ' u ', MAXVAL( fse3u (:,:,:) ), ' u '  , MAXVAL( fse3v (:,:,:) ),   &
1406            &                          ' uw', MAXVAL( fse3uw(:,:,:) ), ' vw'  , MAXVAL( fse3vw(:,:,:) ),   &
1407            &                          ' w ', MAXVAL( fse3w (:,:,:) )
1408      ENDIF
1409      !  END DO
1410      IF(lwp) THEN                                  ! selected vertical profiles
1411         WRITE(numout,*)
1412         WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
1413         WRITE(numout,*) ' ~~~~~~  --------------------'
1414         WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
1415         WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk),     &
1416            &                                 fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk )
1417         DO jj = mj0(20), mj1(20)
1418            DO ji = mi0(20), mi1(20)
1419               WRITE(numout,*)
1420               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
1421               WRITE(numout,*) ' ~~~~~~  --------------------'
1422               WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
1423               WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk),     &
1424                  &                                 fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
1425            END DO
1426         END DO
1427         DO jj = mj0(74), mj1(74)
1428            DO ji = mi0(100), mi1(100)
1429               WRITE(numout,*)
1430               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
1431               WRITE(numout,*) ' ~~~~~~  --------------------'
1432               WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
1433               WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk),     &
1434                  &                                 fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
1435            END DO
1436         END DO
1437      ENDIF
1438
1439!================================================================================
1440! check the coordinate makes sense
1441!================================================================================
1442      DO ji = 1, jpi
1443         DO jj = 1, jpj
1444
1445            IF( hbatt(ji,jj) > 0._wp) THEN
1446               DO jk = 1, mbathy(ji,jj)
1447                 ! check coordinate is monotonically increasing
1448                 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
1449                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
1450                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
1451                    WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
1452                    WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
1453                    CALL ctl_stop( ctmp1 )
1454                 ENDIF
1455                 ! and check it has never gone negative
1456                 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
1457                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
1458                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk
1459                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
1460                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
1461                    CALL ctl_stop( ctmp1 )
1462                 ENDIF
1463                 ! and check it never exceeds the total depth
1464                 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
1465                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
1466                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
1467                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
1468                    CALL ctl_stop( ctmp1 )
1469                 ENDIF
1470               END DO
1471
1472               DO jk = 1, mbathy(ji,jj)-1
1473                 ! and check it never exceeds the total depth
1474                IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
1475                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
1476                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
1477                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
1478                    CALL ctl_stop( ctmp1 )
1479                 ENDIF
1480               END DO
1481
1482            ENDIF
1483
1484         END DO
1485      END DO
1486      !
1487      CALL wrk_dealloc( jpi, jpj,      zenv, ztmp, zmsk, zri, zrj, zhbat                           )
1488      !
1489      IF( nn_timing == 1 )  CALL timing_stop('zgr_sco')
1490      !
1491   END SUBROUTINE zgr_sco
1492
1493!!======================================================================
1494   SUBROUTINE s_sh94()
1495
1496      !!----------------------------------------------------------------------
1497      !!                  ***  ROUTINE s_sh94  ***
1498      !!                     
1499      !! ** Purpose :   stretch the s-coordinate system
1500      !!
1501      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994
1502      !!                mixed S/sigma coordinate
1503      !!
1504      !! Reference : Song and Haidvogel 1994.
1505      !!----------------------------------------------------------------------
1506      !
1507      INTEGER  ::   ji, jj, jk           ! dummy loop argument
1508      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
1509      !
1510      REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3
1511      REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3           
1512
1513      CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      )
1514      CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 )
1515
1516      gsigw3  = 0._wp   ;   gsigt3  = 0._wp   ;   gsi3w3  = 0._wp
1517      esigt3  = 0._wp   ;   esigw3  = 0._wp 
1518      esigtu3 = 0._wp   ;   esigtv3 = 0._wp   ;   esigtf3 = 0._wp
1519      esigwu3 = 0._wp   ;   esigwv3 = 0._wp
1520
1521      DO ji = 1, jpi
1522         DO jj = 1, jpj
1523
1524            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
1525               DO jk = 1, jpk
1526                  gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
1527                  gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
1528               END DO
1529            ELSE ! shallow water, uniform sigma
1530               DO jk = 1, jpk
1531                  gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
1532                  gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
1533                  END DO
1534            ENDIF
1535            !
1536            DO jk = 1, jpkm1
1537               esigt3(ji,jj,jk  ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk)
1538               esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk)
1539            END DO
1540            esigw3(ji,jj,1  ) = 2._wp * ( gsigt3(ji,jj,1  ) - gsigw3(ji,jj,1  ) )
1541            esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) )
1542            !
1543            ! Coefficients for vertical depth as the sum of e3w scale factors
1544            gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1)
1545            DO jk = 2, jpk
1546               gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk)
1547            END DO
1548            !
1549            DO jk = 1, jpk
1550               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1551               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1552               gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft )
1553               gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw )
1554               gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
1555            END DO
1556           !
1557         END DO   ! for all jj's
1558      END DO    ! for all ji's
1559
1560      DO ji = 1, jpim1
1561         DO jj = 1, jpjm1
1562            DO jk = 1, jpk
1563               esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) )   &
1564                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1565               esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) )   &
1566                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1567               esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk)     &
1568                  &                + hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) )   &
1569                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
1570               esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) )   &
1571                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1572               esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) )   &
1573                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1574               !
1575               e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1576               e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1577               e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1578               e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1579               !
1580               e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1581               e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1582               e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1583            END DO
1584        END DO
1585      END DO
1586
1587      CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      )
1588      CALL wrk_dealloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 )
1589
1590   END SUBROUTINE s_sh94
1591
1592   SUBROUTINE s_sf12
1593
1594      !!----------------------------------------------------------------------
1595      !!                  ***  ROUTINE s_sf12 ***
1596      !!                     
1597      !! ** Purpose :   stretch the s-coordinate system
1598      !!
1599      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
1600      !!                mixed S/sigma/Z coordinate
1601      !!
1602      !!                This method allows the maintenance of fixed surface and or
1603      !!                bottom cell resolutions (cf. geopotential coordinates)
1604      !!                within an analytically derived stretched S-coordinate framework.
1605      !!
1606      !!
1607      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
1608      !!----------------------------------------------------------------------
1609      !
1610      INTEGER  ::   ji, jj, jk           ! dummy loop argument
1611      REAL(wp) ::   fsmth          ! smoothing around critical depth
1612      REAL(wp) ::   zss, zbb       ! Surface and bottom cell thickness in sigma space
1613      !
1614      REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3
1615      REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3           
1616
1617      !
1618      CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      )
1619      CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 )
1620
1621      gsigw3  = 0._wp   ;   gsigt3  = 0._wp   ;   gsi3w3  = 0._wp
1622      esigt3  = 0._wp   ;   esigw3  = 0._wp 
1623      esigtu3 = 0._wp   ;   esigtv3 = 0._wp   ;   esigtf3 = 0._wp
1624      esigwu3 = 0._wp   ;   esigwv3 = 0._wp
1625
1626      DO ji = 1, jpi
1627         DO jj = 1, jpj
1628
1629          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
1630             
1631              zbb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
1632                                                     ! could be changed by users but care must be taken to do so carefully
1633              zbb = 1.0_wp-(zbb/hbatt(ji,jj))
1634           
1635              zss = rn_zs / hbatt(ji,jj) 
1636             
1637              IF (rn_efold /= 0.0_wp) THEN
1638                fsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
1639              ELSE
1640                fsmth = 1.0_wp 
1641              ENDIF
1642               
1643              DO jk = 1, jpk
1644                gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
1645                gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
1646              ENDDO
1647              gsigw3(ji,jj,:) = fgamma( gsigw3(ji,jj,:), zbb, zss, fsmth  )
1648              gsigt3(ji,jj,:) = fgamma( gsigt3(ji,jj,:), zbb, zss, fsmth  )
1649 
1650          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
1651
1652            DO jk = 1, jpk
1653              gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
1654              gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
1655            END DO
1656
1657          ELSE  ! shallow water, z coordinates
1658
1659            DO jk = 1, jpk
1660              gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
1661              gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
1662            END DO
1663
1664          ENDIF
1665
1666          DO jk = 1, jpkm1
1667             esigt3(ji,jj,jk) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk)
1668             esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk)
1669          END DO
1670          esigw3(ji,jj,1  ) = 2.0_wp * (gsigt3(ji,jj,1  ) - gsigw3(ji,jj,1  ))
1671          esigt3(ji,jj,jpk) = 2.0_wp * (gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk))
1672
1673          ! Coefficients for vertical depth as the sum of e3w scale factors
1674          gsi3w3(ji,jj,1) = 0.5 * esigw3(ji,jj,1)
1675          DO jk = 2, jpk
1676             gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk)
1677          END DO
1678
1679          DO jk = 1, jpk
1680             gdept (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*gsigt3(ji,jj,jk)
1681             gdepw (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*gsigw3(ji,jj,jk)
1682             gdep3w(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*gsi3w3(ji,jj,jk)
1683          END DO
1684
1685        ENDDO   ! for all jj's
1686      ENDDO    ! for all ji's
1687
1688      DO ji=1,jpi
1689        DO jj=1,jpj
1690
1691          DO jk = 1, jpk
1692                esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) ) / &
1693                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1694                esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) ) / &
1695                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1696                esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) +  &
1697                                      hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) ) / &
1698                                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
1699                esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) ) / &
1700                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1701                esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) ) / &
1702                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1703
1704             e3t(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*esigt3(ji,jj,jk)
1705             e3u(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*esigtu3(ji,jj,jk)
1706             e3v(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*esigtv3(ji,jj,jk)
1707             e3f(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*esigtf3(ji,jj,jk)
1708             !
1709             e3w(ji,jj,jk)=hbatt(ji,jj)*esigw3(ji,jj,jk)
1710             e3uw(ji,jj,jk)=hbatu(ji,jj)*esigwu3(ji,jj,jk)
1711             e3vw(ji,jj,jk)=hbatv(ji,jj)*esigwv3(ji,jj,jk)
1712          END DO
1713
1714        ENDDO
1715      ENDDO
1716
1717      CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      )
1718      CALL wrk_dealloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 )
1719
1720   END SUBROUTINE s_sf12
1721
1722   SUBROUTINE s_tanh()
1723
1724      !!----------------------------------------------------------------------
1725      !!                  ***  ROUTINE s_tanh***
1726      !!                     
1727      !! ** Purpose :   stretch the s-coordinate system
1728      !!
1729      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
1730      !!                mixed S/sigma/Z coordinate
1731      !!
1732      !!
1733      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
1734      !!----------------------------------------------------------------------
1735
1736      INTEGER  ::   ji, jj, jk           ! dummy loop argument
1737      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
1738
1739      DO jk = 1, jpk
1740        gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
1741        gsigt(jk) = -fssig( REAL(jk,wp)        )
1742      END DO
1743      IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk)
1744      !
1745      ! Coefficients for vertical scale factors at w-, t- levels
1746!!gm bug :  define it from analytical function, not like juste bellow....
1747!!gm        or betteroffer the 2 possibilities....
1748      DO jk = 1, jpkm1
1749         esigt(jk  ) = gsigw(jk+1) - gsigw(jk)
1750         esigw(jk+1) = gsigt(jk+1) - gsigt(jk)
1751      END DO
1752      esigw( 1 ) = 2._wp * ( gsigt(1  ) - gsigw(1  ) ) 
1753      esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) )
1754      !
1755      ! Coefficients for vertical depth as the sum of e3w scale factors
1756      gsi3w(1) = 0.5_wp * esigw(1)
1757      DO jk = 2, jpk
1758         gsi3w(jk) = gsi3w(jk-1) + esigw(jk)
1759      END DO
1760!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
1761      DO jk = 1, jpk
1762         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1763         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1764         gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigt(jk) + hift(:,:)*zcoeft )
1765         gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigw(jk) + hift(:,:)*zcoefw )
1766         gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsi3w(jk) + hift(:,:)*zcoeft )
1767      END DO
1768!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
1769      DO jj = 1, jpj
1770         DO ji = 1, jpi
1771            DO jk = 1, jpk
1772              e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1773              e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1774              e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1775              e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
1776              !
1777              e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1778              e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1779              e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1780            END DO
1781         END DO
1782      END DO
1783   END SUBROUTINE s_tanh
1784
1785   FUNCTION fssig( pk ) RESULT( pf )
1786      !!----------------------------------------------------------------------
1787      !!                 ***  ROUTINE fssig ***
1788      !!       
1789      !! ** Purpose :   provide the analytical function in s-coordinate
1790      !!         
1791      !! ** Method  :   the function provide the non-dimensional position of
1792      !!                T and W (i.e. between 0 and 1)
1793      !!                T-points at integer values (between 1 and jpk)
1794      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1795      !!----------------------------------------------------------------------
1796      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
1797      REAL(wp)             ::   pf   ! sigma value
1798      !!----------------------------------------------------------------------
1799      !
1800      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   &
1801         &     - TANH( rn_thetb * rn_theta                                )  )   &
1802         & * (   COSH( rn_theta                           )                      &
1803         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              &
1804         & / ( 2._wp * SINH( rn_theta ) )
1805      !
1806   END FUNCTION fssig
1807
1808
1809   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
1810      !!----------------------------------------------------------------------
1811      !!                 ***  ROUTINE fssig1 ***
1812      !!
1813      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
1814      !!
1815      !! ** Method  :   the function provides the non-dimensional position of
1816      !!                T and W (i.e. between 0 and 1)
1817      !!                T-points at integer values (between 1 and jpk)
1818      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1819      !!----------------------------------------------------------------------
1820      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
1821      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
1822      REAL(wp)             ::   pf1   ! sigma value
1823      !!----------------------------------------------------------------------
1824      !
1825      IF ( rn_theta == 0 ) then      ! uniform sigma
1826         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
1827      ELSE                        ! stretched sigma
1828         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
1829            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
1830            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
1831      ENDIF
1832      !
1833   END FUNCTION fssig1
1834
1835
1836   FUNCTION fgamma( pk1, Zbb, Zss, F ) RESULT( gam )
1837      !!----------------------------------------------------------------------
1838      !!                 ***  ROUTINE fgamma  ***
1839      !!
1840      !! ** Purpose :   provide analytical function for the s-coordinate
1841      !!
1842      !! ** Method  :   the function provides the non-dimensional position of
1843      !!                T and W (i.e. between 0 and 1)
1844      !!                T-points at integer values (between 1 and jpk)
1845      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1846      !!
1847      !!                This method allows the maintenance of fixed surface and or
1848      !!                bottom cell resolutions (cf. geopotential coordinates)
1849      !!                within an analytically derived stretched S-coordinate framework.
1850      !!
1851      !! Reference  :   Siddorn and Furner, in prep
1852      !!----------------------------------------------------------------------
1853      REAL(wp), INTENT(in   ) ::   pk1(jpk)   ! continuous "k" coordinate
1854      REAL(wp)                ::   gam(jpk)   ! stretched coordinate
1855      REAL(wp), INTENT(in   ) ::   Zbb      ! Bottom box depth
1856      REAL(wp), INTENT(in   ) ::   Zss      ! surface box depth
1857      REAL(wp), INTENT(in   ) ::   F        ! Smoothing parameter
1858      REAL(wp)                ::   a1,a2,a3 ! local variables
1859      REAL(wp)                ::   n1,n2    ! local variables
1860      REAL(wp)                ::   A,B,X    ! local variables
1861      integer                 ::   jk
1862      !!----------------------------------------------------------------------
1863      !
1864
1865      n1  =  1./(jpk-1.)
1866      n2  =  1. -  n1
1867
1868      a1 = (rn_alpha+2.0_wp)*n1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*n1**(rn_alpha+2.0_wp) 
1869      a2 = (rn_alpha+2.0_wp)*n2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*n2**(rn_alpha+2.0_wp)
1870      a3 = ( n2**3.0_wp - a2)/( n1**3.0_wp - a1)
1871     
1872      A = Zbb - a3*(Zss-a1)-a2
1873      A = A/( n2-0.5_wp*(a2+n2**2.0_wp) - a3*(n1-0.5_wp*(a1+n1**2.0_wp) ) )
1874      B = (Zss - a1 - A*( n1-0.5_wp*(a1+n1**2.0_wp ) ) ) / (n1**3.0_wp - a1)
1875      X = 1.0_wp-A/2.0_wp-B
1876 
1877      DO jk = 1, jpk
1878        gam(jk) = A*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+B*pk1(jk)**3.0_wp + X*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
1879        gam(jk) = gam(jk)*F+pk1(jk)*(1.0_wp-F)
1880      ENDDO 
1881
1882      !
1883   END FUNCTION fgamma
1884
1885   !!======================================================================
1886END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.