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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 2633

Last change on this file since 2633 was 2633, checked in by trackstand2, 12 years ago

Renamed wrk_use => wrk_in_use and wrk_release => wrk_not_released

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