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

source: branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 6736

Last change on this file since 6736 was 6736, checked in by jamesharle, 8 years ago

FASTNEt code modifications

  • Property svn:keywords set to Id
File size: 117.2 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   !!            3.4  ! 2012-12  (R. Bourdalle-Badie and G. Reffray)  modify C1D case 
18  !!----------------------------------------------------------------------
19
20   !!----------------------------------------------------------------------
21   !!   dom_zgr          : defined the ocean vertical coordinate system
22   !!       zgr_bat      : bathymetry fields (levels and meters)
23   !!       zgr_bat_zoom : modify the bathymetry field if zoom domain
24   !!       zgr_bat_ctl  : check the bathymetry files
25   !!       zgr_bot_level: deepest ocean level for t-, u, and v-points
26   !!       zgr_z        : reference z-coordinate
27   !!       zgr_zco      : z-coordinate
28   !!       zgr_zps      : z-coordinate with partial steps
29   !!       zgr_sco      : s-coordinate
30   !!       fssig        : sigma coordinate non-dimensional function
31   !!       dfssig       : derivative of the sigma coordinate function    !!gm  (currently missing!)
32   !!---------------------------------------------------------------------
33   USE oce               ! ocean variables
34   USE dom_oce           ! ocean domain
35   USE closea            ! closed seas
36   USE c1d               ! 1D vertical configuration
37   USE in_out_manager    ! I/O manager
38   USE iom               ! I/O library
39   USE lbclnk            ! ocean lateral boundary conditions (or mpp link)
40   USE lib_mpp           ! distributed memory computing library
41   USE wrk_nemo          ! Memory allocation
42   USE timing            ! Timing
43
44   IMPLICIT NONE
45   PRIVATE
46
47   PUBLIC   dom_zgr        ! called by dom_init.F90
48
49   !                                       !!* Namelist namzgr_sco *
50   REAL(wp) ::   rn_sbot_min =  300._wp     ! minimum depth of s-bottom surface (>0) (m)
51   REAL(wp) ::   rn_sbot_max = 5250._wp     ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
52   REAL(wp) ::   rn_theta    =    6.00_wp   ! surface control parameter (0<=rn_theta<=20)
53   REAL(wp) ::   rn_thetb    =    0.75_wp   ! bottom control parameter  (0<=rn_thetb<= 1)
54   REAL(wp) ::   rn_rmax     =    0.15_wp   ! maximum cut-off r-value allowed (0<rn_rmax<1)
55!  LOGICAL  ::   ln_s_sigma  = .false.      ! use hybrid s-sigma -coordinate & stretching function fssig1 (ln_sco=T)
56   REAL(wp) ::   rn_bb       =    0.80_wp   ! stretching parameter for song and haidvogel stretching
57   !                                        ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
58   REAL(wp) ::   rn_hc       =  150._wp     ! Critical depth for s-sigma coordinates
59   REAL(wp) ::   rn_zsigma   =  300._wp     ! Maximum  depth for s-sigma layer
60   INTEGER  ::   nn_sig_lev  =  10          ! Maximum number of levels of s-sigma layer
61   REAL(wp) ::   rn_kth      = 15._wp    ! Approximate layer number, beyond which streching will be maximum
62   REAL(wp) ::   rn_acr      = 9.00_wp   !
63
64  !! * Substitutions
65#  include "domzgr_substitute.h90"
66#  include "vectopt_loop_substitute.h90"
67   !!----------------------------------------------------------------------
68   !! NEMO/OPA 3.3.1 , NEMO Consortium (2011)
69   !! $Id$
70   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
71   !!----------------------------------------------------------------------
72CONTAINS       
73
74   SUBROUTINE dom_zgr
75      !!----------------------------------------------------------------------
76      !!                ***  ROUTINE dom_zgr  ***
77      !!                   
78      !! ** Purpose :   set the depth of model levels and the resulting
79      !!              vertical scale factors.
80      !!
81      !! ** Method  : - reference 1D vertical coordinate (gdep._0, e3._0)
82      !!              - read/set ocean depth and ocean levels (bathy, mbathy)
83      !!              - vertical coordinate (gdep., e3.) depending on the
84      !!                coordinate chosen :
85      !!                   ln_zco=T   z-coordinate   
86      !!                   ln_zps=T   z-coordinate with partial steps
87      !!                   ln_zco=T   s-coordinate
88      !!
89      !! ** Action  :   define gdep., e3., mbathy and bathy
90      !!----------------------------------------------------------------------
91      INTEGER ::   ioptio, ibat   ! local integer
92      !
93      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco , ln_hyb
94      !!----------------------------------------------------------------------
95      !
96      IF( nn_timing == 1 )   CALL timing_start('dom_zgr')
97      !
98      REWIND( numnam )                 ! Read Namelist namzgr : vertical coordinate'
99      READ  ( numnam, namzgr )
100
101      IF(lwp) THEN                     ! Control print
102         WRITE(numout,*)
103         WRITE(numout,*) 'dom_zgr : vertical coordinate'
104         WRITE(numout,*) '~~~~~~~'
105         WRITE(numout,*) '          Namelist namzgr : set vertical coordinate'
106         WRITE(numout,*) '             z-coordinate - full steps      ln_zco = ', ln_zco
107         WRITE(numout,*) '             z-coordinate - partial steps   ln_zps = ', ln_zps
108         WRITE(numout,*) '             s- or hybrid z-s-coordinate    ln_sco = ', ln_sco
109         WRITE(numout,*) '             hybrid s-z-coordinates,s at shelf ln_hyb = ', ln_hyb
110
111      ENDIF
112
113      ioptio = 0                       ! Check Vertical coordinate options
114      IF( ln_zco      )   ioptio = ioptio + 1
115      IF( ln_zps      )   ioptio = ioptio + 1
116      IF( ln_sco      )   ioptio = ioptio + 1
117      IF( ioptio /= 1 )   CALL ctl_stop( ' none or several vertical coordinate options used' )
118      !
119      ! Build the vertical coordinate system
120      ! ------------------------------------
121                          CALL zgr_z            ! Reference z-coordinate system (always called)
122                          CALL zgr_bat          ! Bathymetry fields (levels and meters)
123      IF( lk_c1d      )   CALL lbc_lnk( bathy , 'T', 1._wp )   ! 1D config.: same bathy value over the 3x3 domain
124      IF( ln_zco      )   CALL zgr_zco          ! z-coordinate
125      IF( ln_zps      )   CALL zgr_zps          ! Partial step z-coordinate
126      IF( ln_sco.AND. .NOT. ln_hyb )   CALL zgr_sco          ! s-coordinate or hybrid z-s coordinate (z at upper levels )
127      IF( ln_sco  .AND.     ln_hyb )   CALL zgr_hyb          ! hybrid s-sigma z      ( s- at shel
128      !
129      ! final adjustment of mbathy & check
130      ! -----------------------------------
131      IF( lzoom       )   CALL zgr_bat_zoom     ! correct mbathy in case of zoom subdomain
132      IF( .NOT.lk_c1d )   CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isolated ocean points
133                          CALL zgr_bot_level    ! deepest ocean level for t-, u- and v-points
134      !
135      IF( lk_c1d ) THEN                         ! 1D config.: same mbathy value over the 3x3 domain
136         ibat = mbathy(2,2)
137         mbathy(:,:) = ibat
138      END IF
139      !
140      IF( nprint == 1 .AND. lwp )   THEN
141         WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
142         WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ),   &
143            &                   ' w ',   MINVAL( fsdepw(:,:,:) ), '3w ', MINVAL( fsde3w(:,:,:) )
144         WRITE(numout,*) ' MIN val e3    t ', MINVAL( fse3t(:,:,:) ), ' f ', MINVAL( fse3f(:,:,:) ),  &
145            &                   ' u ',   MINVAL( fse3u(:,:,:) ), ' u ', MINVAL( fse3v(:,:,:) ),  &
146            &                   ' uw',   MINVAL( fse3uw(:,:,:)), ' vw', MINVAL( fse3vw(:,:,:)),   &
147            &                   ' w ',   MINVAL( fse3w(:,:,:) )
148
149         WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ),   &
150            &                   ' w ',   MAXVAL( fsdepw(:,:,:) ), '3w ', MAXVAL( fsde3w(:,:,:) )
151         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( fse3t(:,:,:) ), ' f ', MAXVAL( fse3f(:,:,:) ),  &
152            &                   ' u ',   MAXVAL( fse3u(:,:,:) ), ' u ', MAXVAL( fse3v(:,:,:) ),  &
153            &                   ' uw',   MAXVAL( fse3uw(:,:,:)), ' vw', MAXVAL( fse3vw(:,:,:)),   &
154            &                   ' w ',   MAXVAL( fse3w(:,:,:) )
155      ENDIF
156      !
157      IF( nn_timing == 1 )  CALL timing_stop('dom_zgr')
158      !
159   END SUBROUTINE dom_zgr
160
161
162   SUBROUTINE zgr_z
163      !!----------------------------------------------------------------------
164      !!                   ***  ROUTINE zgr_z  ***
165      !!                   
166      !! ** Purpose :   set the depth of model levels and the resulting
167      !!      vertical scale factors.
168      !!
169      !! ** Method  :   z-coordinate system (use in all type of coordinate)
170      !!        The depth of model levels is defined from an analytical
171      !!      function the derivative of which gives the scale factors.
172      !!        both depth and scale factors only depend on k (1d arrays).
173      !!              w-level: gdepw_0  = fsdep(k)
174      !!                       e3w_0(k) = dk(fsdep)(k)     = fse3(k)
175      !!              t-level: gdept_0  = fsdep(k+0.5)
176      !!                       e3t_0(k) = dk(fsdep)(k+0.5) = fse3(k+0.5)
177      !!
178      !! ** Action  : - gdept_0, gdepw_0 : depth of T- and W-point (m)
179      !!              - e3t_0  , e3w_0   : scale factors at T- and W-levels (m)
180      !!
181      !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
182      !!----------------------------------------------------------------------
183      INTEGER  ::   jk                     ! dummy loop indices
184      REAL(wp) ::   zt, zw                 ! temporary scalars
185      REAL(wp) ::   zsur, za0, za1, zkth   ! Values set from parameters in
186      REAL(wp) ::   zacr, zdzmin, zhmax    ! par_CONFIG_Rxx.h90
187      REAL(wp) ::   zrefdep                ! depth of the reference level (~10m)
188      REAL(wp) ::   za2, zkth2, zacr2      ! Values for optional double tanh function set from parameters
189      !!----------------------------------------------------------------------
190      !
191      IF( nn_timing == 1 )  CALL timing_start('zgr_z')
192      !
193      ! Set variables from parameters
194      ! ------------------------------
195       zkth = ppkth       ;   zacr = ppacr
196       zdzmin = ppdzmin   ;   zhmax = pphmax
197       zkth2 = ppkth2     ;   zacr2 = ppacr2   ! optional (ldbletanh=T) double tanh parameters
198
199      ! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
200      !  za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
201      IF(   ppa1  == pp_to_be_computed  .AND.  &
202         &  ppa0  == pp_to_be_computed  .AND.  &
203         &  ppsur == pp_to_be_computed           ) THEN
204         !
205         za1  = (  ppdzmin - pphmax / FLOAT(jpkm1)  )                                                      &
206            & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * (  LOG( COSH( (jpk - ppkth) / ppacr) )      &
207            &                                                   - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
208         za0  = ppdzmin - za1 *              TANH( (1-ppkth) / ppacr )
209         zsur =   - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
210      ELSE
211         za1 = ppa1 ;       za0 = ppa0 ;          zsur = ppsur
212         za2 = ppa2                            ! optional (ldbletanh=T) double tanh parameter
213      ENDIF
214
215      IF(lwp) THEN                         ! Parameter print
216         WRITE(numout,*)
217         WRITE(numout,*) '    zgr_z   : Reference vertical z-coordinates'
218         WRITE(numout,*) '    ~~~~~~~'
219         IF(  ppkth == 0._wp ) THEN             
220              WRITE(numout,*) '            Uniform grid with ',jpk-1,' layers'
221              WRITE(numout,*) '            Total depth    :', zhmax
222              WRITE(numout,*) '            Layer thickness:', zhmax/(jpk-1)
223         ELSE
224            IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN
225               WRITE(numout,*) '         zsur, za0, za1 computed from '
226               WRITE(numout,*) '                 zdzmin = ', zdzmin
227               WRITE(numout,*) '                 zhmax  = ', zhmax
228            ENDIF
229            WRITE(numout,*) '           Value of coefficients for vertical mesh:'
230            WRITE(numout,*) '                 zsur = ', zsur
231            WRITE(numout,*) '                 za0  = ', za0
232            WRITE(numout,*) '                 za1  = ', za1
233            WRITE(numout,*) '                 zkth = ', zkth
234            WRITE(numout,*) '                 zacr = ', zacr
235            IF( ldbletanh ) THEN
236               WRITE(numout,*) ' (Double tanh    za2  = ', za2
237               WRITE(numout,*) '  parameters)    zkth2= ', zkth2
238               WRITE(numout,*) '                 zacr2= ', zacr2
239            ENDIF
240         ENDIF
241      ENDIF
242
243
244      ! Reference z-coordinate (depth - scale factor at T- and W-points)
245      ! ======================
246      IF( ppkth == 0._wp ) THEN            !  uniform vertical grid       
247         za1 = zhmax / FLOAT(jpk-1) 
248         DO jk = 1, jpk
249            zw = FLOAT( jk )
250            zt = FLOAT( jk ) + 0.5_wp
251            gdepw_0(jk) = ( zw - 1 ) * za1
252            gdept_0(jk) = ( zt - 1 ) * za1
253            e3w_0  (jk) =  za1
254            e3t_0  (jk) =  za1
255         END DO
256      ELSE                                ! Madec & Imbard 1996 function
257         IF( .NOT. ldbletanh ) THEN
258            DO jk = 1, jpk
259               zw = REAL( jk , wp )
260               zt = REAL( jk , wp ) + 0.5_wp
261               gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) )  )
262               gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) )  )
263               e3w_0  (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   )
264               e3t_0  (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   )
265            END DO
266         ELSE
267            DO jk = 1, jpk
268               zw = FLOAT( jk )
269               zt = FLOAT( jk ) + 0.5_wp
270               ! Double tanh function
271               gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr  ) )    &
272                  &                            + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) )  )
273               gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr  ) )    &
274                  &                            + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) )  )
275               e3w_0  (jk) =          za0      + za1        * TANH(       (zw-zkth ) / zacr  )    &
276                  &                            + za2        * TANH(       (zw-zkth2) / zacr2 )
277               e3t_0  (jk) =          za0      + za1        * TANH(       (zt-zkth ) / zacr  )    &
278                  &                            + za2        * TANH(       (zt-zkth2) / zacr2 )
279            END DO
280         ENDIF
281         gdepw_0(1) = 0._wp                    ! force first w-level to be exactly at zero
282      ENDIF
283
284!!gm BUG in s-coordinate this does not work!
285      ! deepest/shallowest W level Above/Below ~10m
286      zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_0 )                    ! ref. depth with tolerance (10% of minimum layer thickness)
287      nlb10 = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 )   ! shallowest W level Below ~10m
288      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m
289!!gm end bug
290
291      IF(lwp) THEN                        ! control print
292         WRITE(numout,*)
293         WRITE(numout,*) '              Reference z-coordinate depth and scale factors:'
294         WRITE(numout, "(9x,' level   gdept    gdepw     e3t      e3w  ')" )
295         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_0(jk), gdepw_0(jk), e3t_0(jk), e3w_0(jk), jk = 1, jpk )
296      ENDIF
297      DO jk = 1, jpk                      ! control positivity
298         IF( e3w_0  (jk) <= 0._wp .OR. e3t_0  (jk) <= 0._wp )   CALL ctl_stop( 'dom:zgr_z: e3w or e3t =< 0 '    )
299         IF( gdepw_0(jk) <  0._wp .OR. gdept_0(jk) <  0._wp )   CALL ctl_stop( 'dom:zgr_z: gdepw or gdept < 0 ' )
300      END DO
301      !
302      IF( nn_timing == 1 )  CALL timing_stop('zgr_z')
303      !
304   END SUBROUTINE zgr_z
305
306
307   SUBROUTINE zgr_bat
308      !!----------------------------------------------------------------------
309      !!                    ***  ROUTINE zgr_bat  ***
310      !!
311      !! ** Purpose :   set bathymetry both in levels and meters
312      !!
313      !! ** Method  :   read or define mbathy and bathy arrays
314      !!       * level bathymetry:
315      !!      The ocean basin geometry is given by a two-dimensional array,
316      !!      mbathy, which is defined as follow :
317      !!            mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
318      !!                              at t-point (ji,jj).
319      !!                            = 0  over the continental t-point.
320      !!      The array mbathy is checked to verified its consistency with
321      !!      model option. in particular:
322      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
323      !!                  along closed boundary.
324      !!            mbathy must be cyclic IF jperio=1.
325      !!            mbathy must be lower or equal to jpk-1.
326      !!            isolated ocean grid points are suppressed from mbathy
327      !!                  since they are only connected to remaining
328      !!                  ocean through vertical diffusion.
329      !!      ntopo=-1 :   rectangular channel or bassin with a bump
330      !!      ntopo= 0 :   flat rectangular channel or basin
331      !!      ntopo= 1 :   mbathy is read in 'bathy_level.nc' NetCDF file
332      !!                   bathy  is read in 'bathy_meter.nc' NetCDF file
333      !!
334      !! ** Action  : - mbathy: level bathymetry (in level index)
335      !!              - bathy : meter bathymetry (in meters)
336      !!----------------------------------------------------------------------
337      INTEGER  ::   ji, jj, jl, jk            ! dummy loop indices
338      INTEGER  ::   inum                      ! temporary logical unit
339      INTEGER  ::   ii_bump, ij_bump, ih      ! bump center position
340      INTEGER  ::   ii0, ii1, ij0, ij1, ik    ! local indices
341      REAL(wp) ::   r_bump , h_bump , h_oce   ! bump characteristics
342      REAL(wp) ::   zi, zj, zh, zhmin         ! local scalars
343      INTEGER , POINTER, DIMENSION(:,:) ::   idta   ! global domain integer data
344      REAL(wp), POINTER, DIMENSION(:,:) ::   zdta   ! global domain scalar data
345      !!----------------------------------------------------------------------
346      !
347      IF( nn_timing == 1 )  CALL timing_start('zgr_bat')
348      !
349      CALL wrk_alloc( jpidta, jpjdta, idta )
350      CALL wrk_alloc( jpidta, jpjdta, zdta )
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 == 2 ) THEN    ! ORCA R2 configuration
470               !
471              IF( nn_cla == 0 ) THEN
472                 ii0 = 140   ;   ii1 = 140                   ! Gibraltar Strait open
473                 ij0 = 102   ;   ij1 = 102                   ! (Thomson, Ocean Modelling, 1995)
474                 DO ji = mi0(ii0), mi1(ii1)
475                    DO jj = mj0(ij0), mj1(ij1)
476                       bathy(ji,jj) = 284._wp
477                    END DO
478                 END DO
479                 IF(lwp) WRITE(numout,*)     
480                 IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
481                 !
482                 ii0 = 160   ;   ii1 = 160                   ! Bab el mandeb Strait open
483                 ij0 = 88    ;   ij1 = 88                    ! (Thomson, Ocean Modelling, 1995)
484                 DO ji = mi0(ii0), mi1(ii1)
485                    DO jj = mj0(ij0), mj1(ij1)
486                       bathy(ji,jj) = 137._wp
487                    END DO
488                 END DO
489                 IF(lwp) WRITE(numout,*)
490                 IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
491              ENDIF
492              !
493           ENDIF
494            !
495        ENDIF
496         !                                            ! =============== !
497      ELSE                                            !      error      !
498         !                                            ! =============== !
499         WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
500         CALL ctl_stop( '    zgr_bat : '//trim(ctmp1) )
501      ENDIF
502      !
503      IF( nn_closea == 0 )   CALL clo_bat( bathy, mbathy )    !==  NO closed seas or lakes  ==!
504      !                       
505      IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==!
506         IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level
507         ELSE                          ;   ik = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
508         ENDIF
509         zhmin = gdepw_0(ik+1)                                                         ! minimum depth = ik+1 w-levels
510         WHERE( bathy(:,:) <= 0._wp )   ;   bathy(:,:) = 0._wp                         ! min=0     over the lands
511         ELSE WHERE                     ;   bathy(:,:) = MAX(  zhmin , bathy(:,:)  )   ! min=zhmin over the oceans
512         END WHERE
513         IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
514      ENDIF
515      !
516      !
517      CALL wrk_dealloc( jpidta, jpjdta, idta )
518      CALL wrk_dealloc( jpidta, jpjdta, zdta )
519      !
520      IF( nn_timing == 1 )  CALL timing_stop('zgr_bat')
521      !
522   END SUBROUTINE zgr_bat
523
524
525   SUBROUTINE zgr_bat_zoom
526      !!----------------------------------------------------------------------
527      !!                    ***  ROUTINE zgr_bat_zoom  ***
528      !!
529      !! ** Purpose : - Close zoom domain boundary if necessary
530      !!              - Suppress Med Sea from ORCA R2 and R05 arctic zoom
531      !!
532      !! ** Method  :
533      !!
534      !! ** Action  : - update mbathy: level bathymetry (in level index)
535      !!----------------------------------------------------------------------
536      INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integers
537      !!----------------------------------------------------------------------
538      !
539      IF(lwp) WRITE(numout,*)
540      IF(lwp) WRITE(numout,*) '    zgr_bat_zoom : modify the level bathymetry for zoom domain'
541      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~'
542      !
543      ! Zoom domain
544      ! ===========
545      !
546      ! Forced closed boundary if required
547      IF( lzoom_s )   mbathy(  : , mj0(jpjzoom):mj1(jpjzoom) )      = 0
548      IF( lzoom_w )   mbathy(      mi0(jpizoom):mi1(jpizoom) , :  ) = 0
549      IF( lzoom_e )   mbathy(      mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , :  ) = 0
550      IF( lzoom_n )   mbathy(  : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) )      = 0
551      !
552      ! Configuration specific domain modifications
553      ! (here, ORCA arctic configuration: suppress Med Sea)
554      IF( cp_cfg == "orca" .AND. lzoom_arct ) THEN
555         SELECT CASE ( jp_cfg )
556         !                                        ! =======================
557         CASE ( 2 )                               !  ORCA_R2 configuration
558            !                                     ! =======================
559            IF(lwp) WRITE(numout,*) '                   ORCA R2 arctic zoom: suppress the Med Sea'
560            ii0 = 141   ;   ii1 = 162      ! Sea box i,j indices
561            ij0 =  98   ;   ij1 = 110
562            !                                     ! =======================
563         CASE ( 05 )                              !  ORCA_R05 configuration
564            !                                     ! =======================
565            IF(lwp) WRITE(numout,*) '                   ORCA R05 arctic zoom: suppress the Med Sea'
566            ii0 = 563   ;   ii1 = 642      ! zero over the Med Sea boxe
567            ij0 = 314   ;   ij1 = 370 
568         END SELECT
569         !
570         mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0   ! zero over the Med Sea boxe
571         !
572      ENDIF
573      !
574   END SUBROUTINE zgr_bat_zoom
575
576
577   SUBROUTINE zgr_bat_ctl
578      !!----------------------------------------------------------------------
579      !!                    ***  ROUTINE zgr_bat_ctl  ***
580      !!
581      !! ** Purpose :   check the bathymetry in levels
582      !!
583      !! ** Method  :   The array mbathy is checked to verified its consistency
584      !!      with the model options. in particular:
585      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
586      !!                  along closed boundary.
587      !!            mbathy must be cyclic IF jperio=1.
588      !!            mbathy must be lower or equal to jpk-1.
589      !!            isolated ocean grid points are suppressed from mbathy
590      !!                  since they are only connected to remaining
591      !!                  ocean through vertical diffusion.
592      !!      C A U T I O N : mbathy will be modified during the initializa-
593      !!      tion phase to become the number of non-zero w-levels of a water
594      !!      column, with a minimum value of 1.
595      !!
596      !! ** Action  : - update mbathy: level bathymetry (in level index)
597      !!              - update bathy : meter bathymetry (in meters)
598      !!----------------------------------------------------------------------
599      !!
600      INTEGER ::   ji, jj, jl                    ! dummy loop indices
601      INTEGER ::   icompt, ibtest, ikmax         ! temporary integers
602      REAL(wp), POINTER, DIMENSION(:,:) ::  zbathy
603      !!----------------------------------------------------------------------
604      !
605      IF( nn_timing == 1 )  CALL timing_start('zgr_bat_ctl')
606      !
607      CALL wrk_alloc( jpi, jpj, zbathy )
608      !
609      IF(lwp) WRITE(numout,*)
610      IF(lwp) WRITE(numout,*) '    zgr_bat_ctl : check the bathymetry'
611      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
612
613      !                                          ! Suppress isolated ocean grid points
614      IF(lwp) WRITE(numout,*)
615      IF(lwp) WRITE(numout,*)'                   suppress isolated ocean grid points'
616      IF(lwp) WRITE(numout,*)'                   -----------------------------------'
617      icompt = 0
618      DO jl = 1, 2
619         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
620            mbathy( 1 ,:) = mbathy(jpim1,:)           ! local domain is cyclic east-west
621            mbathy(jpi,:) = mbathy(  2  ,:)
622         ENDIF
623         DO jj = 2, jpjm1
624            DO ji = 2, jpim1
625               ibtest = MAX(  mbathy(ji-1,jj), mbathy(ji+1,jj),   &
626                  &           mbathy(ji,jj-1), mbathy(ji,jj+1)  )
627               IF( ibtest < mbathy(ji,jj) ) THEN
628                  IF(lwp) WRITE(numout,*) ' the number of ocean level at ',   &
629                     &   'grid-point (i,j) =  ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
630                  mbathy(ji,jj) = ibtest
631                  icompt = icompt + 1
632               ENDIF
633            END DO
634         END DO
635      END DO
636      IF( lk_mpp )   CALL mpp_sum( icompt )
637      IF( icompt == 0 ) THEN
638         IF(lwp) WRITE(numout,*)'     no isolated ocean grid points'
639      ELSE
640         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed'
641      ENDIF
642      IF( lk_mpp ) THEN
643         zbathy(:,:) = FLOAT( mbathy(:,:) )
644         CALL lbc_lnk( zbathy, 'T', 1._wp )
645         mbathy(:,:) = INT( zbathy(:,:) )
646      ENDIF
647
648      !                                          ! East-west cyclic boundary conditions
649      IF( nperio == 0 ) THEN
650         IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio
651         IF( lk_mpp ) THEN
652            IF( nbondi == -1 .OR. nbondi == 2 ) THEN
653               IF( jperio /= 1 )   mbathy(1,:) = 0
654            ENDIF
655            IF( nbondi == 1 .OR. nbondi == 2 ) THEN
656               IF( jperio /= 1 )   mbathy(nlci,:) = 0
657            ENDIF
658         ELSE
659            IF( ln_zco .OR. ln_zps ) THEN
660               mbathy( 1 ,:) = 0
661               mbathy(jpi,:) = 0
662            ELSE
663               mbathy( 1 ,:) = jpkm1
664               mbathy(jpi,:) = jpkm1
665            ENDIF
666         ENDIF
667      ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio ==  6 ) THEN
668         IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio
669         mbathy( 1 ,:) = mbathy(jpim1,:)
670         mbathy(jpi,:) = mbathy(  2  ,:)
671      ELSEIF( nperio == 2 ) THEN
672         IF(lwp) WRITE(numout,*) '   equatorial boundary conditions on mbathy: nperio = ', nperio
673      ELSE
674         IF(lwp) WRITE(numout,*) '    e r r o r'
675         IF(lwp) WRITE(numout,*) '    parameter , nperio = ', nperio
676         !         STOP 'dom_mba'
677      ENDIF
678
679      !  Boundary condition on mbathy
680      IF( .NOT.lk_mpp ) THEN 
681!!gm     !!bug ???  think about it !
682         !   ... mono- or macro-tasking: T-point, >0, 2D array, no slab
683         zbathy(:,:) = FLOAT( mbathy(:,:) )
684         CALL lbc_lnk( zbathy, 'T', 1._wp )
685         mbathy(:,:) = INT( zbathy(:,:) )
686      ENDIF
687
688      ! Number of ocean level inferior or equal to jpkm1
689      ikmax = 0
690      DO jj = 1, jpj
691         DO ji = 1, jpi
692            ikmax = MAX( ikmax, mbathy(ji,jj) )
693         END DO
694      END DO
695!!gm  !!! test to do:   ikmax = MAX( mbathy(:,:) )   ???
696      IF( ikmax > jpkm1 ) THEN
697         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' >  jpk-1'
698         IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
699      ELSE IF( ikmax < jpkm1 ) THEN
700         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1' 
701         IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
702      ENDIF
703
704      IF( lwp .AND. nprint == 1 ) THEN      ! control print
705         WRITE(numout,*)
706         WRITE(numout,*) ' bathymetric field :   number of non-zero T-levels '
707         WRITE(numout,*) ' ------------------'
708         CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )
709         WRITE(numout,*)
710      ENDIF
711      !
712      CALL wrk_dealloc( jpi, jpj, zbathy )
713      !
714      IF( nn_timing == 1 )  CALL timing_stop('zgr_bat_ctl')
715      !
716   END SUBROUTINE zgr_bat_ctl
717
718
719   SUBROUTINE zgr_bot_level
720      !!----------------------------------------------------------------------
721      !!                    ***  ROUTINE zgr_bot_level  ***
722      !!
723      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
724      !!
725      !! ** Method  :   computes from mbathy with a minimum value of 1 over land
726      !!
727      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest
728      !!                                     ocean level at t-, u- & v-points
729      !!                                     (min value = 1 over land)
730      !!----------------------------------------------------------------------
731      !!
732      INTEGER ::   ji, jj   ! dummy loop indices
733      REAL(wp), POINTER, DIMENSION(:,:) ::  zmbk
734      !!----------------------------------------------------------------------
735      !
736      IF( nn_timing == 1 )  CALL timing_start('zgr_bot_level')
737      !
738      CALL wrk_alloc( jpi, jpj, zmbk )
739      !
740      IF(lwp) WRITE(numout,*)
741      IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
742      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
743      !
744      mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
745 
746      !                                     ! bottom k-index of W-level = mbkt+1
747      DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level
748         DO ji = 1, jpim1
749            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
750            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
751         END DO
752      END DO
753      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
754      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
755      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
756      !
757      CALL wrk_dealloc( jpi, jpj, zmbk )
758      !
759      IF( nn_timing == 1 )  CALL timing_stop('zgr_bot_level')
760      !
761   END SUBROUTINE zgr_bot_level
762
763
764   SUBROUTINE zgr_zco
765      !!----------------------------------------------------------------------
766      !!                  ***  ROUTINE zgr_zco  ***
767      !!
768      !! ** Purpose :   define the z-coordinate system
769      !!
770      !! ** Method  :   set 3D coord. arrays to reference 1D array
771      !!----------------------------------------------------------------------
772      INTEGER  ::   jk
773      !!----------------------------------------------------------------------
774      !
775      IF( nn_timing == 1 )  CALL timing_start('zgr_zco')
776      !
777      DO jk = 1, jpk
778            gdept(:,:,jk) = gdept_0(jk)
779            gdepw(:,:,jk) = gdepw_0(jk)
780            gdep3w(:,:,jk) = gdepw_0(jk)
781            e3t (:,:,jk) = e3t_0(jk)
782            e3u (:,:,jk) = e3t_0(jk)
783            e3v (:,:,jk) = e3t_0(jk)
784            e3f (:,:,jk) = e3t_0(jk)
785            e3w (:,:,jk) = e3w_0(jk)
786            e3uw(:,:,jk) = e3w_0(jk)
787            e3vw(:,:,jk) = e3w_0(jk)
788      END DO
789      !
790      IF( nn_timing == 1 )  CALL timing_stop('zgr_zco')
791      !
792   END SUBROUTINE zgr_zco
793
794
795   SUBROUTINE zgr_zps
796      !!----------------------------------------------------------------------
797      !!                  ***  ROUTINE zgr_zps  ***
798      !!                     
799      !! ** Purpose :   the depth and vertical scale factor in partial step
800      !!      z-coordinate case
801      !!
802      !! ** Method  :   Partial steps : computes the 3D vertical scale factors
803      !!      of T-, U-, V-, W-, UW-, VW and F-points that are associated with
804      !!      a partial step representation of bottom topography.
805      !!
806      !!        The reference depth of model levels is defined from an analytical
807      !!      function the derivative of which gives the reference vertical
808      !!      scale factors.
809      !!        From  depth and scale factors reference, we compute there new value
810      !!      with partial steps  on 3d arrays ( i, j, k ).
811      !!
812      !!              w-level: gdepw(i,j,k)  = fsdep(k)
813      !!                       e3w(i,j,k) = dk(fsdep)(k)     = fse3(i,j,k)
814      !!              t-level: gdept(i,j,k)  = fsdep(k+0.5)
815      !!                       e3t(i,j,k) = dk(fsdep)(k+0.5) = fse3(i,j,k+0.5)
816      !!
817      !!        With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
818      !!      we find the mbathy index of the depth at each grid point.
819      !!      This leads us to three cases:
820      !!
821      !!              - bathy = 0 => mbathy = 0
822      !!              - 1 < mbathy < jpkm1   
823      !!              - bathy > gdepw(jpk) => mbathy = jpkm1 
824      !!
825      !!        Then, for each case, we find the new depth at t- and w- levels
826      !!      and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
827      !!      and f-points.
828      !!
829      !!        This routine is given as an example, it must be modified
830      !!      following the user s desiderata. nevertheless, the output as
831      !!      well as the way to compute the model levels and scale factors
832      !!      must be respected in order to insure second order accuracy
833      !!      schemes.
834      !!
835      !!         c a u t i o n : gdept_0, gdepw_0 and e3._0 are positives
836      !!         - - - - - - -   gdept, gdepw and e3. are positives
837      !!     
838      !!  Reference :   Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
839      !!----------------------------------------------------------------------
840      !!
841      INTEGER  ::   ji, jj, jk       ! dummy loop indices
842      INTEGER  ::   ik, it           ! temporary integers
843      LOGICAL  ::   ll_print         ! Allow  control print for debugging
844      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
845      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
846      REAL(wp) ::   zmax             ! Maximum depth
847      REAL(wp) ::   zdiff            ! temporary scalar
848      REAL(wp) ::   zrefdep          ! temporary scalar
849      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt
850      !!---------------------------------------------------------------------
851      !
852      IF( nn_timing == 1 )  CALL timing_start('zgr_zps')
853      !
854      CALL wrk_alloc( jpi, jpj, jpk, zprt )
855      !
856      IF(lwp) WRITE(numout,*)
857      IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps'
858      IF(lwp) WRITE(numout,*) '    ~~~~~~~ '
859      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used'
860
861      ll_print = .FALSE.                   ! Local variable for debugging
862     
863      IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth
864         WRITE(numout,*)
865         WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)'
866         CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
867      ENDIF
868
869
870      ! bathymetry in level (from bathy_meter)
871      ! ===================
872      zmax = gdepw_0(jpk) + e3t_0(jpk)          ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) )
873      bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat)
874      WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0
875      ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level
876      END WHERE
877
878      ! Compute mbathy for ocean points (i.e. the number of ocean levels)
879      ! find the number of ocean levels such that the last level thickness
880      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_0 (where
881      ! e3t_0 is the reference level thickness
882      DO jk = jpkm1, 1, -1
883         zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat )
884         WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1
885      END DO
886
887      ! Scale factors and depth at T- and W-points
888      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate
889         gdept(:,:,jk) = gdept_0(jk)
890         gdepw(:,:,jk) = gdepw_0(jk)
891         e3t  (:,:,jk) = e3t_0  (jk)
892         e3w  (:,:,jk) = e3w_0  (jk)
893      END DO
894      !
895      DO jj = 1, jpj
896         DO ji = 1, jpi
897            ik = mbathy(ji,jj)
898            IF( ik > 0 ) THEN               ! ocean point only
899               ! max ocean level case
900               IF( ik == jpkm1 ) THEN
901                  zdepwp = bathy(ji,jj)
902                  ze3tp  = bathy(ji,jj) - gdepw_0(ik)
903                  ze3wp = 0.5_wp * e3w_0(ik) * ( 1._wp + ( ze3tp/e3t_0(ik) ) )
904                  e3t(ji,jj,ik  ) = ze3tp
905                  e3t(ji,jj,ik+1) = ze3tp
906                  e3w(ji,jj,ik  ) = ze3wp
907                  e3w(ji,jj,ik+1) = ze3tp
908                  gdepw(ji,jj,ik+1) = zdepwp
909                  gdept(ji,jj,ik  ) = gdept_0(ik-1) + ze3wp
910                  gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp
911                  !
912               ELSE                         ! standard case
913                  IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN   ;   gdepw(ji,jj,ik+1) = bathy(ji,jj)
914                  ELSE                                       ;   gdepw(ji,jj,ik+1) = gdepw_0(ik+1)
915                  ENDIF
916!gm Bug?  check the gdepw_0
917                  !       ... on ik
918                  gdept(ji,jj,ik) = gdepw_0(ik) + ( gdepw  (ji,jj,ik+1) - gdepw_0(ik) )   &
919                     &                          * ((gdept_0(      ik  ) - gdepw_0(ik) )   &
920                     &                          / ( gdepw_0(      ik+1) - gdepw_0(ik) ))
921                  e3t  (ji,jj,ik) = e3t_0  (ik) * ( gdepw  (ji,jj,ik+1) - gdepw_0(ik) )   & 
922                     &                          / ( gdepw_0(      ik+1) - gdepw_0(ik) ) 
923                  e3w  (ji,jj,ik) = 0.5_wp * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2._wp * gdepw_0(ik) )   &
924                     &                     * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) )
925                  !       ... on ik+1
926                  e3w  (ji,jj,ik+1) = e3t  (ji,jj,ik)
927                  e3t  (ji,jj,ik+1) = e3t  (ji,jj,ik)
928                  gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)
929               ENDIF
930            ENDIF
931         END DO
932      END DO
933      !
934      it = 0
935      DO jj = 1, jpj
936         DO ji = 1, jpi
937            ik = mbathy(ji,jj)
938            IF( ik > 0 ) THEN               ! ocean point only
939               e3tp (ji,jj) = e3t(ji,jj,ik  )
940               e3wp (ji,jj) = e3w(ji,jj,ik  )
941               ! test
942               zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik  )
943               IF( zdiff <= 0._wp .AND. lwp ) THEN
944                  it = it + 1
945                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj
946                  WRITE(numout,*) ' bathy = ', bathy(ji,jj)
947                  WRITE(numout,*) ' gdept = ', gdept(ji,jj,ik), ' gdepw = ', gdepw(ji,jj,ik+1), ' zdiff = ', zdiff
948                  WRITE(numout,*) ' e3tp  = ', e3t  (ji,jj,ik), ' e3wp  = ', e3w  (ji,jj,ik  )
949               ENDIF
950            ENDIF
951         END DO
952      END DO
953
954      ! Scale factors and depth at U-, V-, UW and VW-points
955      DO jk = 1, jpk                        ! initialisation to z-scale factors
956         e3u (:,:,jk) = e3t_0(jk)
957         e3v (:,:,jk) = e3t_0(jk)
958         e3uw(:,:,jk) = e3w_0(jk)
959         e3vw(:,:,jk) = e3w_0(jk)
960      END DO
961      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
962         DO jj = 1, jpjm1
963            DO ji = 1, fs_jpim1   ! vector opt.
964               e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk) )
965               e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk) )
966               e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) )
967               e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) )
968            END DO
969         END DO
970      END DO
971      CALL lbc_lnk( e3u , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw, 'U', 1._wp )   ! lateral boundary conditions
972      CALL lbc_lnk( e3v , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw, 'V', 1._wp )
973      !
974      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
975         WHERE( e3u (:,:,jk) == 0._wp )   e3u (:,:,jk) = e3t_0(jk)
976         WHERE( e3v (:,:,jk) == 0._wp )   e3v (:,:,jk) = e3t_0(jk)
977         WHERE( e3uw(:,:,jk) == 0._wp )   e3uw(:,:,jk) = e3w_0(jk)
978         WHERE( e3vw(:,:,jk) == 0._wp )   e3vw(:,:,jk) = e3w_0(jk)
979      END DO
980     
981      ! Scale factor at F-point
982      DO jk = 1, jpk                        ! initialisation to z-scale factors
983         e3f(:,:,jk) = e3t_0(jk)
984      END DO
985      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors
986         DO jj = 1, jpjm1
987            DO ji = 1, fs_jpim1   ! vector opt.
988               e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) )
989            END DO
990         END DO
991      END DO
992      CALL lbc_lnk( e3f, 'F', 1._wp )       ! Lateral boundary conditions
993      !
994      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
995         WHERE( e3f(:,:,jk) == 0._wp )   e3f(:,:,jk) = e3t_0(jk)
996      END DO
997!!gm  bug ? :  must be a do loop with mj0,mj1
998      !
999      e3t(:,mj0(1),:) = e3t(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2
1000      e3w(:,mj0(1),:) = e3w(:,mj0(2),:) 
1001      e3u(:,mj0(1),:) = e3u(:,mj0(2),:) 
1002      e3v(:,mj0(1),:) = e3v(:,mj0(2),:) 
1003      e3f(:,mj0(1),:) = e3f(:,mj0(2),:) 
1004
1005      ! Control of the sign
1006      IF( MINVAL( e3t  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t   <= 0' )
1007      IF( MINVAL( e3w  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w   <= 0' )
1008      IF( MINVAL( gdept(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw <  0' )
1009      IF( MINVAL( gdepw(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw <  0' )
1010     
1011      ! Compute gdep3w (vertical sum of e3w)
1012      gdep3w(:,:,1) = 0.5_wp * e3w(:,:,1)
1013      DO jk = 2, jpk
1014         gdep3w(:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk) 
1015      END DO
1016       
1017      !                                               ! ================= !
1018      IF(lwp .AND. ll_print) THEN                     !   Control print   !
1019         !                                            ! ================= !
1020         DO jj = 1,jpj
1021            DO ji = 1, jpi
1022               ik = MAX( mbathy(ji,jj), 1 )
1023               zprt(ji,jj,1) = e3t   (ji,jj,ik)
1024               zprt(ji,jj,2) = e3w   (ji,jj,ik)
1025               zprt(ji,jj,3) = e3u   (ji,jj,ik)
1026               zprt(ji,jj,4) = e3v   (ji,jj,ik)
1027               zprt(ji,jj,5) = e3f   (ji,jj,ik)
1028               zprt(ji,jj,6) = gdep3w(ji,jj,ik)
1029            END DO
1030         END DO
1031         WRITE(numout,*)
1032         WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1033         WRITE(numout,*)
1034         WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1035         WRITE(numout,*)
1036         WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1037         WRITE(numout,*)
1038         WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1039         WRITE(numout,*)
1040         WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1041         WRITE(numout,*)
1042         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1043      ENDIF 
1044      !
1045      CALL wrk_dealloc( jpi, jpj, jpk, zprt )
1046      !
1047      IF( nn_timing == 1 )  CALL timing_stop('zgr_zps')
1048      !
1049   END SUBROUTINE zgr_zps
1050
1051
1052   FUNCTION fssig( pk ) RESULT( pf )
1053      !!----------------------------------------------------------------------
1054      !!                 ***  ROUTINE eos_init  ***
1055      !!       
1056      !! ** Purpose :   provide the analytical function in s-coordinate
1057      !!         
1058      !! ** Method  :   the function provide the non-dimensional position of
1059      !!                T and W (i.e. between 0 and 1)
1060      !!                T-points at integer values (between 1 and jpk)
1061      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1062      !!----------------------------------------------------------------------
1063      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
1064      REAL(wp)             ::   pf   ! sigma value
1065      !!----------------------------------------------------------------------
1066      !
1067      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   &
1068         &     - TANH( rn_thetb * rn_theta                                )  )   &
1069         & * (   COSH( rn_theta                           )                      &
1070         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              &
1071         & / ( 2._wp * SINH( rn_theta ) )
1072      !
1073   END FUNCTION fssig
1074
1075
1076   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
1077      !!----------------------------------------------------------------------
1078      !!                 ***  ROUTINE eos_init  ***
1079      !!
1080      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
1081      !!
1082      !! ** Method  :   the function provides the non-dimensional position of
1083      !!                T and W (i.e. between 0 and 1)
1084      !!                T-points at integer values (between 1 and jpk)
1085      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1086      !!----------------------------------------------------------------------
1087      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
1088      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
1089      REAL(wp)             ::   pf1   ! sigma value
1090      !!----------------------------------------------------------------------
1091      !
1092      IF ( rn_theta == 0._wp ) then      ! uniform sigma
1093         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
1094      ELSE                        ! stretched sigma
1095         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
1096            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
1097            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
1098      ENDIF
1099      !
1100   END FUNCTION fssig1
1101   FUNCTION fssig2 ( pk, kmax ) RESULT( pf2 )
1102      !!----------------------------------------------------------------------
1103      !!                 ***  ROUTINE eos_init  ***
1104      !!       
1105      !! ** Purpose :   provide the analytical function in s-coordinate
1106      !!         
1107      !! ** Method  :   the function provide the non-dimensional position of
1108      !!                T and W (i.e. between 0 and 1)
1109      !!                T-points at integer values (between 1 and kmax )
1110      !!                W-points at integer values - 1/2 (between 0.5 and kmax-0.5)
1111      !!
1112      !! Reference  :   ???
1113      !!----------------------------------------------------------------------
1114      REAL(wp), INTENT(in   ) ::   pk   ! continuous "k" coordinate
1115      REAL(wp)                ::   pf2   ! sigma value
1116      INTEGER, INTENT (in)    ::   kmax ! max of sigma)level
1117      !!----------------------------------------------------------------------
1118      !
1119      pf2 =   (   TANH( rn_theta * ( -(pk-0.5) / REAL(kmax-1,wp) + rn_thetb )  )      &
1120         &     - TANH( rn_thetb * rn_theta                                )  )   &
1121         & * (   COSH( rn_theta                           )                   &
1122         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )                &
1123         & / ( 2._wp * SINH( rn_theta ) )
1124      !
1125   END FUNCTION fssig2
1126
1127   FUNCTION fssig3( pk1, pbb ,kmax ) RESULT( pf3 )
1128      !!----------------------------------------------------------------------
1129      !!                 ***  ROUTINE eos_init  ***
1130      !!
1131      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
1132      !!
1133      !! ** Method  :   the function provides the non-dimensional position of
1134      !!                T and W (i.e. between 0 and 1)
1135      !!                T-points at integer values (between 1 and jpk)
1136      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1137      !!
1138      !! Reference  :   ???
1139      !!----------------------------------------------------------------------
1140      REAL(wp), INTENT(in   ) ::   pk1   ! continuous "k" coordinate
1141      REAL(wp), INTENT(in   ) ::   pbb   ! Stretching coefficient
1142      REAL(wp)                ::   pf3   ! sigma value
1143      INTEGER, INTENT (in)    ::   kmax  ! max number of s-sigma levels
1144      !!----------------------------------------------------------------------
1145      !
1146      IF ( rn_theta == 0 ) then      ! uniform sigma
1147         pf3 = -(pk1-0.5_wp) / REAL( kmax-1,wp )
1148      ELSE                        ! stretched sigma
1149         pf3 =   (1.0-pbb) * (sinh( rn_theta*(-(pk1-0.5_wp)/REAL(kmax-1,wp)) ) ) / sinh(rn_theta) + &
1150            &    pbb * ( (tanh( rn_theta*( (-(pk1-0.5_wp)/REAL(kmax-1,wp)) + 0.5_wp) ) - tanh(0.5*rn_theta) ) / &
1151            &    (2._wp*tanh(0.5_wp*rn_theta) ) )
1152      ENDIF
1153   END FUNCTION fssig3
1154   
1155    SUBROUTINE fszref (zkth, zdzmin, zacr, zhmax,jpup,zhsigm )
1156      INTEGER  ::   jk                             ! dummy loop indices
1157      REAL(wp) ::   zt, zw                         ! temporary scalars
1158      REAL(wp) ::   zsur, za0, za1, zkth           ! Values set from parameters in
1159      REAL(wp) ::   zacr, zdzmin, zhmax, zhmax_r    ! read from namelist or par_XXX.h90
1160      REAL(wp) ::   zhsigm                         ! depth of sigma layer
1161      INTEGER  ::   jpup, jpkmax                   ! the last sigma level and number of z-levels
1162      !!----------------------------------------------------------------------
1163! compute reference  depth  leveles 
1164      ! Set variables from parameters
1165      ! ------------------------------
1166      ! zkth = rn_kth       ;   zacr = rn_acr
1167      ! zdzmin = rn_dzmin   ;   zhmax_r = rn_hmax
1168
1169      !  za0, za1, zsur are computed from zdzmin , zhmax, zkth, zacr
1170      !
1171       jpkmax= jpk - jpup
1172       zhmax_r = zhmax - zhsigm
1173
1174         za1  = (  zdzmin - zhmax_r / REAL(jpkmax,wp)  )                                                      &
1175            & / ( TANH((1-zkth)/zacr) - zacr/REAL(jpkmax,wp) * (  LOG( COSH( (jpkmax + 1 - zkth) / zacr) )      &
1176            &                                                   - LOG( COSH( ( 1  - zkth) / zacr) )  )  )
1177         za0  =     zdzmin - za1 *              TANH( (1-zkth) / zacr )
1178         zsur =   - za0 - za1 * zacr * LOG( COSH( (1-zkth) / zacr )  )
1179      ! Reference z-coordinate (depth - scale factor at T- and W-points)
1180      ! ======================
1181      IF( zkth == 0.e0 ) THEN            !  uniform vertical grid
1182         za1 = zhmax_r / REAL(jpkmax-1,wp)
1183         DO jk = 1, jpkmax+1
1184            zw = REAL( jk,wp )
1185            zt = REAL( jk,wp ) + 0.5_wp
1186            gdepw_0(jk+jpup-1 ) = ( zw - 1 ) * za1 + zhsigm
1187            gdept_0(jk+jpup-1 ) = ( zt - 1 ) * za1 + zhsigm
1188            e3w_0  (jk+jpup-1 ) =  za1
1189            e3t_0  (jk+jpup-1 ) =  za1
1190
1191         END DO
1192      ELSE                                ! Madec & Imbard 1996 function
1193         DO jk = 1, jpkmax+1
1194            zw = REAL( jk,wp)
1195            zt = REAL( jk,wp ) + 0.5_wp
1196            gdepw_0(jk+jpup-1) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) )  )+zhsigm
1197            gdept_0(jk+jpup-1) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) )  )+zhsigm
1198            e3w_0  (jk+jpup-1) =          za0      + za1        * TANH(       (zw-zkth) / zacr   )
1199            e3t_0  (jk+jpup-1) =          za0      + za1        * TANH(       (zt-zkth) / zacr   )
1200
1201         END DO
1202         gdepw_0(jpup) = zhsigm                     ! force first w-level to be exactly at zhsigm
1203      ENDIF
1204      IF(lwp) WRITE (numout,*) " max and min z-vertical level",jpkmax+1,jpup
1205
1206   END SUBROUTINE  fszref
1207
1208
1209   SUBROUTINE zgr_sco
1210      !!----------------------------------------------------------------------
1211      !!                  ***  ROUTINE zgr_sco  ***
1212      !!                     
1213      !! ** Purpose :   define the s-coordinate system
1214      !!
1215      !! ** Method  :   s-coordinate
1216      !!         The depth of model levels is defined as the product of an
1217      !!      analytical function by the local bathymetry, while the vertical
1218      !!      scale factors are defined as the product of the first derivative
1219      !!      of the analytical function by the bathymetry.
1220      !!      (this solution save memory as depth and scale factors are not
1221      !!      3d fields)
1222      !!          - Read bathymetry (in meters) at t-point and compute the
1223      !!         bathymetry at u-, v-, and f-points.
1224      !!            hbatu = mi( hbatt )
1225      !!            hbatv = mj( hbatt )
1226      !!            hbatf = mi( mj( hbatt ) )
1227      !!          - Compute gsigt, gsigw, esigt, esigw from an analytical
1228      !!         function and its derivative given as function.
1229      !!            gsigt(k) = fssig (k    )
1230      !!            gsigw(k) = fssig (k-0.5)
1231      !!            esigt(k) = fsdsig(k    )
1232      !!            esigw(k) = fsdsig(k-0.5)
1233      !!      This routine is given as an example, it must be modified
1234      !!      following the user s desiderata. nevertheless, the output as
1235      !!      well as the way to compute the model levels and scale factors
1236      !!      must be respected in order to insure second order a!!uracy
1237      !!      schemes.
1238      !!
1239      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
1240      !!----------------------------------------------------------------------
1241      !
1242      INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
1243      INTEGER  ::   inum                      ! temporary logical unit
1244      INTEGER  ::   iip1, ijp1, iim1, ijm1, kdep   ! temporary integers
1245      REAL(wp) ::   zcoeft, zcoefw, zrmax, ztaper, maxzenv   ! temporary scalars
1246#if defined key_melange
1247      REAL(wp) ::   rn_hc_bak   ! temporary scalars
1248#endif
1249      REAL(wp) ::   zrfact   ! temporary scalars
1250      REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
1251      !
1252#if defined key_fudge
1253      REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, zri, zrj, zhbat, fenv
1254#else
1255      REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, zri, zrj, zhbat
1256#endif
1257#if defined key_smsh
1258      REAL(wp), POINTER, DIMENSION(:,:,:) :: esigtu3, esigtv3, esigtf3, esigwu3, esigwv3           
1259#else
1260      REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3
1261      REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3           
1262#endif
1263      NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc
1264#if defined key_melange
1265      NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc, nn_sig_lev
1266#endif
1267      !!----------------------------------------------------------------------
1268      !
1269      IF( nn_timing == 1 )  CALL timing_start('zgr_sco')
1270      !
1271      CALL wrk_alloc( jpi, jpj,      ztmpi1, ztmpi2, ztmpj1, ztmpj2         )
1272#if defined key_fudge
1273      CALL wrk_alloc( jpi, jpj,      zenv, zri, zrj, zhbat, fenv     )
1274#else
1275      CALL wrk_alloc( jpi, jpj,      zenv, zri, zrj, zhbat           )
1276#endif
1277#if defined key_smsh
1278      CALL wrk_alloc( jpi, jpj, jpk, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 )
1279#else
1280      CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      )
1281      CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 )
1282#endif     
1283      REWIND( numnam )                       ! Read Namelist namzgr_sco : sigma-stretching parameters
1284      READ  ( numnam, namzgr_sco )
1285
1286      IF(lwp) THEN                           ! control print
1287         WRITE(numout,*)
1288         WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate'
1289         WRITE(numout,*) '~~~~~~~~~~~'
1290         WRITE(numout,*) '   Namelist namzgr_sco'
1291         WRITE(numout,*) '      sigma-stretching coeffs '
1292         WRITE(numout,*) '      maximum depth of s-bottom surface (>0)       rn_sbot_max   = ' ,rn_sbot_max
1293         WRITE(numout,*) '      minimum depth of s-bottom surface (>0)       rn_sbot_min   = ' ,rn_sbot_min
1294         WRITE(numout,*) '      surface control parameter (0<=rn_theta<=20)  rn_theta      = ', rn_theta
1295         WRITE(numout,*) '      bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ', rn_thetb
1296         WRITE(numout,*) '      maximum cut-off r-value allowed              rn_rmax       = ', rn_rmax
1297         WRITE(numout,*) '      Hybrid s-sigma-coordinate                    ln_s_sigma    = ', ln_s_sigma
1298         WRITE(numout,*) '      stretching parameter (song and haidvogel)    rn_bb         = ', rn_bb
1299         WRITE(numout,*) '      Critical depth                               rn_hc         = ', rn_hc
1300      ENDIF
1301     
1302#if defined key_melange
1303      CALL zgr_zps          ! Partial step z-coordinate
1304      ! Scale factors and depth at U-, V-, UW and VW-points
1305      DO jk = 1, nn_sig_lev                        ! initialisation to z-scale factors above ln_s_sigma to remove any zps
1306         e3u (:,:,jk) = e3t_0(jk)
1307         e3v (:,:,jk) = e3t_0(jk)
1308         e3uw(:,:,jk) = e3w_0(jk)
1309         e3vw(:,:,jk) = e3w_0(jk)
1310      END DO
1311#endif
1312
1313      gsigw3  = 0._wp   ;   gsigt3  = 0._wp   ;   gsi3w3  = 0._wp
1314      esigt3  = 0._wp   ;   esigw3  = 0._wp 
1315      esigtu3 = 0._wp   ;   esigtv3 = 0._wp   ;   esigtf3 = 0._wp
1316      esigwu3 = 0._wp   ;   esigwv3 = 0._wp
1317
1318      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate
1319      hifu(:,:) = rn_sbot_min
1320      hifv(:,:) = rn_sbot_min
1321      hiff(:,:) = rn_sbot_min
1322
1323      !                                        ! set maximum ocean depth
1324      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1325
1326      DO jj = 1, jpj
1327         DO ji = 1, jpi
1328           IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1329         END DO
1330      END DO
1331      !                                        ! =============================
1332      !                                        ! Define the envelop bathymetry   (hbatt)
1333      !                                        ! =============================
1334      ! use r-value to create hybrid coordinates
1335      zenv(:,:) = bathy(:,:)
1336#if defined key_melange
1337      DO jj = 1, jpj
1338         DO ji = 1, jpi
1339           zenv(ji,jj) = MIN( bathy(ji,jj), gdepw_0(nn_sig_lev + 1) )
1340         ENDDO
1341      ENDDO
1342#endif
1343#if defined key_fudge
1344            CALL iom_open ( 'fudge.nc', inum ) 
1345            CALL iom_get  ( inum, jpdom_data, 'zenv', fenv )
1346            CALL iom_close( inum )
1347      DO jj = 1, jpj
1348         DO ji = 1, jpi
1349              zenv(ji,jj) = MAX( zenv(ji,jj), fenv(ji,jj) )
1350         ENDDO
1351      ENDDO
1352#endif
1353      !
1354      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
1355      DO jj = 1, jpj
1356         DO ji = 1, jpi
1357           IF( bathy(ji,jj) == 0._wp ) THEN
1358             iip1 = MIN( ji+1, jpi )
1359             ijp1 = MIN( jj+1, jpj )
1360             iim1 = MAX( ji-1, 1 )
1361             ijm1 = MAX( jj-1, 1 )
1362             IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) +              &
1363        &         bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
1364               zenv(ji,jj) = rn_sbot_min
1365             ENDIF
1366           ENDIF
1367         END DO
1368      END DO
1369      ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1370      CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1371      !
1372      ! smooth the bathymetry (if required)
1373      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
1374      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
1375      !
1376      jl = 0
1377      zrmax = 1._wp
1378      !     
1379      ! set scaling factor used in reducing vertical gradients
1380      zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax ) 
1381      !
1382      ! initialise temporary evelope depth arrays
1383      ztmpi1(:,:) = zenv(:,:)
1384      ztmpi2(:,:) = zenv(:,:)
1385      ztmpj1(:,:) = zenv(:,:)
1386      ztmpj2(:,:) = zenv(:,:)
1387      !
1388      ! initialise temporary r-value arrays
1389      zri(:,:) = 1._wp
1390      zrj(:,:) = 1._wp
1391      !                                                            ! ================ !
1392      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) !  Iterative loop  !
1393         !                                                         ! ================ !
1394         jl = jl + 1
1395         zrmax = 0._wp
1396         ! we set zrmax from previous r-values (zri and zrj) first
1397         ! if set after current r-value calculation (as previously)
1398         ! we could exit DO WHILE prematurely before checking r-value
1399         ! of current zenv
1400         DO jj = 1, nlcj
1401            DO ji = 1, nlci
1402               zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
1403            END DO
1404         END DO
1405         zri(:,:) = 0._wp
1406         zrj(:,:) = 0._wp
1407         DO jj = 1, nlcj
1408            DO ji = 1, nlci
1409               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
1410               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
1411               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
1412                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
1413               END IF
1414               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
1415                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
1416               END IF
1417               IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact
1418               IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact 
1419               IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact
1420               IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact
1421            END DO
1422         END DO
1423         IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
1424         !
1425         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax
1426         !
1427         DO jj = 1, nlcj
1428            DO ji = 1, nlci
1429               zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
1430            END DO
1431         END DO
1432         ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1433         CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1434         !                                                  ! ================ !
1435      END DO                                                !     End loop     !
1436      !                                                     ! ================ !
1437      DO jj = 1, jpj
1438         DO ji = 1, jpi
1439            zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
1440         END DO
1441      END DO
1442      !
1443      ! Envelope bathymetry saved in hbatt
1444      hbatt(:,:) = zenv(:,:) 
1445
1446      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
1447         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1448         DO jj = 1, jpj
1449            DO ji = 1, jpi
1450               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
1451               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
1452            END DO
1453         END DO
1454      ENDIF
1455      !
1456      IF(lwp) THEN                             ! Control print
1457         WRITE(numout,*)
1458         WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
1459         WRITE(numout,*)
1460         CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )
1461         IF( nprint == 1 )   THEN       
1462            WRITE(numout,*) ' bathy  MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
1463            WRITE(numout,*) ' hbatt  MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
1464         ENDIF
1465      ENDIF
1466
1467      !                                        ! ==============================
1468      !                                        !   hbatu, hbatv, hbatf fields
1469      !                                        ! ==============================
1470      IF(lwp) THEN
1471         WRITE(numout,*)
1472         WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
1473      ENDIF
1474      hbatu(:,:) = rn_sbot_min
1475      hbatv(:,:) = rn_sbot_min
1476      hbatf(:,:) = rn_sbot_min
1477      DO jj = 1, jpjm1
1478        DO ji = 1, jpim1   ! NO vector opt.
1479           hbatu(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
1480           hbatv(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
1481           hbatf(ji,jj) = 0.25_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
1482              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
1483        END DO
1484      END DO
1485      !
1486      ! Apply lateral boundary condition
1487!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
1488      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk( hbatu, 'U', 1._wp )
1489      DO jj = 1, jpj
1490         DO ji = 1, jpi
1491            IF( hbatu(ji,jj) == 0._wp ) THEN
1492               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
1493               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
1494            ENDIF
1495         END DO
1496      END DO
1497      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk( hbatv, 'V', 1._wp )
1498      DO jj = 1, jpj
1499         DO ji = 1, jpi
1500            IF( hbatv(ji,jj) == 0._wp ) THEN
1501               IF( zhbat(ji,jj) == 0._wp )   hbatv(ji,jj) = rn_sbot_min
1502               IF( zhbat(ji,jj) /= 0._wp )   hbatv(ji,jj) = zhbat(ji,jj)
1503            ENDIF
1504         END DO
1505      END DO
1506      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk( hbatf, 'F', 1._wp )
1507      DO jj = 1, jpj
1508         DO ji = 1, jpi
1509            IF( hbatf(ji,jj) == 0._wp ) THEN
1510               IF( zhbat(ji,jj) == 0._wp )   hbatf(ji,jj) = rn_sbot_min
1511               IF( zhbat(ji,jj) /= 0._wp )   hbatf(ji,jj) = zhbat(ji,jj)
1512            ENDIF
1513         END DO
1514      END DO
1515
1516!!bug:  key_helsinki a verifer
1517      hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1518      hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1519      hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1520      hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1521
1522      IF( nprint == 1 .AND. lwp )   THEN
1523         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
1524            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
1525         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
1526            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
1527         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
1528            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1529         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
1530            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1531      ENDIF
1532!! helsinki
1533
1534      !                                            ! =======================
1535      !                                            !   s-ccordinate fields     (gdep., e3.)
1536      !                                            ! =======================
1537      !
1538      ! non-dimensional "sigma" for model level depth at w- and t-levels
1539
1540      IF( ln_s_sigma ) THEN        ! Song and Haidvogel style stretched sigma for depths
1541         !                         ! below rn_hc, with uniform sigma in shallower waters
1542         DO ji = 1, jpi
1543            DO jj = 1, jpj
1544
1545               IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
1546                  DO jk = 1, jpk
1547#if defined key_melange
1548                     gsigw3(ji,jj,jk) = gdepw_0(jk)/gdepw_0(nn_sig_lev + 1)
1549                     gsigt3(ji,jj,jk) = gdept_0(jk)/gdepw_0(nn_sig_lev + 1)
1550#else
1551                     gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
1552                     gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
1553#endif
1554                  END DO
1555               ELSE ! shallow water, uniform sigma
1556                  DO jk = 1, jpk
1557#if defined key_melange
1558                     gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(nn_sig_lev,wp)
1559                     gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(nn_sig_lev,wp)
1560#else
1561                     gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
1562                     gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
1563#endif
1564                  END DO
1565               ENDIF
1566               IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'gsigw3 1 jpk    ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk)
1567               !
1568               DO jk = 1, jpkm1
1569                  esigt3(ji,jj,jk  ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk)
1570                  esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk)
1571               END DO
1572               esigw3(ji,jj,1  ) = 2._wp * ( gsigt3(ji,jj,1  ) - gsigw3(ji,jj,1  ) )
1573               esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) )
1574               !
1575               ! Coefficients for vertical depth as the sum of e3w scale factors
1576               gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1)
1577               DO jk = 2, jpk
1578                  gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk)
1579               END DO
1580               !
1581#if defined key_melange
1582               DO jk = 1, nn_sig_lev+1
1583!              DO jk = 1, jpk
1584               IF( bathy(ji,jj) < gdepw_0(nn_sig_lev + 1) ) THEN ! should this be bathy or hbatt?
1585#else
1586               DO jk = 1, jpk
1587#endif
1588#if defined key_melange
1589                  zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(nn_sig_lev,wp)
1590                  zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(nn_sig_lev,wp)
1591!                 zcoeft = ( REAL(MIN(jk,nn_sig_lev),wp) - 0.5_wp ) / REAL(nn_sig_lev-1,wp)
1592!                 zcoefw = ( REAL(MIN(jk,nn_sig_lev),wp) - 1.0_wp ) / REAL(nn_sig_lev-1,wp)
1593#else
1594                  zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1595                  zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1596#endif
1597#if defined key_melange
1598                  rn_hc_bak = rn_hc
1599                  rn_hc = MIN( MAX ( &
1600                &            (hbatt(ji,jj)-gdepw_0(nn_sig_lev + 1)) / (1._wp - (gdepw_0(nn_sig_lev + 1)/rn_hc)) &
1601                &                   ,0._wp) ,rn_hc)
1602#endif
1603                  gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft )
1604                  gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw )
1605                  gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
1606#if defined key_melange
1607                  rn_hc = rn_hc_bak
1608#endif
1609               IF( gdepw(ji,jj,jk) < 0._wp ) THEN
1610                  WRITE(*,*) 'zgr_sco :   gdepw  at point (i,j,k)= ', ji, jj, jk, (gsigw3(ji,jj,jk)*10000._wp-zcoefw*10000._wp)
1611               ENDIF
1612#if defined key_melange
1613               ENDIF
1614#endif
1615               END DO
1616               !
1617            END DO   ! for all jj's
1618         END DO    ! for all ji's
1619
1620         DO ji = 1, jpim1
1621            DO jj = 1, jpjm1
1622#if defined key_melange
1623               IF( bathy(ji,jj) < gdepw_0(nn_sig_lev + 1) ) THEN ! should this be bathy or hbatt?
1624               DO jk = 1, nn_sig_lev+1      ! scale factors should be the same in both zps and sco when H > Hcrit??
1625!              DO jk = 1, jpk
1626#else
1627               DO jk = 1, jpk
1628#endif
1629                  esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) )   &
1630                     &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1631                  esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) )   &
1632                     &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1633                  esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk)     &
1634                     &                + hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) )   &
1635                     &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
1636                  esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) )   &
1637                     &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1638                  esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) )   &
1639                     &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1640                  !
1641#if defined key_melange
1642                  rn_hc_bak = rn_hc
1643                  rn_hc = MIN( MAX( &
1644                &            (hbatt(ji,jj)-gdepw_0(nn_sig_lev + 1)) / (1._wp - (gdepw_0(nn_sig_lev + 1)/rn_hc)) &
1645                &                   ,0._wp) ,rn_hc)
1646!                 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) )
1647!                 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) )
1648                  e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) )
1649                  e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) )
1650#else
1651                  e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1652                  e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1653#endif
1654#if defined key_melange
1655                  rn_hc = MIN( MAX( &
1656                &            (hbatu(ji,jj)-gdepw_0(nn_sig_lev + 1)) / (1._wp - (gdepw_0(nn_sig_lev + 1)/rn_hc_bak)) &
1657                &                   ,0._wp) ,rn_hc_bak)
1658!                 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) )
1659!                 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) )
1660                  e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) )
1661                  e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) )
1662#else
1663                  e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1664                  e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1665#endif
1666#if defined key_melange
1667                  rn_hc = MIN( MAX( &
1668                &            (hbatv(ji,jj)-gdepw_0(nn_sig_lev + 1)) / (1._wp - (gdepw_0(nn_sig_lev + 1)/rn_hc_bak)) &
1669                &                   ,0._wp) ,rn_hc_bak)
1670!                 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) )
1671!                 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) )
1672                  e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) )
1673                  e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) )
1674#else
1675                  e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1676                  e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1677#endif
1678#if defined key_melange
1679                  rn_hc = MIN( MAX( &
1680                &            (hbatf(ji,jj)-gdepw_0(nn_sig_lev + 1)) / (1._wp - (gdepw_0(nn_sig_lev + 1)/rn_hc_bak)) &
1681                &                   ,0._wp), rn_hc_bak)
1682!                 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) )
1683                  e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) )
1684#else
1685                  e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp))
1686#endif
1687                  !
1688#if defined key_melange
1689                  rn_hc = rn_hc_bak
1690#endif
1691               END DO
1692#if defined key_melange
1693               ENDIF
1694#endif
1695            END DO
1696         END DO
1697
1698         CALL lbc_lnk( e3t , 'T', 1._wp )
1699         CALL lbc_lnk( e3u , 'U', 1._wp )
1700         CALL lbc_lnk( e3v , 'V', 1._wp )
1701         CALL lbc_lnk( e3f , 'F', 1._wp )
1702         CALL lbc_lnk( e3w , 'W', 1._wp )
1703         CALL lbc_lnk( e3uw, 'U', 1._wp )
1704         CALL lbc_lnk( e3vw, 'V', 1._wp )
1705
1706         !
1707      ELSE   ! not ln_s_sigma
1708         !
1709         DO jk = 1, jpk
1710           gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
1711           gsigt(jk) = -fssig( REAL(jk,wp)        )
1712         END DO
1713         IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk)
1714         !
1715         ! Coefficients for vertical scale factors at w-, t- levels
1716!!gm bug :  define it from analytical function, not like juste bellow....
1717!!gm        or betteroffer the 2 possibilities....
1718         DO jk = 1, jpkm1
1719            esigt(jk  ) = gsigw(jk+1) - gsigw(jk)
1720            esigw(jk+1) = gsigt(jk+1) - gsigt(jk)
1721         END DO
1722         esigw( 1 ) = 2._wp * ( gsigt(1  ) - gsigw(1  ) ) 
1723         esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) )
1724
1725!!gm  original form
1726!!org DO jk = 1, jpk
1727!!org    esigt(jk)=fsdsig( FLOAT(jk)     )
1728!!org    esigw(jk)=fsdsig( FLOAT(jk)-0.5 )
1729!!org END DO
1730!!gm
1731         !
1732         ! Coefficients for vertical depth as the sum of e3w scale factors
1733         gsi3w(1) = 0.5_wp * esigw(1)
1734         DO jk = 2, jpk
1735            gsi3w(jk) = gsi3w(jk-1) + esigw(jk)
1736         END DO
1737!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
1738         DO jk = 1, jpk
1739            zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1740            zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1741            gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigt(jk) + hift(:,:)*zcoeft )
1742            gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigw(jk) + hift(:,:)*zcoefw )
1743            gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsi3w(jk) + hift(:,:)*zcoeft )
1744         END DO
1745!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
1746         DO jj = 1, jpj
1747            DO ji = 1, jpi
1748               DO jk = 1, jpk
1749                 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1750                 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1751                 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1752                 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
1753                 !
1754                 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1755                 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1756                 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1757               END DO
1758            END DO
1759         END DO
1760         !
1761      ENDIF ! ln_s_sigma
1762
1763      where (e3t   (:,:,:).eq.0._wp)  e3t(:,:,:) = 1._wp
1764      where (e3u   (:,:,:).eq.0._wp)  e3u(:,:,:) = 1._wp
1765      where (e3v   (:,:,:).eq.0._wp)  e3v(:,:,:) = 1._wp
1766      where (e3f   (:,:,:).eq.0._wp)  e3f(:,:,:) = 1._wp
1767      where (e3w   (:,:,:).eq.0._wp)  e3w(:,:,:) = 1._wp
1768      where (e3uw  (:,:,:).eq.0._wp)  e3uw(:,:,:) = 1._wp
1769      where (e3vw  (:,:,:).eq.0._wp)  e3vw(:,:,:) = 1._wp
1770
1771
1772      fsdept(:,:,:) = gdept (:,:,:)
1773      fsdepw(:,:,:) = gdepw (:,:,:)
1774      fsde3w(:,:,:) = gdep3w(:,:,:)
1775      fse3t (:,:,:) = e3t   (:,:,:)
1776      fse3u (:,:,:) = e3u   (:,:,:)
1777      fse3v (:,:,:) = e3v   (:,:,:)
1778      fse3f (:,:,:) = e3f   (:,:,:)
1779      fse3w (:,:,:) = e3w   (:,:,:)
1780      fse3uw(:,:,:) = e3uw  (:,:,:)
1781      fse3vw(:,:,:) = e3vw  (:,:,:)
1782!!
1783      ! HYBRID :
1784      DO jj = 1, jpj
1785         DO ji = 1, jpi
1786#if defined key_melange
1787               IF( bathy(ji,jj) < gdepw_0(nn_sig_lev + 1) ) THEN ! should this be hbatt or bathy
1788               DO jk = 1, nn_sig_lev 
1789!              DO jk = 1, jpkm1
1790#else
1791               DO jk = 1, jpkm1
1792#endif
1793                   IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX(  2, jk )
1794               END DO
1795#if defined key_melange
1796               ENDIF
1797#endif
1798               IF( scobot(ji,jj) == 0._wp            )   mbathy(ji,jj) = 0
1799         END DO
1800      END DO
1801      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
1802         &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
1803
1804      !                                               ! =============
1805      IF(lwp) THEN                                    ! Control print
1806         !                                            ! =============
1807         WRITE(numout,*) 
1808         WRITE(numout,*) ' domzgr: vertical coefficients for model level'
1809         WRITE(numout, "(9x,'  level    gsigt      gsigw      esigt      esigw      gsi3w')" )
1810         WRITE(numout, "(10x,i4,5f11.4)" ) ( jk, gsigt(jk), gsigw(jk), esigt(jk), esigw(jk), gsi3w(jk), jk=1,jpk )
1811      ENDIF
1812      IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain
1813         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)   ), ' MAX ', MAXVAL( mbathy(:,:) )
1814         WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ),   &
1815            &                          ' w ', MINVAL( fsdepw(:,:,:) ), '3w '  , MINVAL( fsde3w(:,:,:) )
1816         WRITE(numout,*) ' MIN val e3    t ', MINVAL( fse3t (:,:,:) ), ' f '  , MINVAL( fse3f (:,:,:) ),   &
1817            &                          ' u ', MINVAL( fse3u (:,:,:) ), ' u '  , MINVAL( fse3v (:,:,:) ),   &
1818            &                          ' uw', MINVAL( fse3uw(:,:,:) ), ' vw'  , MINVAL( fse3vw(:,:,:) ),   &
1819            &                          ' w ', MINVAL( fse3w (:,:,:) )
1820
1821         WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ),   &
1822            &                          ' w ', MAXVAL( fsdepw(:,:,:) ), '3w '  , MAXVAL( fsde3w(:,:,:) )
1823         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( fse3t (:,:,:) ), ' f '  , MAXVAL( fse3f (:,:,:) ),   &
1824            &                          ' u ', MAXVAL( fse3u (:,:,:) ), ' u '  , MAXVAL( fse3v (:,:,:) ),   &
1825            &                          ' uw', MAXVAL( fse3uw(:,:,:) ), ' vw'  , MAXVAL( fse3vw(:,:,:) ),   &
1826            &                          ' w ', MAXVAL( fse3w (:,:,:) )
1827      ENDIF
1828      !
1829      IF(lwp) THEN                                  ! selected vertical profiles
1830         WRITE(numout,*)
1831         WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
1832         WRITE(numout,*) ' ~~~~~~  --------------------'
1833         WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
1834         WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk),     &
1835            &                                 fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk )
1836         DO jj = mj0(20), mj1(20)
1837            DO ji = mi0(20), mi1(20)
1838               WRITE(numout,*)
1839               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
1840               WRITE(numout,*) ' ~~~~~~  --------------------'
1841               WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
1842               WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk),     &
1843                  &                                 fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
1844            END DO
1845         END DO
1846         DO jj = mj0(74), mj1(74)
1847            DO ji = mi0(100), mi1(100)
1848               WRITE(numout,*)
1849               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
1850               WRITE(numout,*) ' ~~~~~~  --------------------'
1851               WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
1852               WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk),     &
1853                  &                                 fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
1854            END DO
1855         END DO
1856      ENDIF
1857
1858!!gm bug?  no more necessary?  if ! defined key_helsinki
1859      DO jk = 1, jpk
1860         DO jj = 1, jpj
1861            DO ji = 1, jpi
1862               IF( fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
1863                  WRITE(*,*) 'zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk, fse3w(ji,jj,jk), fse3t(ji,jj,jk)
1864!                 WRITE(ctmp1,*) 'zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
1865!                 CALL ctl_stop( ctmp1 )
1866               ENDIF
1867               IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
1868                  WRITE(*,*) 'zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk, fsdepw(ji,jj,jk), fsdept(ji,jj,jk)
1869!                 WRITE(ctmp1,*) 'zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
1870!                 CALL ctl_stop( ctmp1 )
1871               ENDIF
1872            END DO
1873         END DO
1874      END DO
1875!!gm bug    #endif
1876      !
1877
1878      CALL wrk_dealloc( jpi, jpj,      zenv, ztmpi1, ztmpi2, ztmpj1, ztmpj2, zri, zrj, zhbat                           )
1879
1880#if defined key_smsh
1881      CALL wrk_dealloc( jpi, jpj, jpk, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 )
1882#else
1883      CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      )
1884      CALL wrk_dealloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 )
1885#endif
1886     !
1887      IF( nn_timing == 1 )  CALL timing_stop('zgr_sco')
1888      !
1889   END SUBROUTINE zgr_sco
1890      SUBROUTINE zgr_hyb
1891      !!----------------------------------------------------------------------
1892      !!                  ***  ROUTINE zgr_sco  ***
1893      !!    Combination of zgr_sco in upper layers ( shelf ) and zgr_zps in abyss      !!                     
1894      !! ** Purpose :   define the s-z coordinate system
1895      !!
1896      !! ** Method  :   s-coordinate in upper layers and z-coordinates below
1897      !!         The depth of model levels is defined as the product of an
1898      !!      analytical function by the local bathymetry, while the vertical
1899      !!      scale factors are defined as the product of the first derivative
1900      !!      of the analytical function by the bathymetry.
1901      !!      (this solution save memory as depth and scale factors are not
1902      !!      3d fields)
1903      !!          - Read bathymetry (in meters) at t-point and compute the
1904      !!         bathymetry at u-, v-, and f-points.
1905      !!            hbatu = mi( hbatt )
1906      !!            hbatv = mj( hbatt )
1907      !!            hbatf = mi( mj( hbatt ) )
1908      !!          - Compute gsigt, gsigw, esigt, esigw from an analytical
1909      !!         function
1910      !!            gsigt(k) = fssig (k    )
1911      !!            gsigw(k) = fssig (k-0.5)
1912      !!      This routine is given as an example, it must be modified
1913      !!      following the user s desiderata. nevertheless, the output as
1914      !!      well as the way to compute the model levels and scale factors
1915      !!      must be respected in order to insure second order a!!uracy
1916      !!      schemes.
1917      !!
1918
1919  !!======================================================================
1920      INTEGER  ::   ji, jj, jk, jl, ik           ! dummy loop argument
1921      INTEGER  ::   iip1, ijp1, iim1, ijm1       ! temporary integers
1922      INTEGER  ::   jpksigm                      ! temporary integer for maxnumber of s-levels
1923      REAL(wp) ::   zcoeft, zcoefw, zrmax, ztaper,zrmin,e3t_t,e3w_t  ! temporary scalars
1924      REAL(wp) ::   ze3tp , ze3wp           ! Last ocean level thickness at T- and W-points
1925      REAL(wp) ::   zdepwp, zdepth          ! Ajusted ocean depth to avoid too small e3t
1926      REAL(wp) ::   zmax, zmin ,zsigma      ! Maximum and minimum depth and depth of sigma layer
1927      REAL(wp) ::   zacr , zkth ,za1         ! parameters for z- layer (as ppacr , ppzkth )
1928
1929
1930      !
1931      REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat ,zrpt
1932
1933      NAMELIST/namzgr_hyb/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, &
1934                 ln_s_sigma, rn_bb, rn_hc,rn_zsigma , nn_sig_lev , rn_kth , rn_acr
1935
1936
1937      !!----------------------------------------------------------------------
1938      !
1939      IF( nn_timing == 1 )  CALL timing_start('zgr_hyb')
1940      !
1941      CALL wrk_alloc( jpi, jpj,      zenv, ztmp, zmsk, zri, zrj, zhbat                           )
1942
1943!      CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3                                              )
1944      !
1945      REWIND( numnam )                       ! Read Namelist namzgr_sco : sigma-stretching parameters
1946      READ  ( numnam, namzgr_hyb )
1947      IF(lwp) THEN                           ! control print
1948         WRITE(numout,*)
1949         WRITE(numout,*) 'dom:zgr_hyb : s-coordinate or hybrid z-s-coordinate'
1950         WRITE(numout,*) '~~~~~~~~~~~'
1951         WRITE(numout,*) '   Namelist namzgr_hyb'
1952         WRITE(numout,*) '      sigma-stretching coeffs '
1953         WRITE(numout,*) '      maximum depth of s-bottom surface (>0)       rn_sbot_max   = ' ,rn_sbot_max
1954         WRITE(numout,*) '      minimum depth of s-bottom surface (>0)       rn_sbot_min   = ' ,rn_sbot_min
1955         WRITE(numout,*) '      surface control parameter (0<=rn_theta<=20)  rn_theta      = ', rn_theta
1956         WRITE(numout,*) '      bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ', rn_thetb
1957         WRITE(numout,*) '      maximum cut-off r-value allowed              rn_rmax       = ', rn_rmax
1958         WRITE(numout,*) '      Hybrid s-sigma-coordinate                    ln_s_sigma    = ', ln_s_sigma
1959         WRITE(numout,*) '      stretching parameter (song and haidvogel)    rn_bb         = ', rn_bb
1960         WRITE(numout,*) '      Critical depth                               rn_hc         = ', rn_hc
1961         WRITE(numout,*) '      Sigma  depth                                rn_zsigma     = ', rn_zsigma
1962         WRITE(numout,*) '      The same as pp_arc                           rn_arc        = ', rn_acr
1963         WRITE(numout,*) '      Number of sigma levels                       rn_arc        = ', nn_sig_lev
1964         WRITE(numout,*) '      Number of levels for stretching z-coord      rn_kth        = ', rn_kth
1965
1966
1967      ENDIF
1968      zsigma   =    rn_zsigma
1969      jpksigm =  nn_sig_lev
1970      zmax    =  rn_sbot_max
1971      zacr    =  rn_acr
1972      zkth    =  rn_kth
1973      e3t(:,:,:) = 1._wp     
1974      e3w(:,:,:) = 1._wp
1975      e3u(:,:,:) = 1._wp
1976      e3v(:,:,:) = 1._wp
1977      e3f(:,:,:) = 1._wp
1978      e3uw(:,:,:)= 1._wp
1979      e3vw(:,:,:)= 1._wp
1980
1981
1982
1983
1984
1985      DO jj = 1, jpj
1986         DO ji= 1, jpi
1987            IF( bathy(ji,jj) <= 0._wp ) THEN   ;   bathy(ji,jj) = 0.e0_wp
1988            ELSE                               ;   bathy(ji,jj) = MIN( rn_sbot_max,  MAX( bathy(ji,jj),rn_sbot_min )  )
1989            ENDIF
1990         END DO
1991      END DO
1992
1993! create bathymetry for enveloping
1994      DO jj = 1, jpj
1995         DO ji = 1, jpi
1996            zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min )
1997            zenv(ji,jj) = MIN (zenv(ji,jj),  zsigma )
1998            hbatt(ji,jj) = zenv(ji,jj)
1999         END DO
2000      END DO
2001      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
2002      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
2003      jl = 0
2004      zrmax = 1._wp
2005      !                                                     ! ================ !
2006      DO WHILE( jl <= 10000 .AND. zrmax > rn_rmax )         !  Iterative loop  !
2007         !                                                  ! ================ !
2008         jl = jl + 1
2009         zrmax = 0._wp
2010         zmsk(:,:) = 0._wp
2011         DO jj = 1, nlcj
2012            DO ji = 1, nlci
2013               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
2014               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
2015               zri(ji,jj) = ABS( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
2016               zrj(ji,jj) = ABS( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
2017               zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) )
2018               IF( zri(ji,jj) > rn_rmax )   zmsk(ji  ,jj  ) = 1._wp
2019               IF( zri(ji,jj) > rn_rmax )   zmsk(iip1,jj  ) = 1._wp
2020               IF( zrj(ji,jj) > rn_rmax )   zmsk(ji  ,jj  ) = 1._wp
2021               IF( zrj(ji,jj) > rn_rmax )   zmsk(ji  ,ijp1) = 1._wp
2022            END DO
2023         END DO
2024         IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
2025         ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)
2026         ztmp(:,:) = zmsk(:,:)   ;   CALL lbc_lnk( zmsk, 'T', 1._wp )
2027         DO jj = 1, nlcj
2028            DO ji = 1, nlci
2029                zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )
2030            END DO
2031         END DO
2032         IF(lwp)WRITE(numout,*) 'zgr_hyb :   iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) )
2033         !
2034         DO jj = 1, nlcj
2035            DO ji = 1, nlci
2036               iip1 = MIN( ji+1, nlci )     ! last  line (ji=nlci)
2037               ijp1 = MIN( jj+1, nlcj )     ! last  raw  (jj=nlcj)
2038               iim1 = MAX( ji-1,  1  )      ! first line (ji=nlci)
2039               ijm1 = MAX( jj-1,  1  )      ! first raw  (jj=nlcj)
2040               IF( zmsk(ji,jj) == 1._wp ) THEN
2041                  ztmp(ji,jj) =   (                                                                                   &
2042             &      zenv(iim1,ijp1)*zmsk(iim1,ijp1) + zenv(ji,ijp1)*zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1)   &
2043             &    + zenv(iim1,jj  )*zmsk(iim1,jj  ) + zenv(ji,jj  )*    2._wp     + zenv(iip1,jj  )*zmsk(iip1,jj  )   &
2044             &    + zenv(iim1,ijm1)*zmsk(iim1,ijm1) + zenv(ji,ijm1)*zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1)   &
2045             &                    ) / (                                                                               &
2046             &                      zmsk(iim1,ijp1) +               zmsk(ji,ijp1) +                 zmsk(iip1,ijp1)   &
2047             &    +                 zmsk(iim1,jj  ) +                   2._wp     +                 zmsk(iip1,jj  )   &
2048             &    +                 zmsk(iim1,ijm1) +               zmsk(ji,ijm1) +                 zmsk(iip1,ijm1)   &
2049             &                        )
2050               ENDIF
2051            END DO
2052         END DO
2053         !
2054         DO jj = 1, nlcj
2055            DO ji = 1, nlci
2056               IF( zmsk(ji,jj) == 1._wp )   zenv(ji,jj) = MAX( ztmp(ji,jj), hbatt(ji,jj) )
2057            END DO
2058         END DO
2059         !
2060         ! Apply lateral boundary condition   CAUTION: kept the value when the lbc field is zero
2061         ztmp(:,:) = zenv(:,:)   ;   CALL lbc_lnk( zenv, 'T', 1._wp )
2062         DO jj = 1, nlcj
2063            DO ji = 1, nlci
2064               IF( zenv(ji,jj) == 0._wp )   zenv(ji,jj) = ztmp(ji,jj)
2065            END DO
2066         END DO
2067         !                                                  ! ================ !
2068      END DO                                                !     End loop     !
2069      !                                                     ! ================ !
2070      !
2071      !                                        ! envelop bathymetry saved in hbatt
2072      hbatt(:,:) = zenv(:,:) 
2073   IF( lk_mpp )   CALL mpp_max( nstop )
2074   IF (lwp) write(numout,*)"after envelope", nstop
2075
2076! define new reference levels
2077! for s- levels, allows stretching at surface and bottom layer
2078IF( jpksigm > 1 )THEN
2079    IF(ln_s_sigma)THEN
2080           DO jk = 1, jpksigm
2081            gsigw(jk) = -fssig3( REAL(jk,wp)-0.5_wp     , rn_bb,jpksigm )
2082            gsigt(jk) = -fssig3( REAL(jk,wp)            , rn_bb,jpksigm )
2083           END DO
2084    ELSE
2085          DO jk = 1, jpksigm
2086           gsigw(jk) = -fssig2( REAL(jk,wp)-0.5_wp ,jpksigm )
2087           gsigt(jk) = -fssig2( REAL(jk,wp)        ,jpksigm )
2088
2089          END DO
2090   ENDIF 
2091    gsigw(1)=0._wp   ! set gsigw exactly to zero     
2092    IF( lk_mpp )   CALL mpp_max( nstop )
2093    IF (lwp) THEN
2094    write(numout,*)"after fssig", nstop,"gsigw,gsigt="
2095    do jk=1,jpksigm
2096    write(numout,*)jk,gsigw(jk),gsigt(jk)
2097    enddo
2098    ENDIF   
2099     
2100       DO jk=1,jpksigm
2101          DO jj=1,jpj
2102            DO ji=1,jpi
2103                 zrmin= min( hbatt(ji,jj), zsigma )
2104                   IF(hbatt(ji,jj).lt.rn_hc)THEN
2105                   zcoefw=REAL(jk-1,wp)      / REAL(jpksigm-1,wp)
2106                   zcoeft=(REAL(jk-1,wp)+0.5)/ REAL(jpksigm-1,wp)
2107                   ELSE
2108                   zcoefw=gsigw(jk)
2109                   zcoeft=gsigt(jk)
2110
2111                   ENDIF
2112
2113                   gdept(ji,jj,jk)=scosrf(ji,jj)+(zrmin-rn_hc)*zcoeft &
2114                       +rn_hc* (REAL(jk,wp)- 0.5_wp)  / REAL(jpksigm-1,wp)
2115                   gdepw(ji,jj,jk)=scosrf(ji,jj)+(zrmin-rn_hc)*zcoefw &
2116                       +rn_hc*( REAL(jk,wp)- 1.0_wp ) / REAL(jpksigm-1,wp)
2117
2118           ENDDO
2119          ENDDO
2120        ENDDO
2121        gdepw(:,:,1)= scosrf(:,:)
2122 ! redefine gdept_0, gdepw_0 which will be used in diawri.F90
2123    DO jk=1,jpksigm
2124                   IF(zsigma.lt.rn_hc)THEN
2125                   zcoefw=REAL(jk-1,wp)      / REAL(jpksigm-1,wp)
2126                   zcoeft=(REAL(jk-1,wp)+0.5)/ REAL(jpksigm-1,wp)
2127                   ELSE
2128                   zcoefw=gsigw(jk)
2129                   zcoeft=gsigt(jk)
2130                   ENDIF
2131
2132         gdept_0(jk)=  zcoeft * (zsigma-rn_hc)+rn_hc* (REAL(jk,wp)- 0.5_wp )/REAL(jpksigm-1,wp)
2133         gdepw_0(jk)=  zcoefw * (zsigma-rn_hc)+rn_hc* (REAL(jk,wp)- 1.0_wp )/REAL(jpksigm-1,wp)
2134    ENDDO
2135
2136    DO jk=1,jpksigm-1
2137         e3t_0(jk)  = gdepw_0(jk+1)-gdepw_0(jk)
2138         e3w_0(jk+1)= gdept_0(jk+1)-gdept_0(jk)
2139    ENDDO
2140         e3w_0(1) = 2._wp * ( gdept_0(1  ) - gdepw_0(1  ) ) 
2141
2142
2143 ! now for lower z-levels :
2144         zmin = e3t_0 (jpksigm -1 ) ! min layer width in z- zone is the same as lowest in s- layer
2145 ELSE
2146  zsigma = 0._wp
2147  hbatt(:,:) = zsigma
2148  zmin = 5._wp         
2149 ENDIF
2150 IF(lwp) write(numout,*) ": last vertical level of sigma-coordinates",zmin
2151        CALL fszref (  zkth, zmin, zacr, zmax, jpksigm , zsigma  )
2152
2153
2154    DO jk=jpksigm,jpkm1
2155         
2156         e3t_0(jk)  = gdepw_0(jk+1)-gdepw_0(jk)
2157         e3w_0(jk+1)= gdept_0(jk+1)-gdept_0(jk)
2158    ENDDO
2159         e3w_0(1) = 2._wp * ( gdept_0(1  ) - gdepw_0(1  ) )
2160         e3t_0(jpk) = 2._wp * ( gdept_0(jpk) - gdepw_0(jpk) )
2161         
2162   IF( lk_mpp )   CALL mpp_max( nstop )
2163   IF (lwp) write(numout,*)"e3t0" ,nstop
2164
2165      IF(lwp) THEN                        ! control print
2166         WRITE(numout,*)
2167         WRITE(numout,*) '  zhyb          Reference z-coordinate depth and scale factors:'
2168         WRITE(numout, "(9x,' level   gdept    gdepw     e3t      e3w  ')" )
2169         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_0(jk), gdepw_0(jk), e3t_0(jk), e3w_0(jk), jk = 1, jpk )
2170         open(333,file='zmesh.dat')
2171         WRITE(333,*)'initial zsigma  =', zsigma, jpk
2172         WRITE(333,*)'initial jpksigm =', jpksigm
2173         WRITE(333,*)'rn_bb =', rn_bb,'rn_theta=',rn_theta
2174         do jk=1,jpk
2175         WRITE(333,'(i4,1x,4(1x,e13.6))') jk, gdept_0(jk), gdepw_0(jk), e3t_0(jk),e3w_0(jk) 
2176         enddo
2177         close(333)
2178      ENDIF
2179
2180        DO jk=jpksigm,jpk
2181          DO jj=1,jpj
2182           DO ji=1,jpi
2183           IF(jpksigm>1)THEN
2184                   gdept(ji,jj,jk)=scosrf(ji,jj) + gdept_0(jk) *hbatt(ji,jj)/zsigma ! differ from gdept0+scorf only at land
2185                   gdepw(ji,jj,jk)=scosrf(ji,jj) + gdepw_0(jk) *hbatt(ji,jj)/zsigma ! as hbatt=zsigma over deep part of basin
2186            ELSE
2187                   gdept(ji,jj,jk)=scosrf(ji,jj) + gdept_0(jk) 
2188                   gdepw(ji,jj,jk)=scosrf(ji,jj) + gdepw_0(jk) 
2189           
2190            ENDIF
2191           ENDDO
2192          ENDDO
2193        ENDDO
2194
2195! define e3t, e3w for general levels
2196        DO jk=1,jpkm1
2197            e3t(:,:,jk)  = gdepw(:,:,jk+1)-gdepw(:,:,jk)
2198            e3w(:,:,jk+1)= gdept(:,:,jk+1)-gdept(:,:,jk)
2199        ENDDO
2200        e3w(:,:,1)   = 2._wp * ( gdept(:,:,1  ) - gdepw(:,:,1  ) )
2201        e3t(:,:,jpk) = 2._wp * ( gdept(:,:,jpk) - gdepw(:,:,jpk) )
2202
2203! and surface :
2204!      ! HYBRID mbathy :
2205     
2206      IF(lwp) THEN                        ! control print
2207         WRITE(numout,*)
2208         WRITE(numout,*) '  zhyb          centre of basin s-z-coordinate depth and scale factors:'
2209         WRITE(numout, "(9x,' level   gdept    gdepw     e3t      e3w  ')" )
2210         write(numout,*)"bathy" ,"min e3t"
2211         do jk=1,jpk
2212         WRITE(numout, "(10x, i4, 4f9.2)" )  jk, gdept(20,20,jk), gdepw(20,20,jk), &
2213  &       e3t(20,20,jk), e3w(20,20,jk) 
2214         enddo
2215      ENDIF
2216
2217      mbathy(:,:)=0
2218!      WHERE( 0._wp < bathy(:,:)) mbathy(:,:)=jpkm1
2219      DO jj = 1, jpj
2220         DO ji = 1, jpi
2221             DO jk = 1, jpkm1
2222                IF( bathy(ji,jj) >= gdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
2223               
2224         END DO
2225        END DO
2226      END DO
2227
2228 !     DO jk = jpkm1, jpksigm+1, -1
2229 !        zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat )
2230 !        WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1
2231 !     END DO
2232
2233
2234
2235      ! z-partial steps :goto 20     
2236       DO jj = 1, jpj
2237         DO ji = 1, jpi
2238            ik = mbathy(ji,jj)
2239            IF( ik > jpksigm ) THEN               ! ocean point only
2240               ! max ocean level case
2241               IF( ik == jpkm1 ) THEN
2242                  zdepwp = bathy(ji,jj)
2243                  ze3tp  = bathy(ji,jj) - gdepw_0(ik)
2244                  ze3wp = 0.5_wp * e3w_0(ik) * ( 1._wp + ( ze3tp/e3t_0(ik) ) )
2245                  e3t(ji,jj,ik  ) = ze3tp
2246                  e3t(ji,jj,ik+1) = ze3tp
2247                  e3w(ji,jj,ik  ) = ze3wp
2248                  e3w(ji,jj,ik+1) = ze3tp
2249                  gdepw(ji,jj,ik+1) = zdepwp
2250                  gdept(ji,jj,ik  ) = gdept_0(ik-1) + ze3wp
2251                  gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp
2252                  !
2253               ELSE                         ! standard case
2254                  IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN   ;   gdepw(ji,jj,ik+1) = bathy(ji,jj)
2255                  ELSE                                       ;   gdepw(ji,jj,ik+1) = gdepw_0(ik+1)
2256                  ENDIF
2257!gm Bug?  check the gdepw_0
2258                  !       ... on ik
2259                  gdept(ji,jj,ik) = gdepw_0(ik) + ( gdepw  (ji,jj,ik+1) - gdepw_0(ik) )   &
2260                     &                          * ((gdept_0(      ik  ) - gdepw_0(ik) )   &
2261                     &                          / ( gdepw_0(      ik+1) - gdepw_0(ik) ))
2262                  e3t  (ji,jj,ik) = e3t_0  (ik) * ( gdepw  (ji,jj,ik+1) - gdepw_0(ik) )   & 
2263                     &                          / ( gdepw_0(      ik+1) - gdepw_0(ik) ) 
2264                  e3w  (ji,jj,ik) = 0.5_wp * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2._wp * gdepw_0(ik) )   &
2265                     &                     * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) )
2266                  !       ... on ik+1
2267                  e3w  (ji,jj,ik+1) = e3t  (ji,jj,ik)
2268                  e3t  (ji,jj,ik+1) = e3t  (ji,jj,ik)
2269                  gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)
2270               ENDIF
2271            ENDIF
2272         END DO
2273      END DO
2274      !
2275      jl = 0
2276      DO jj = 1, jpj
2277         DO ji = 1, jpi
2278            ik = mbathy(ji,jj)
2279            IF( ik > jpksigm ) THEN               ! ocean point only
2280               e3tp (ji,jj) = e3t(ji,jj,ik  )
2281               e3wp (ji,jj) = e3w(ji,jj,ik  )
2282               ! test
2283               zmin= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik  )
2284               IF( zmin <= 0._wp .AND. lwp ) THEN
2285                  jl = jl + 1
2286                  WRITE(numout,*) ' it      = ', jl, ' ik      = ', ik, ' (i,j) = ', ji, jj
2287                  WRITE(numout,*) ' bathy = ', bathy(ji,jj)
2288                  WRITE(numout,*) ' gdept = ', gdept(ji,jj,ik), ' gdepw = ', gdepw(ji,jj,ik+1), ' zdiff = ', zmin
2289                  WRITE(numout,*) ' e3tp  = ', e3t  (ji,jj,ik), ' e3wp  = ', e3w  (ji,jj,ik  )
2290               ENDIF
2291            ENDIF
2292         END DO
2293      END DO
2294
2295      ! Scale factors and depth at U-, V-, UW and VW-points
2296      DO jk = 1, jpk                        ! initialisation to z-scale factors
2297         e3u (:,:,jk) = e3t(:,:,jk)
2298         e3v (:,:,jk) = e3t(:,:,jk)
2299         e3uw(:,:,jk) = e3w(:,:,jk)
2300         e3vw(:,:,jk) = e3w(:,:,jk)
2301         e3f (:,:,jk) = e3t(:,:,jk)
2302      END DO
2303          DO jk = 1, jpksigm-1
2304           DO jj = 1, jpjm1
2305            DO ji = 1, fs_jpim1
2306                e3u(ji,jj,jk)=( REAL(MIN(1,mbathy(ji,jj)),wp)* e3t(ji,jj,jk)     +          &
2307                                REAL(MIN(1,mbathy(ji+1,jj)),wp)*e3t(ji+1,jj,jk) )           &
2308                /MAX(  1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj))   )
2309
2310                e3uw(ji,jj,jk)=(REAL(MIN(1,mbathy(ji,jj)),wp)* e3w(ji,jj,jk)    +          &
2311                                REAL(MIN(1,mbathy(ji+1,jj)),wp)*e3w(ji+1,jj,jk) )           &
2312                  /REAL(MAX(  1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj))),wp)
2313
2314                e3v(ji,jj,jk)=(REAL(MIN(1,mbathy(ji,jj)),wp)* e3t(ji,jj,jk)     +          &
2315                               REAL(MIN(1,mbathy(ji,jj+1)),wp)*e3t(ji,jj+1,jk) )           &
2316                /REAL (MAX(  1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1))),wp)
2317
2318                e3vw(ji,jj,jk)=(REAL(MIN(1,mbathy(ji,jj)),wp)* e3w(ji,jj,jk)    +          &
2319                                REAL(MIN(1,mbathy(ji,jj+1)),wp)*e3w(ji,jj+1,jk) )                 &
2320                /REAL(MAX(  1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1))),wp)
2321               
2322                e3f(ji,jj,jk)=(REAL(MIN(1,mbathy(ji,jj)),wp)* e3t(ji,jj,jk)     +          &
2323                               REAL(MIN(1,mbathy(ji+1,jj)),wp)*e3t(ji+1,jj,jk)  +          &
2324                               REAL(MIN(1,mbathy(ji+1,jj+1)),wp)* e3t(ji+1,jj+1,jk)+       &
2325                               REAL(MIN(1,mbathy(ji,jj+1)),wp)*e3t(ji,jj+1,jk) )           &
2326                /REAL(MAX(  1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1))             &
2327                      +   MIN(1,mbathy(ji+1,jj))+MIN(1,mbathy(ji+1,jj+1))),wp)
2328
2329               ENDDO
2330            END DO
2331         END DO
2332
2333      DO jk = jpksigm,jpk                         ! Computed as the minimum of neighbooring scale factors
2334               DO jj = 1, jpjm1
2335            DO ji = 1, fs_jpim1   ! vector opt.
2336               e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk) )
2337               e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk) )
2338               e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) )
2339               e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) )
2340            END DO
2341         END DO
2342      END DO
2343     
2344      CALL lbc_lnk( e3u , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw, 'U', 1._wp )   ! lateral boundary conditions
2345      CALL lbc_lnk( e3v , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw, 'V', 1._wp )
2346      !
2347      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
2348         WHERE( e3u (:,:,jk) == 0._wp )   e3u (:,:,jk) = e3t_0(jk)
2349         WHERE( e3v (:,:,jk) == 0._wp )   e3v (:,:,jk) = e3t_0(jk)
2350         WHERE( e3uw(:,:,jk) == 0._wp )   e3uw(:,:,jk) = e3w_0(jk)
2351         WHERE( e3vw(:,:,jk) == 0._wp )   e3vw(:,:,jk) = e3w_0(jk)
2352      END DO
2353
2354      DO jk = jpksigm, jpk                        ! Computed as the minimum of neighbooring V-scale factors
2355         DO jj = 1, jpjm1
2356            DO ji = 1, fs_jpim1   ! vector opt.
2357               e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) )
2358            END DO
2359         END DO
2360      END DO
2361      CALL lbc_lnk( e3f, 'F', 1._wp )       ! Lateral boundary conditions
2362      !
2363      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
2364         WHERE( e3f(:,:,jk) == 0._wp )   e3f(:,:,jk) = e3t_0(jk)
2365      END DO
2366!!gm  bug ? :  must be a do loop with mj0,mj1
2367      !
2368      e3t(:,mj0(1),:) = e3t(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2
2369      e3w(:,mj0(1),:) = e3w(:,mj0(2),:) 
2370      e3u(:,mj0(1),:) = e3u(:,mj0(2),:) 
2371      e3v(:,mj0(1),:) = e3v(:,mj0(2),:) 
2372      e3f(:,mj0(1),:) = e3f(:,mj0(2),:) 
2373
2374      ! Control of the sign
2375      Do jk=1,jpk
2376       do jj=1,jpj
2377        do ji=1,jpi
2378      IF( ( e3t  (ji,jj,jk) ) <= 0._wp )then
2379      write(numout,*)'    zgr_hyb :   e r r o r   e3t   <= 0',ji,jj,jk,e3t  (ji,jj,jk); endif
2380      IF( ( e3w  (ji,jj,jk) ) <= 0._wp )then
2381      write(numout,*)'    zgr_hyb :   e r r o r   e3t   <= 0',ji,jj,jk,e3w  (ji,jj,jk); endif
2382     
2383     
2384      IF( ( gdept(ji,jj,jk) ) <  0._wp )THEN
2385       write (numout,*)'   zgr_hyb :   e r r o r   gdept <  0',ji,jj,jk ,gdept(ji,jj,jj);    endif
2386      IF( ( gdepw(ji,jj,jk) ) <  0._wp )then
2387      write (numout,*)'   zgr_hyb :   e r r o r   gdepw <  0',ji,jj,jk , gdepw(ji,jj,jj);   endif
2388        enddo
2389        enddo
2390        enddo
2391     
2392      ! Compute gdep3w (vertical sum of e3w)
2393      gdep3w(:,:,1) = 0.5_wp * e3w(:,:,1)
2394      DO jk = 2, jpk
2395         gdep3w(:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk) 
2396      END DO
2397
2398
2399     IF( lk_mpp )   CALL mpp_max( nstop )
2400   IF (lwp) write(numout,*)"zpartial" ,nstop
2401
2402
2403
2404      CALL lbc_lnk( e3f, 'F', 1._wp )       ! Lateral boundary conditions
2405
2406     
2407
2408
2409      CALL wrk_dealloc( jpi, jpj,      zenv, ztmp, zmsk, zri, zrj, zhbat )
2410      !
2411      IF( nn_timing == 1 )  CALL timing_stop('zgr_hyb')
2412
2413   END SUBROUTINE zgr_hyb
2414
2415
2416   !!======================================================================
2417END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.