source: branches/2015/dev_r5918_nemo_v3_6_STABLE_XIOS2/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 6064

Last change on this file since 6064 was 6064, checked in by cetlod, 5 years ago

dev_xios2 : update to the head of v3.6 stable branch

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