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

Last change on this file since 2380 was 2380, checked in by acc, 10 years ago

nemo_v3_3beta. ORCA_R1 settings (step 2, see ticket #758). Introduces key_orca_r1 (46 level default, 75 level if key_orca_r1=75)

  • 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( n_cla == 0 ) THEN
415                  !
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                  !
436               ENDIF
437               !
438            ENDIF
439            !
440         ENDIF
441         IF( ln_zps .OR. ln_sco )   THEN              ! zps or sco : read meter bathymetry
442            CALL iom_open( 'bathy_meter.nc', inum ) 
443            CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy )
444            CALL iom_close (inum)
445!                                                            ! =====================
446            IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
447               ii0 = 142   ;   ii1 = 142                     ! Close Halmera Strait 
448               ij0 =  51   ;   ij1 =  53                     ! =====================
449               DO ji = mi0(ii0), mi1(ii1)
450                  DO jj = mj0(ij0), mj1(ij1)
451                     bathy(ji,jj) = 0.0 
452                  END DO
453               END DO
454               IF(lwp) WRITE(numout,*)
455               IF(lwp) WRITE(numout,*) '             orca_r1: Halmera strait closed at i=',ii0,' j=',ij0,'->',ij1
456            ENDIF
457            !                                                ! =====================
458            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
459               !                                             ! =====================
460              IF( n_cla == 0 ) THEN
461                 !
462                 ii0 = 140   ;   ii1 = 140                  ! Gibraltar Strait open
463                 ij0 = 102   ;   ij1 = 102                  ! (Thomson, Ocean Modelling, 1995)
464                 DO ji = mi0(ii0), mi1(ii1)
465                    DO jj = mj0(ij0), mj1(ij1)
466                       bathy(ji,jj) = 284.e0
467                    END DO
468                 END DO
469                 IF(lwp) WRITE(numout,*)
470                 IF(lwp) WRITE(numout,*) '             orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
471                 !
472                 ii0 = 160   ;   ii1 = 160                  ! Bab el mandeb Strait open
473                 ij0 = 88    ;   ij1 = 88                   ! (Thomson, Ocean Modelling, 1995)
474                 DO ji = mi0(ii0), mi1(ii1)
475                    DO jj = mj0(ij0), mj1(ij1)
476                       bathy(ji,jj) = 137.e0
477                    END DO
478                 END DO
479                 IF(lwp) WRITE(numout,*)
480                 IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
481                 !
482              ENDIF
483              !
484           ENDIF
485            !
486        ENDIF
487         !                                            ! =============== !
488      ELSE                                            !      error      !
489         !                                            ! =============== !
490         WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
491         CALL ctl_stop( '    zgr_bat : '//trim(ctmp1) )
492      ENDIF
493      !
494      !                                               ! =========================== !
495      IF( nclosea == 0 ) THEN                         !   NO closed seas or lakes   !
496         DO jl = 1, jpncs                             ! =========================== !
497            DO jj = ncsj1(jl), ncsj2(jl)
498               DO ji = ncsi1(jl), ncsi2(jl)
499                  mbathy(ji,jj) = 0                   ! suppress closed seas and lakes from bathymetry
500                  bathy (ji,jj) = 0.e0               
501               END DO
502            END DO
503         END DO
504      ENDIF
505#if defined key_orca_lev10
506      !                                               ! ================= !
507      mbathy(:,:) = 10 * mbathy(:,:)                  !   key_orca_lev10  !
508      !                                               ! ================= !
509      IF( ln_zps .OR. ln_sco )   CALL ctl_stop( ' CAUTION: 300 levels only with level bathymetry' )
510#endif
511      !                                               ! =============== !
512      IF( lzoom        )   CALL zgr_bat_zoom          !   Zoom domain   !
513      !                                               ! =============== !
514
515#if ! defined key_c1d
516      !                          ! =================== !
517      CALL zgr_bat_ctl           !   Bathymetry check  !
518      !                          ! =================== !
519#endif
520   END SUBROUTINE zgr_bat
521
522
523   SUBROUTINE zgr_bat_zoom
524      !!----------------------------------------------------------------------
525      !!                    ***  ROUTINE zgr_bat_zoom  ***
526      !!
527      !! ** Purpose : - Close zoom domain boundary if necessary
528      !!              - Suppress Med Sea from ORCA R2 and R05 arctic zoom
529      !!
530      !! ** Method  :
531      !!
532      !! ** Action  : - update mbathy: level bathymetry (in level index)
533      !!----------------------------------------------------------------------
534      INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integers
535      !!----------------------------------------------------------------------
536      !
537      IF(lwp) WRITE(numout,*)
538      IF(lwp) WRITE(numout,*) '    zgr_bat_zoom : modify the level bathymetry for zoom domain'
539      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~'
540      !
541      ! Zoom domain
542      ! ===========
543      !
544      ! Forced closed boundary if required
545      IF( lzoom_s )   mbathy(  : , mj0(jpjzoom):mj1(jpjzoom) )      = 0
546      IF( lzoom_w )   mbathy(      mi0(jpizoom):mi1(jpizoom) , :  ) = 0
547      IF( lzoom_e )   mbathy(      mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , :  ) = 0
548      IF( lzoom_n )   mbathy(  : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) )      = 0
549      !
550      ! Configuration specific domain modifications
551      ! (here, ORCA arctic configuration: suppress Med Sea)
552      IF( cp_cfg == "orca" .AND. lzoom_arct ) THEN
553         SELECT CASE ( jp_cfg )
554         !                                        ! =======================
555         CASE ( 2 )                               !  ORCA_R2 configuration
556            !                                     ! =======================
557            IF(lwp) WRITE(numout,*) '                   ORCA R2 arctic zoom: suppress the Med Sea'
558            ii0 = 141   ;   ii1 = 162      ! Sea box i,j indices
559            ij0 =  98   ;   ij1 = 110
560            !                                     ! =======================
561         CASE ( 05 )                              !  ORCA_R05 configuration
562            !                                     ! =======================
563            IF(lwp) WRITE(numout,*) '                   ORCA R05 arctic zoom: suppress the Med Sea'
564            ii0 = 563   ;   ii1 = 642      ! zero over the Med Sea boxe
565            ij0 = 314   ;   ij1 = 370 
566         END SELECT
567         !
568         mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0   ! zero over the Med Sea boxe
569         !
570      ENDIF
571      !
572   END SUBROUTINE zgr_bat_zoom
573
574
575   SUBROUTINE zgr_bat_ctl
576      !!----------------------------------------------------------------------
577      !!                    ***  ROUTINE zgr_bat_ctl  ***
578      !!
579      !! ** Purpose :   check the bathymetry in levels
580      !!
581      !! ** Method  :   The array mbathy is checked to verified its consistency
582      !!      with the model options. in particular:
583      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
584      !!                  along closed boundary.
585      !!            mbathy must be cyclic IF jperio=1.
586      !!            mbathy must be lower or equal to jpk-1.
587      !!            isolated ocean grid points are suppressed from mbathy
588      !!                  since they are only connected to remaining
589      !!                  ocean through vertical diffusion.
590      !!      C A U T I O N : mbathy will be modified during the initializa-
591      !!      tion phase to become the number of non-zero w-levels of a water
592      !!      column, with a minimum value of 1.
593      !!
594      !! ** Action  : - update mbathy: level bathymetry (in level index)
595      !!              - update bathy : meter bathymetry (in meters)
596      !!----------------------------------------------------------------------
597      INTEGER ::   ji, jj, jl                    ! dummy loop indices
598      INTEGER ::   icompt, ibtest, ikmax         ! temporary integers
599      REAL(wp), DIMENSION(jpi,jpj) ::   zbathy   ! temporary workspace
600      !!----------------------------------------------------------------------
601
602      IF(lwp) WRITE(numout,*)
603      IF(lwp) WRITE(numout,*) '    zgr_bat_ctl : check the bathymetry'
604      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
605
606      !                                          ! Suppress isolated ocean grid points
607      IF(lwp) WRITE(numout,*)
608      IF(lwp) WRITE(numout,*)'                   suppress isolated ocean grid points'
609      IF(lwp) WRITE(numout,*)'                   -----------------------------------'
610      icompt = 0
611      DO jl = 1, 2
612         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
613            mbathy( 1 ,:) = mbathy(jpim1,:)           ! local domain is cyclic east-west
614            mbathy(jpi,:) = mbathy(  2  ,:)
615         ENDIF
616         DO jj = 2, jpjm1
617            DO ji = 2, jpim1
618               ibtest = MAX(  mbathy(ji-1,jj), mbathy(ji+1,jj),   &
619                  &           mbathy(ji,jj-1), mbathy(ji,jj+1)  )
620               IF( ibtest < mbathy(ji,jj) ) THEN
621                  IF(lwp) WRITE(numout,*) ' the number of ocean level at ',   &
622                     &   'grid-point (i,j) =  ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
623                  mbathy(ji,jj) = ibtest
624                  icompt = icompt + 1
625               ENDIF
626            END DO
627         END DO
628      END DO
629      IF( icompt == 0 ) THEN
630         IF(lwp) WRITE(numout,*)'     no isolated ocean grid points'
631      ELSE
632         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed'
633      ENDIF
634      IF( lk_mpp ) THEN
635         zbathy(:,:) = FLOAT( mbathy(:,:) )
636         CALL lbc_lnk( zbathy, 'T', 1. )
637         mbathy(:,:) = INT( zbathy(:,:) )
638      ENDIF
639
640      !                                          ! East-west cyclic boundary conditions
641      IF( nperio == 0 ) THEN
642         IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio
643         IF( lk_mpp ) THEN
644            IF( nbondi == -1 .OR. nbondi == 2 ) THEN
645               IF( jperio /= 1 )   mbathy(1,:) = 0
646            ENDIF
647            IF( nbondi == 1 .OR. nbondi == 2 ) THEN
648               IF( jperio /= 1 )   mbathy(nlci,:) = 0
649            ENDIF
650         ELSE
651            IF( ln_zco .OR. ln_zps ) THEN
652               mbathy( 1 ,:) = 0
653               mbathy(jpi,:) = 0
654            ELSE
655               mbathy( 1 ,:) = jpkm1
656               mbathy(jpi,:) = jpkm1
657            ENDIF
658         ENDIF
659      ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio ==  6 ) THEN
660         IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio
661         mbathy( 1 ,:) = mbathy(jpim1,:)
662         mbathy(jpi,:) = mbathy(  2  ,:)
663      ELSEIF( nperio == 2 ) THEN
664         IF(lwp) WRITE(numout,*) '   equatorial boundary conditions on mbathy: nperio = ', nperio
665      ELSE
666         IF(lwp) WRITE(numout,*) '    e r r o r'
667         IF(lwp) WRITE(numout,*) '    parameter , nperio = ', nperio
668         !         STOP 'dom_mba'
669      ENDIF
670
671      !  Boundary condition on mbathy
672      IF( .NOT.lk_mpp ) THEN 
673!!gm     !!bug ???  think about it !
674         !   ... mono- or macro-tasking: T-point, >0, 2D array, no slab
675         zbathy(:,:) = FLOAT( mbathy(:,:) )
676         CALL lbc_lnk( zbathy, 'T', 1. )
677         mbathy(:,:) = INT( zbathy(:,:) )
678      ENDIF
679
680      ! Number of ocean level inferior or equal to jpkm1
681      ikmax = 0
682      DO jj = 1, jpj
683         DO ji = 1, jpi
684            ikmax = MAX( ikmax, mbathy(ji,jj) )
685         END DO
686      END DO
687!!gm  !!! test to do:   ikmax = MAX( mbathy(:,:) )   ???
688      IF( ikmax > jpkm1 ) THEN
689         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' >  jpk-1'
690         IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
691      ELSE IF( ikmax < jpkm1 ) THEN
692         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1' 
693         IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
694      ENDIF
695
696      IF( lwp .AND. nprint == 1 ) THEN      ! control print
697         WRITE(numout,*)
698         WRITE(numout,*) ' bathymetric field :   number of non-zero T-levels '
699         WRITE(numout,*) ' ------------------'
700         CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )
701         WRITE(numout,*)
702      ENDIF
703      !
704   END SUBROUTINE zgr_bat_ctl
705
706
707   SUBROUTINE zgr_zco
708      !!----------------------------------------------------------------------
709      !!                  ***  ROUTINE zgr_zco  ***
710      !!
711      !! ** Purpose :   define the z-coordinate system
712      !!
713      !! ** Method  :   set 3D coord. arrays to reference 1D array
714      !!----------------------------------------------------------------------
715      INTEGER  ::   jk
716      !!----------------------------------------------------------------------
717      !
718      DO jk = 1, jpk
719         fsdept(:,:,jk) = gdept_0(jk)
720         fsdepw(:,:,jk) = gdepw_0(jk)
721         fsde3w(:,:,jk) = gdepw_0(jk)
722         fse3t (:,:,jk) = e3t_0(jk)
723         fse3u (:,:,jk) = e3t_0(jk)
724         fse3v (:,:,jk) = e3t_0(jk)
725         fse3f (:,:,jk) = e3t_0(jk)
726         fse3w (:,:,jk) = e3w_0(jk)
727         fse3uw(:,:,jk) = e3w_0(jk)
728         fse3vw(:,:,jk) = e3w_0(jk)
729      END DO
730      !
731   END SUBROUTINE zgr_zco
732
733   !!----------------------------------------------------------------------
734   !!   Default option :                      zco, zps and/or sco available (gedp & e3 are 3D arrays)
735   !!----------------------------------------------------------------------
736
737   SUBROUTINE zgr_zps
738      !!----------------------------------------------------------------------
739      !!                  ***  ROUTINE zgr_zps  ***
740      !!                     
741      !! ** Purpose :   the depth and vertical scale factor in partial step
742      !!      z-coordinate case
743      !!
744      !! ** Method  :   Partial steps : computes the 3D vertical scale factors
745      !!      of T-, U-, V-, W-, UW-, VW and F-points that are associated with
746      !!      a partial step representation of bottom topography.
747      !!
748      !!        The reference depth of model levels is defined from an analytical
749      !!      function the derivative of which gives the reference vertical
750      !!      scale factors.
751      !!        From  depth and scale factors reference, we compute there new value
752      !!      with partial steps  on 3d arrays ( i, j, k ).
753      !!
754      !!              w-level: gdepw(i,j,k)  = fsdep(k)
755      !!                       e3w(i,j,k) = dk(fsdep)(k)     = fse3(i,j,k)
756      !!              t-level: gdept(i,j,k)  = fsdep(k+0.5)
757      !!                       e3t(i,j,k) = dk(fsdep)(k+0.5) = fse3(i,j,k+0.5)
758      !!
759      !!        With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
760      !!      we find the mbathy index of the depth at each grid point.
761      !!      This leads us to three cases:
762      !!
763      !!              - bathy = 0 => mbathy = 0
764      !!              - 1 < mbathy < jpkm1   
765      !!              - bathy > gdepw(jpk) => mbathy = jpkm1 
766      !!
767      !!        Then, for each case, we find the new depth at t- and w- levels
768      !!      and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
769      !!      and f-points.
770      !!
771      !!        This routine is given as an example, it must be modified
772      !!      following the user s desiderata. nevertheless, the output as
773      !!      well as the way to compute the model levels and scale factors
774      !!      must be respected in order to insure second order accuracy
775      !!      schemes.
776      !!
777      !!         c a u t i o n : gdept_0, gdepw_0 and e3._0 are positives
778      !!         - - - - - - -   gdept, gdepw and e3. are positives
779      !!     
780      !!  Reference :   Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
781      !!----------------------------------------------------------------------
782      INTEGER  ::   ji, jj, jk       ! dummy loop indices
783      INTEGER  ::   ik, it           ! temporary integers
784      INTEGER  ::   nlbelow          ! temporary integer
785      LOGICAL  ::   ll_print         ! Allow  control print for debugging
786      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
787      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
788      REAL(wp) ::   zmax, zmin       ! Maximum and minimum depth
789      REAL(wp) ::   zdiff            ! temporary scalar
790      REAL(wp) ::   zrefdep          ! temporary scalar
791      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zprt   ! 3D workspace
792      !!---------------------------------------------------------------------
793
794      IF(lwp) WRITE(numout,*)
795      IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps'
796      IF(lwp) WRITE(numout,*) '    ~~~~~~~ '
797      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used'
798
799      ll_print=.FALSE.                     ! Local variable for debugging
800!!    ll_print=.TRUE.
801     
802      IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth
803         WRITE(numout,*)
804         WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)'
805         CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
806      ENDIF
807
808
809      ! bathymetry in level (from bathy_meter)
810      ! ===================
811      zmax = gdepw_0(jpk) + e3t_0(jpk)     ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) )
812
813      zrefdep = 25.
814      nlbelow = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 )   ! shallowest W level Below zrefdep
815      IF(lwp) write(numout,*) 'Minimum depth level selected: ', nlbelow
816      zmin = gdepw_0(nlbelow)              ! minimum depth = nlbelow-1 levels
817
818      mbathy(:,:) = jpkm1                  ! initialize mbathy to the maximum ocean level available
819
820      !                                    ! storage of land and island's number (zera and negative values) in mbathy
821      WHERE( bathy(:,:) <= 0. )   mbathy(:,:) = INT( bathy(:,:) )
822
823      ! bounded value of bathy
824!!gm  bathy(:,:) = MIN(  zmax,  MAX( bathy(:,:), zmin )  )     !!gm : simpler   as zmin is > 0
825      DO jj = 1, jpj
826         DO ji= 1, jpi
827            IF( bathy(ji,jj) <= 0. ) THEN   ;   bathy(ji,jj) = 0.e0
828            ELSE                            ;   bathy(ji,jj) = MIN(  zmax,  MAX( bathy(ji,jj), zmin )  )
829            ENDIF
830         END DO
831      END DO
832
833      ! Compute mbathy for ocean points (i.e. the number of ocean levels)
834      ! find the number of ocean levels such that the last level thickness
835      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_0 (where
836      ! e3t_0 is the reference level thickness
837      DO jk = jpkm1, 1, -1
838         zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat )
839         DO jj = 1, jpj
840            DO ji = 1, jpi
841               IF( 0. < bathy(ji,jj) .AND. bathy(ji,jj) <= zdepth )   mbathy(ji,jj) = jk-1
842            END DO
843         END DO
844      END DO
845
846      ! Scale factors and depth at T- and W-points
847      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate
848         gdept(:,:,jk) = gdept_0(jk)
849         gdepw(:,:,jk) = gdepw_0(jk)
850         e3t  (:,:,jk) = e3t_0  (jk)
851         e3w  (:,:,jk) = e3w_0  (jk)
852      END DO
853      hdept(:,:) = gdept(:,:,2 )
854      hdepw(:,:) = gdepw(:,:,3 )     
855      !
856      DO jj = 1, jpj
857         DO ji = 1, jpi
858            ik = mbathy(ji,jj)
859            IF( ik > 0 ) THEN               ! ocean point only
860               ! max ocean level case
861               IF( ik == jpkm1 ) THEN
862                  zdepwp = bathy(ji,jj)
863                  ze3tp  = bathy(ji,jj) - gdepw_0(ik)
864                  ze3wp = 0.5 * e3w_0(ik) * ( 1. + ( ze3tp/e3t_0(ik) ) )
865                  e3t(ji,jj,ik  ) = ze3tp
866                  e3t(ji,jj,ik+1) = ze3tp
867                  e3w(ji,jj,ik  ) = ze3wp
868                  e3w(ji,jj,ik+1) = ze3tp
869                  gdepw(ji,jj,ik+1) = zdepwp
870                  gdept(ji,jj,ik  ) = gdept_0(ik-1) + ze3wp
871                  gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp
872                  !
873               ELSE                         ! standard case
874                  IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN   ;   gdepw(ji,jj,ik+1) = bathy(ji,jj)
875                  ELSE                                       ;   gdepw(ji,jj,ik+1) = gdepw_0(ik+1)
876                  ENDIF
877!gm Bug?  check the gdepw_0
878                  !       ... on ik
879                  gdept(ji,jj,ik) = gdepw_0(ik) + ( gdepw  (ji,jj,ik+1) - gdepw_0(ik) )   &
880                     &                          * ((gdept_0(      ik  ) - gdepw_0(ik) )   &
881                     &                          / ( gdepw_0(      ik+1) - gdepw_0(ik) ))
882                  e3t  (ji,jj,ik) = e3t_0  (ik) * ( gdepw  (ji,jj,ik+1) - gdepw_0(ik) )   & 
883                     &                          / ( gdepw_0(      ik+1) - gdepw_0(ik) ) 
884                  e3w  (ji,jj,ik) = 0.5 * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2.*gdepw_0(ik) )   &
885                     &                  * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) )
886                  !       ... on ik+1
887                  e3w  (ji,jj,ik+1) = e3t  (ji,jj,ik)
888                  e3t  (ji,jj,ik+1) = e3t  (ji,jj,ik)
889                  gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)
890               ENDIF
891            ENDIF
892         END DO
893      END DO
894      !
895      it = 0
896      DO jj = 1, jpj
897         DO ji = 1, jpi
898            ik = mbathy(ji,jj)
899            IF( ik > 0 ) THEN               ! ocean point only
900               hdept(ji,jj) = gdept(ji,jj,ik  )
901               hdepw(ji,jj) = gdepw(ji,jj,ik+1)
902               e3tp (ji,jj) = e3t(ji,jj,ik  )
903               e3wp (ji,jj) = e3w(ji,jj,ik  )
904               ! test
905               zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik  )
906               IF( zdiff <= 0. .AND. lwp ) THEN
907                  it = it + 1
908                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj
909                  WRITE(numout,*) ' bathy = ', bathy(ji,jj)
910                  WRITE(numout,*) ' gdept = ', gdept(ji,jj,ik), ' gdepw = ', gdepw(ji,jj,ik+1), ' zdiff = ', zdiff
911                  WRITE(numout,*) ' e3tp  = ', e3t  (ji,jj,ik), ' e3wp  = ', e3w  (ji,jj,ik  )
912               ENDIF
913            ENDIF
914         END DO
915      END DO
916
917      ! Scale factors and depth at U-, V-, UW and VW-points
918      DO jk = 1, jpk                        ! initialisation to z-scale factors
919         e3u (:,:,jk)  = e3t_0(jk)
920         e3v (:,:,jk)  = e3t_0(jk)
921         e3uw(:,:,jk)  = e3w_0(jk)
922         e3vw(:,:,jk)  = e3w_0(jk)
923      END DO
924      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
925         DO jj = 1, jpjm1
926            DO ji = 1, fs_jpim1   ! vector opt.
927               e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk))
928               e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk))
929               e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) )
930               e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) )
931            END DO
932         END DO
933      END DO
934      !                                     ! Boundary conditions
935      CALL lbc_lnk( e3u , 'U', 1. )   ;   CALL lbc_lnk( e3uw, 'U', 1. )
936      CALL lbc_lnk( e3v , 'V', 1. )   ;   CALL lbc_lnk( e3vw, 'V', 1. )
937      !
938      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
939         WHERE( e3u (:,:,jk) == 0.e0 )   e3u (:,:,jk) = e3t_0(jk)
940         WHERE( e3v (:,:,jk) == 0.e0 )   e3v (:,:,jk) = e3t_0(jk)
941         WHERE( e3uw(:,:,jk) == 0.e0 )   e3uw(:,:,jk) = e3w_0(jk)
942         WHERE( e3vw(:,:,jk) == 0.e0 )   e3vw(:,:,jk) = e3w_0(jk)
943      END DO
944     
945      ! Scale factor at F-point
946      DO jk = 1, jpk                        ! initialisation to z-scale factors
947         e3f (:,:,jk) = e3t_0(jk)
948      END DO
949      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors
950         DO jj = 1, jpjm1
951            DO ji = 1, fs_jpim1   ! vector opt.
952               e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) )
953            END DO
954         END DO
955      END DO
956      CALL lbc_lnk( e3f, 'F', 1. )          ! Boundary conditions
957      !
958      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
959         WHERE( e3f(:,:,jk) == 0.e0 )   e3f(:,:,jk) = e3t_0(jk)
960      END DO
961!!gm  bug ? :  must be a do loop with mj0,mj1
962      !
963      e3t(:,mj0(1),:) = e3t(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2
964      e3w(:,mj0(1),:) = e3w(:,mj0(2),:) 
965      e3u(:,mj0(1),:) = e3u(:,mj0(2),:) 
966      e3v(:,mj0(1),:) = e3v(:,mj0(2),:) 
967      e3f(:,mj0(1),:) = e3f(:,mj0(2),:) 
968
969      ! Control of the sign
970      IF( MINVAL( e3t  (:,:,:) ) <= 0.e0 )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t   <= 0' )
971      IF( MINVAL( e3w  (:,:,:) ) <= 0.e0 )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w   <= 0' )
972      IF( MINVAL( gdept(:,:,:) ) <  0.e0 )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw <  0' )
973      IF( MINVAL( gdepw(:,:,:) ) <  0.e0 )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw <  0' )
974     
975      ! Compute gdep3w (vertical sum of e3w)
976      gdep3w(:,:,1) = 0.5 * e3w(:,:,1)
977      DO jk = 2, jpk
978         gdep3w(:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk) 
979      END DO
980         
981      !                                               ! ================= !
982      IF(lwp .AND. ll_print) THEN                     !   Control print   !
983         !                                            ! ================= !
984         DO jj = 1,jpj
985            DO ji = 1, jpi
986               ik = MAX( mbathy(ji,jj), 1 )
987               zprt(ji,jj,1) = e3t   (ji,jj,ik)
988               zprt(ji,jj,2) = e3w   (ji,jj,ik)
989               zprt(ji,jj,3) = e3u   (ji,jj,ik)
990               zprt(ji,jj,4) = e3v   (ji,jj,ik)
991               zprt(ji,jj,5) = e3f   (ji,jj,ik)
992               zprt(ji,jj,6) = gdep3w(ji,jj,ik)
993            END DO
994         END DO
995         WRITE(numout,*)
996         WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
997         WRITE(numout,*)
998         WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
999         WRITE(numout,*)
1000         WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1001         WRITE(numout,*)
1002         WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1003         WRITE(numout,*)
1004         WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1005         WRITE(numout,*)
1006         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1007      ENDIF 
1008       
1009      !                                               ! =============== !
1010      IF( lzoom )   CALL zgr_bat_zoom                 !   Zoom domain   !
1011      !                                               ! =============== !
1012
1013#if ! defined key_c1d
1014      !                          ! =================== !
1015      CALL zgr_bat_ctl           !   Bathymetry check  !
1016      !                          ! =================== !
1017#endif
1018   END SUBROUTINE zgr_zps
1019
1020
1021   FUNCTION fssig( pk ) RESULT( pf )
1022      !!----------------------------------------------------------------------
1023      !!                 ***  ROUTINE eos_init  ***
1024      !!       
1025      !! ** Purpose :   provide the analytical function in s-coordinate
1026      !!         
1027      !! ** Method  :   the function provide the non-dimensional position of
1028      !!                T and W (i.e. between 0 and 1)
1029      !!                T-points at integer values (between 1 and jpk)
1030      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1031      !!
1032      !! Reference  :   ???
1033      !!----------------------------------------------------------------------
1034      REAL(wp), INTENT(in   ) ::   pk   ! continuous "k" coordinate
1035      REAL(wp)                ::   pf   ! sigma value
1036      !!----------------------------------------------------------------------
1037      !
1038      pf =   (   TANH( rn_theta * ( -(pk-0.5) / REAL(jpkm1) + rn_thetb )  )      &
1039         &     - TANH( rn_thetb * rn_theta                                )  )   &
1040         & * (   COSH( rn_theta                           )                   &
1041         &     + COSH( rn_theta * ( 2.e0 * rn_thetb - 1.e0 ) )  )                &
1042         & / ( 2.e0 * SINH( rn_theta ) )
1043      !
1044   END FUNCTION fssig
1045
1046
1047   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
1048      !!----------------------------------------------------------------------
1049      !!                 ***  ROUTINE eos_init  ***
1050      !!
1051      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
1052      !!
1053      !! ** Method  :   the function provides the non-dimensional position of
1054      !!                T and W (i.e. between 0 and 1)
1055      !!                T-points at integer values (between 1 and jpk)
1056      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1057      !!
1058      !! Reference  :   ???
1059      !!----------------------------------------------------------------------
1060      REAL(wp), INTENT(in   ) ::   pk1   ! continuous "k" coordinate
1061      REAL(wp), INTENT(in   ) ::   pbb   ! Stretching coefficient
1062      REAL(wp)                ::   pf1   ! sigma value
1063      !!----------------------------------------------------------------------
1064      !
1065      IF ( rn_theta == 0 ) then      ! uniform sigma
1066         pf1 = -(pk1-0.5) / REAL( jpkm1 )
1067      ELSE                        ! stretched sigma
1068         pf1 =   (1.0-pbb) * (sinh( rn_theta*(-(pk1-0.5)/REAL(jpkm1)) ) ) / sinh(rn_theta) + &
1069            &    pbb * ( (tanh( rn_theta*( (-(pk1-0.5)/REAL(jpkm1)) + 0.5) ) - tanh(0.5*rn_theta) ) / &
1070            &    (2*tanh(0.5*rn_theta) ) )
1071      ENDIF
1072      !
1073   END FUNCTION fssig1
1074
1075
1076   SUBROUTINE zgr_sco
1077      !!----------------------------------------------------------------------
1078      !!                  ***  ROUTINE zgr_sco  ***
1079      !!                     
1080      !! ** Purpose :   define the s-coordinate system
1081      !!
1082      !! ** Method  :   s-coordinate
1083      !!         The depth of model levels is defined as the product of an
1084      !!      analytical function by the local bathymetry, while the vertical
1085      !!      scale factors are defined as the product of the first derivative
1086      !!      of the analytical function by the bathymetry.
1087      !!      (this solution save memory as depth and scale factors are not
1088      !!      3d fields)
1089      !!          - Read bathymetry (in meters) at t-point and compute the
1090      !!         bathymetry at u-, v-, and f-points.
1091      !!            hbatu = mi( hbatt )
1092      !!            hbatv = mj( hbatt )
1093      !!            hbatf = mi( mj( hbatt ) )
1094      !!          - Compute gsigt, gsigw, esigt, esigw from an analytical
1095      !!         function and its derivative given as function.
1096      !!            gsigt(k) = fssig (k    )
1097      !!            gsigw(k) = fssig (k-0.5)
1098      !!            esigt(k) = fsdsig(k    )
1099      !!            esigw(k) = fsdsig(k-0.5)
1100      !!      This routine is given as an example, it must be modified
1101      !!      following the user s desiderata. nevertheless, the output as
1102      !!      well as the way to compute the model levels and scale factors
1103      !!      must be respected in order to insure second order a!!uracy
1104      !!      schemes.
1105      !!
1106      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
1107      !!----------------------------------------------------------------------
1108      INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
1109      INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
1110      REAL(wp) ::   zcoeft, zcoefw, zrmax, ztaper   ! temporary scalars
1111      REAL(wp), DIMENSION(jpi,jpj) ::   zenv, ztmp, zmsk    ! 2D workspace
1112      REAL(wp), DIMENSION(jpi,jpj) ::   zri , zrj , zhbat   !  -     -
1113      !!
1114      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   gsigw3
1115      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   gsigt3
1116      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   gsi3w3
1117      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   esigt3
1118      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   esigw3
1119      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   esigtu3
1120      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   esigtv3
1121      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   esigtf3
1122      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   esigwu3
1123      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   esigwv3
1124      !!
1125      NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc
1126      !!----------------------------------------------------------------------
1127
1128      REWIND( numnam )                        ! Read Namelist namzgr_sco : sigma-stretching parameters
1129      READ  ( numnam, namzgr_sco )
1130
1131      IF(lwp) THEN                            ! control print
1132         WRITE(numout,*)
1133         WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate'
1134         WRITE(numout,*) '~~~~~~~~~~~'
1135         WRITE(numout,*) '   Namelist namzgr_sco'
1136         WRITE(numout,*) '      sigma-stretching coeffs '
1137         WRITE(numout,*) '      maximum depth of s-bottom surface (>0)       rn_sbot_max   = ' ,rn_sbot_max
1138         WRITE(numout,*) '      minimum depth of s-bottom surface (>0)       rn_sbot_min   = ' ,rn_sbot_min
1139         WRITE(numout,*) '      surface control parameter (0<=rn_theta<=20)  rn_theta      = ', rn_theta
1140         WRITE(numout,*) '      bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ', rn_thetb
1141         WRITE(numout,*) '      maximum cut-off r-value allowed              rn_rmax       = ', rn_rmax
1142         WRITE(numout,*) '      Hybrid s-sigma-coordinate                    ln_s_sigma    = ', ln_s_sigma
1143         WRITE(numout,*) '      stretching parameter (song and haidvogel)    rn_bb         = ', rn_bb
1144         WRITE(numout,*) '      Critical depth                               rn_hc         = ', rn_hc
1145      ENDIF
1146
1147      gsigw3  = 0.0d0   ;   gsigt3  = 0.0d0   ; gsi3w3 = 0.0d0
1148      esigt3  = 0.0d0   ;   esigw3  = 0.0d0 
1149      esigtu3 = 0.0d0   ;   esigtv3 = 0.0d0   ; esigtf3 = 0.0d0
1150      esigwu3 = 0.0d0   ;   esigwv3 = 0.0d0
1151
1152      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate
1153      hifu(:,:) = rn_sbot_min
1154      hifv(:,:) = rn_sbot_min
1155      hiff(:,:) = rn_sbot_min
1156
1157      !                                        ! set maximum ocean depth
1158      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1159
1160      DO jj = 1, jpj
1161         DO ji = 1, jpi
1162           IF (bathy(ji,jj) .gt. 0.0) THEN
1163              bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1164           ENDIF
1165         END DO
1166      END DO
1167      !                                        ! =============================
1168      !                                        ! Define the envelop bathymetry   (hbatt)
1169      !                                        ! =============================
1170      ! use r-value to create hybrid coordinates
1171      DO jj = 1, jpj
1172         DO ji = 1, jpi
1173            zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min )
1174         END DO
1175      END DO
1176      !
1177      ! Smooth the bathymetry (if required)
1178      scosrf(:,:) = 0.e0              ! ocean surface depth (here zero: no under ice-shelf sea)
1179      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
1180      !
1181      jl = 0
1182      zrmax = 1.e0
1183      !                                                     ! ================ !
1184      DO WHILE ( jl <= 10000 .AND. zrmax > rn_rmax )          !  Iterative loop  !
1185         !                                                  ! ================ !
1186         jl = jl + 1
1187         zrmax = 0.e0
1188         zmsk(:,:) = 0.e0
1189         DO jj = 1, nlcj
1190            DO ji = 1, nlci
1191               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
1192               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
1193               zri(ji,jj) = ABS( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
1194               zrj(ji,jj) = ABS( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
1195               zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) )
1196               IF( zri(ji,jj) > rn_rmax )   zmsk(ji  ,jj  ) = 1.0
1197               IF( zri(ji,jj) > rn_rmax )   zmsk(iip1,jj  ) = 1.0
1198               IF( zrj(ji,jj) > rn_rmax )   zmsk(ji  ,jj  ) = 1.0
1199               IF( zrj(ji,jj) > rn_rmax )   zmsk(ji  ,ijp1) = 1.0
1200            END DO
1201         END DO
1202         IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
1203         ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)
1204         ztmp(:,:) = zmsk(:,:)   ;   CALL lbc_lnk( zmsk, 'T', 1. )
1205         DO jj = 1, nlcj
1206            DO ji = 1, nlci
1207                zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )
1208            END DO
1209         END DO
1210         !
1211         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) )
1212         !
1213         DO jj = 1, nlcj
1214            DO ji = 1, nlci
1215               iip1 = MIN( ji+1, nlci )     ! last  line (ji=nlci)
1216               ijp1 = MIN( jj+1, nlcj )     ! last  raw  (jj=nlcj)
1217               iim1 = MAX( ji-1,  1  )      ! first line (ji=nlci)
1218               ijm1 = MAX( jj-1,  1  )      ! first raw  (jj=nlcj)
1219               IF( zmsk(ji,jj) == 1.0 ) THEN
1220                  ztmp(ji,jj) =   (                                                                                   &
1221             &      zenv(iim1,ijp1)*zmsk(iim1,ijp1) + zenv(ji,ijp1)*zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1)   &
1222             &    + zenv(iim1,jj  )*zmsk(iim1,jj  ) + zenv(ji,jj  )*    2.e0      + zenv(iip1,jj  )*zmsk(iip1,jj  )   &
1223             &    + zenv(iim1,ijm1)*zmsk(iim1,ijm1) + zenv(ji,ijm1)*zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1)   &
1224             &                    ) / (                                                                               &
1225             &                      zmsk(iim1,ijp1) +               zmsk(ji,ijp1) +                 zmsk(iip1,ijp1)   &
1226             &    +                 zmsk(iim1,jj  ) +                   2.e0      +                 zmsk(iip1,jj  )   &
1227             &    +                 zmsk(iim1,ijm1) +               zmsk(ji,ijm1) +                 zmsk(iip1,ijm1)   &
1228             &                        )
1229               ENDIF
1230            END DO
1231         END DO
1232         !
1233         DO jj = 1, nlcj
1234            DO ji = 1, nlci
1235               IF( zmsk(ji,jj) == 1.0 )   zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )
1236            END DO
1237         END DO
1238         !
1239         ! Apply lateral boundary condition   CAUTION: kept the value when the lbc field is zero
1240         ztmp(:,:) = zenv(:,:)   ;   CALL lbc_lnk( zenv, 'T', 1. )
1241         DO jj = 1, nlcj
1242            DO ji = 1, nlci
1243               IF( zenv(ji,jj) == 0.e0 )   zenv(ji,jj) = ztmp(ji,jj)
1244            END DO
1245         END DO
1246         !                                                  ! ================ !
1247      END DO                                                !     End loop     !
1248      !                                                     ! ================ !
1249      !
1250      !                                        ! envelop bathymetry saved in hbatt
1251      hbatt(:,:) = zenv(:,:) 
1252      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0.e0 ) THEN
1253         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1254         DO jj = 1, jpj
1255            DO ji = 1, jpi
1256               ztaper = EXP( -(gphit(ji,jj)/8)**2 )
1257               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * (1.-ztaper)
1258            END DO
1259         END DO
1260      ENDIF
1261      !
1262      IF(lwp) THEN                             ! Control print
1263         WRITE(numout,*)
1264         WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
1265         WRITE(numout,*)
1266         CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0., numout )
1267         IF( nprint == 1 )   THEN       
1268            WRITE(numout,*) ' bathy  MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
1269            WRITE(numout,*) ' hbatt  MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
1270         ENDIF
1271      ENDIF
1272
1273      !                                        ! ==============================
1274      !                                        !   hbatu, hbatv, hbatf fields
1275      !                                        ! ==============================
1276      IF(lwp) THEN
1277         WRITE(numout,*)
1278         WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
1279      ENDIF
1280      hbatu(:,:) = rn_sbot_min
1281      hbatv(:,:) = rn_sbot_min
1282      hbatf(:,:) = rn_sbot_min
1283      DO jj = 1, jpjm1
1284        DO ji = 1, jpim1   ! NO vector opt.
1285           hbatu(ji,jj) = 0.5 * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
1286           hbatv(ji,jj) = 0.5 * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
1287           hbatf(ji,jj) = 0.25* ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
1288              &                 + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
1289        END DO
1290      END DO
1291      !
1292      ! Apply lateral boundary condition
1293!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
1294      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk( hbatu, 'U', 1. )
1295      DO jj = 1, jpj
1296         DO ji = 1, jpi
1297            IF( hbatu(ji,jj) == 0.e0 ) THEN
1298               IF( zhbat(ji,jj) == 0.e0 )   hbatu(ji,jj) = rn_sbot_min
1299               IF( zhbat(ji,jj) /= 0.e0 )   hbatu(ji,jj) = zhbat(ji,jj)
1300            ENDIF
1301         END DO
1302      END DO
1303      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk( hbatv, 'V', 1. )
1304      DO jj = 1, jpj
1305         DO ji = 1, jpi
1306            IF( hbatv(ji,jj) == 0.e0 ) THEN
1307               IF( zhbat(ji,jj) == 0.e0 )   hbatv(ji,jj) = rn_sbot_min
1308               IF( zhbat(ji,jj) /= 0.e0 )   hbatv(ji,jj) = zhbat(ji,jj)
1309            ENDIF
1310         END DO
1311      END DO
1312      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk( hbatf, 'F', 1. )
1313      DO jj = 1, jpj
1314         DO ji = 1, jpi
1315            IF( hbatf(ji,jj) == 0.e0 ) THEN
1316               IF( zhbat(ji,jj) == 0.e0 )   hbatf(ji,jj) = rn_sbot_min
1317               IF( zhbat(ji,jj) /= 0.e0 )   hbatf(ji,jj) = zhbat(ji,jj)
1318            ENDIF
1319         END DO
1320      END DO
1321
1322!!bug:  key_helsinki a verifer
1323      hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1324      hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1325      hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1326      hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1327
1328      IF( nprint == 1 .AND. lwp )   THEN
1329         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
1330            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
1331         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
1332            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
1333         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
1334            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1335         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
1336            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1337      ENDIF
1338!! helsinki
1339
1340      !                                            ! =======================
1341      !                                            !   s-ccordinate fields     (gdep., e3.)
1342      !                                            ! =======================
1343      !
1344      ! non-dimensional "sigma" for model level depth at w- and t-levels
1345
1346      IF( ln_s_sigma ) THEN        ! Song and Haidvogel style stretched sigma for depths
1347         !                         ! below rn_hc, with uniform sigma in shallower waters
1348         DO ji = 1, jpi
1349            DO jj = 1, jpj
1350
1351             IF (hbatt(ji,jj).GT.rn_hc) THEN !deep water, stretched sigma
1352               DO jk = 1, jpk
1353                  gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
1354                  gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
1355               END DO
1356             ELSE ! shallow water, uniform sigma
1357               DO jk = 1, jpk
1358                 gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
1359                 gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
1360               END DO
1361             ENDIF
1362             IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw3 1 jpk    ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk)
1363
1364
1365             DO jk = 1, jpkm1
1366                esigt3(ji,jj,jk) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk)
1367                esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk)
1368             END DO
1369             esigw3(ji,jj,1  ) = 2.0_wp * (gsigt3(ji,jj,1  ) - gsigw3(ji,jj,1  ))
1370             esigt3(ji,jj,jpk) = 2.0_wp * (gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk))
1371
1372             ! Coefficients for vertical depth as the sum of e3w scale factors
1373             gsi3w3(ji,jj,1) = 0.5 * esigw3(ji,jj,1)
1374             DO jk = 2, jpk
1375                gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk)
1376             END DO
1377
1378             DO jk = 1, jpk
1379                zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpkm1,wp)
1380                zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpkm1,wp)
1381                gdept (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft)
1382                gdepw (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw)
1383                gdep3w(ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft)
1384             END DO
1385
1386           ENDDO   ! for all jj's
1387         ENDDO    ! for all ji's
1388
1389
1390         DO ji=1,jpi
1391           DO jj=1,jpj
1392
1393             DO jk = 1, jpk
1394                esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) ) / &
1395                                   ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1396                esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) ) / &
1397                                   ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1398                esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) +  &
1399                                     hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) ) / &
1400                                   ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
1401                esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) ) / &
1402                                   ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1403                esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) ) / &
1404                                   ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1405
1406                e3t(ji,jj,jk)=((hbatt(ji,jj)-rn_hc)*esigt3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
1407                e3u(ji,jj,jk)=((hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
1408                e3v(ji,jj,jk)=((hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
1409                e3f(ji,jj,jk)=((hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
1410                !
1411                e3w (ji,jj,jk)=((hbatt(ji,jj)-rn_hc)*esigw3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
1412                e3uw(ji,jj,jk)=((hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
1413                e3vw(ji,jj,jk)=((hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
1414             END DO
1415
1416           ENDDO
1417         ENDDO
1418
1419      ELSE   ! not ln_s_sigma
1420
1421         DO jk = 1, jpk
1422           gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
1423           gsigt(jk) = -fssig( REAL(jk,wp)     )
1424         END DO
1425         IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk)
1426         !
1427     ! Coefficients for vertical scale factors at w-, t- levels
1428!!gm bug :  define it from analytical function, not like juste bellow....
1429!!gm        or betteroffer the 2 possibilities....
1430         DO jk = 1, jpkm1
1431            esigt(jk  ) = gsigw(jk+1) - gsigw(jk)
1432            esigw(jk+1) = gsigt(jk+1) - gsigt(jk)
1433         END DO
1434         esigw( 1 ) = 2.0_wp * (gsigt(1) - gsigw(1)) 
1435         esigt(jpk) = 2.0_wp * (gsigt(jpk) - gsigw(jpk))
1436
1437!!gm  original form
1438!!org DO jk = 1, jpk
1439!!org    esigt(jk)=fsdsig( FLOAT(jk)     )
1440!!org    esigw(jk)=fsdsig( FLOAT(jk)-0.5 )
1441!!org END DO
1442!!gm
1443         !
1444         ! Coefficients for vertical depth as the sum of e3w scale factors
1445         gsi3w(1) = 0.5 * esigw(1)
1446         DO jk = 2, jpk
1447            gsi3w(jk) = gsi3w(jk-1) + esigw(jk)
1448         END DO
1449!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
1450         DO jk = 1, jpk
1451            zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1)
1452            zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1)
1453            gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigt(jk)+hift(:,:)*zcoeft)
1454            gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigw(jk)+hift(:,:)*zcoefw)
1455            gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsi3w(jk)+hift(:,:)*zcoeft)
1456         END DO
1457!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
1458         DO jj = 1, jpj
1459            DO ji = 1, jpi
1460               DO jk = 1, jpk
1461                 e3t(ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/FLOAT(jpkm1))
1462                 e3u(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/FLOAT(jpkm1))
1463                 e3v(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/FLOAT(jpkm1))
1464                 e3f(ji,jj,jk)=((hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/FLOAT(jpkm1))
1465                 !
1466                 e3w (ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/FLOAT(jpkm1))
1467                 e3uw(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/FLOAT(jpkm1))
1468                 e3vw(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/FLOAT(jpkm1))
1469               END DO
1470            END DO
1471         END DO
1472
1473      ENDIF ! ln_s_sigma
1474
1475
1476      !
1477!!    H. Liu, POL. April 2009. Added for passing the scale check for the new released vvl code.
1478
1479      fsdept(:,:,:) = gdept (:,:,:)
1480      fsdepw(:,:,:) = gdepw (:,:,:)
1481      fsde3w(:,:,:) = gdep3w(:,:,:)
1482      fse3t (:,:,:) = e3t   (:,:,:)
1483      fse3u (:,:,:) = e3u   (:,:,:)
1484      fse3v (:,:,:) = e3v   (:,:,:)
1485      fse3f (:,:,:) = e3f   (:,:,:)
1486      fse3w (:,:,:) = e3w   (:,:,:)
1487      fse3uw(:,:,:) = e3uw  (:,:,:)
1488      fse3vw(:,:,:) = e3vw  (:,:,:)
1489!!
1490      ! HYBRID :
1491      DO jj = 1, jpj
1492         DO ji = 1, jpi
1493            DO jk = 1, jpkm1
1494               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
1495               IF( scobot(ji,jj) == 0.e0             )   mbathy(ji,jj) = 0
1496            END DO
1497         END DO
1498      END DO
1499      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
1500         &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
1501
1502
1503      !                                               ! ===========
1504      IF( lzoom )   CALL zgr_bat_zoom                 ! Zoom domain
1505      !                                               ! ===========
1506
1507#if ! defined key_c1d
1508      !                          ! =================== !
1509      CALL zgr_bat_ctl           !   Bathymetry check  !
1510      !                          ! =================== !
1511#endif
1512
1513      !                                               ! =============
1514      IF(lwp) THEN                                    ! Control print
1515         !                                            ! =============
1516         WRITE(numout,*) 
1517         WRITE(numout,*) ' domzgr: vertical coefficients for model level'
1518         WRITE(numout, "(9x,'  level    gsigt      gsigw      esigt      esigw      gsi3w')" )
1519         WRITE(numout, "(10x,i4,5f11.4)" ) ( jk, gsigt(jk), gsigw(jk), esigt(jk), esigw(jk), gsi3w(jk), jk=1,jpk )
1520      ENDIF
1521      IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain
1522         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)   ), ' MAX ', MAXVAL( mbathy(:,:) )
1523         WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ),   &
1524            &                          ' w ', MINVAL( fsdepw(:,:,:) ), '3w '  , MINVAL( fsde3w(:,:,:) )
1525         WRITE(numout,*) ' MIN val e3    t ', MINVAL( fse3t (:,:,:) ), ' f '  , MINVAL( fse3f (:,:,:) ),   &
1526            &                          ' u ', MINVAL( fse3u (:,:,:) ), ' u '  , MINVAL( fse3v (:,:,:) ),   &
1527            &                          ' uw', MINVAL( fse3uw(:,:,:) ), ' vw'  , MINVAL( fse3vw(:,:,:) ),   &
1528            &                          ' w ', MINVAL( fse3w (:,:,:) )
1529
1530         WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ),   &
1531            &                          ' w ', MAXVAL( fsdepw(:,:,:) ), '3w '  , MAXVAL( fsde3w(:,:,:) )
1532         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( fse3t (:,:,:) ), ' f '  , MAXVAL( fse3f (:,:,:) ),   &
1533            &                          ' u ', MAXVAL( fse3u (:,:,:) ), ' u '  , MAXVAL( fse3v (:,:,:) ),   &
1534            &                          ' uw', MAXVAL( fse3uw(:,:,:) ), ' vw'  , MAXVAL( fse3vw(:,:,:) ),   &
1535            &                          ' w ', MAXVAL( fse3w (:,:,:) )
1536      ENDIF
1537      !
1538      IF(lwp) THEN                                  ! selected vertical profiles
1539         WRITE(numout,*)
1540         WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
1541         WRITE(numout,*) ' ~~~~~~  --------------------'
1542         WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
1543         WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk),     &
1544            &                                 fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk )
1545         DO jj = mj0(20), mj1(20)
1546            DO ji = mi0(20), mi1(20)
1547               WRITE(numout,*)
1548               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
1549               WRITE(numout,*) ' ~~~~~~  --------------------'
1550               WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
1551               WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk),     &
1552                  &                                 fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
1553            END DO
1554         END DO
1555         DO jj = mj0(74), mj1(74)
1556            DO ji = mi0(100), mi1(100)
1557               WRITE(numout,*)
1558               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
1559               WRITE(numout,*) ' ~~~~~~  --------------------'
1560               WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
1561               WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk),     &
1562                  &                                 fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
1563            END DO
1564         END DO
1565      ENDIF
1566
1567!!gm bug?  no more necessary?  if ! defined key_helsinki
1568      DO jk = 1, jpk
1569         DO jj = 1, jpj
1570            DO ji = 1, jpi
1571               IF( fse3w(ji,jj,jk) <= 0. .OR. fse3t(ji,jj,jk) <= 0. ) THEN
1572                  WRITE(ctmp1,*) 'zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
1573                  CALL ctl_stop( ctmp1 )
1574               ENDIF
1575               IF( fsdepw(ji,jj,jk) < 0. .OR. fsdept(ji,jj,jk) < 0. ) THEN
1576                  WRITE(ctmp1,*) 'zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
1577                  CALL ctl_stop( ctmp1 )
1578               ENDIF
1579            END DO
1580         END DO
1581      END DO
1582!!gm bug    #endif
1583      !
1584   END SUBROUTINE zgr_sco
1585
1586
1587   !!======================================================================
1588END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.