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

Last change on this file since 2392 was 2392, checked in by gm, 10 years ago

v3.3beta: Cross Land Advection (ticket #127) full rewriting + MPP bug corrections

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