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

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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 2435

Last change on this file since 2435 was 2435, checked in by cetlod, 13 years ago

Improve the 1D vertical configuration in v3.3beta

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