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 @ 2683

Last change on this file since 2683 was 2683, checked in by gm, 12 years ago

dynamic mem: #785 ; add the allocation of s-coord arrays which were missing

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