Ticket #671: domzgr.F90

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