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 NEMO/branches/NERC/dev_release-3.4_NEMOTAM_consolidated/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: NEMO/branches/NERC/dev_release-3.4_NEMOTAM_consolidated/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 11809

Last change on this file since 11809 was 3720, checked in by cbricaud, 11 years ago

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