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

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

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

Last change on this file since 2436 was 2436, checked in by gm, 13 years ago

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