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

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

Last change on this file since 3764 was 3764, checked in by smasson, 11 years ago

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