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 trunk/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 3294

Last change on this file since 3294 was 3294, checked in by rblod, 12 years ago

Merge of 3.4beta into the trunk

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