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.2.F90 on Ticket #671 – Attachment – NEMO

Ticket #671: domzgr.2.F90

File domzgr.2.F90, 76.4 KB (added by gm, 13 years ago)

corrected domzgr.F90 by Gurvan

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