New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domzgr.F90 in branches/UKMO/dev_r5518_nemo2cice_prints/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/UKMO/dev_r5518_nemo2cice_prints/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 9816

Last change on this file since 9816 was 9816, checked in by dancopsey, 6 years ago

Strip out SVN keywords.

File size: 129.8 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(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) +              &
1873        &         bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
1874               zenv(ji,jj) = rn_sbot_min
1875             ENDIF
1876           ENDIF
1877         END DO
1878      END DO
1879      ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1880      CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1881      !
1882      ! smooth the bathymetry (if required)
1883      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
1884      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
1885      !
1886      jl = 0
1887      zrmax = 1._wp
1888      !   
1889      !     
1890      ! set scaling factor used in reducing vertical gradients
1891      zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
1892      !
1893      ! initialise temporary evelope depth arrays
1894      ztmpi1(:,:) = zenv(:,:)
1895      ztmpi2(:,:) = zenv(:,:)
1896      ztmpj1(:,:) = zenv(:,:)
1897      ztmpj2(:,:) = zenv(:,:)
1898      !
1899      ! initialise temporary r-value arrays
1900      zri(:,:) = 1._wp
1901      zrj(:,:) = 1._wp
1902      !                                                            ! ================ !
1903      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) !  Iterative loop  !
1904         !                                                         ! ================ !
1905         jl = jl + 1
1906         zrmax = 0._wp
1907         ! we set zrmax from previous r-values (zri and zrj) first
1908         ! if set after current r-value calculation (as previously)
1909         ! we could exit DO WHILE prematurely before checking r-value
1910         ! of current zenv
1911         DO jj = 1, nlcj
1912            DO ji = 1, nlci
1913               zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
1914            END DO
1915         END DO
1916         zri(:,:) = 0._wp
1917         zrj(:,:) = 0._wp
1918         DO jj = 1, nlcj
1919            DO ji = 1, nlci
1920               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
1921               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
1922               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
1923                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
1924               END IF
1925               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
1926                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
1927               END IF
1928               IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact
1929               IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact
1930               IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact
1931               IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact
1932            END DO
1933         END DO
1934         IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
1935         !
1936         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax
1937         !
1938         DO jj = 1, nlcj
1939            DO ji = 1, nlci
1940               zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
1941            END DO
1942         END DO
1943         ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1944         CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1945         !                                                  ! ================ !
1946      END DO                                                !     End loop     !
1947      !                                                     ! ================ !
1948      DO jj = 1, jpj
1949         DO ji = 1, jpi
1950            zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
1951         END DO
1952      END DO
1953      !
1954      ! Envelope bathymetry saved in hbatt
1955      hbatt(:,:) = zenv(:,:) 
1956      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
1957         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1958         DO jj = 1, jpj
1959            DO ji = 1, jpi
1960               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
1961               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
1962            END DO
1963         END DO
1964      ENDIF
1965      !
1966      IF(lwp) THEN                             ! Control print
1967         WRITE(numout,*)
1968         WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
1969         WRITE(numout,*)
1970         CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )
1971         IF( nprint == 1 )   THEN       
1972            WRITE(numout,*) ' bathy  MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
1973            WRITE(numout,*) ' hbatt  MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
1974         ENDIF
1975      ENDIF
1976
1977      !                                        ! ==============================
1978      !                                        !   hbatu, hbatv, hbatf fields
1979      !                                        ! ==============================
1980      IF(lwp) THEN
1981         WRITE(numout,*)
1982         WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
1983      ENDIF
1984      hbatu(:,:) = rn_sbot_min
1985      hbatv(:,:) = rn_sbot_min
1986      hbatf(:,:) = rn_sbot_min
1987      DO jj = 1, jpjm1
1988        DO ji = 1, jpim1   ! NO vector opt.
1989           hbatu(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
1990           hbatv(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
1991           hbatf(ji,jj) = 0.25_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
1992              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
1993        END DO
1994      END DO
1995      !
1996      ! Apply lateral boundary condition
1997!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
1998      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk( hbatu, 'U', 1._wp )
1999      DO jj = 1, jpj
2000         DO ji = 1, jpi
2001            IF( hbatu(ji,jj) == 0._wp ) THEN
2002               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
2003               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
2004            ENDIF
2005         END DO
2006      END DO
2007      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk( hbatv, 'V', 1._wp )
2008      DO jj = 1, jpj
2009         DO ji = 1, jpi
2010            IF( hbatv(ji,jj) == 0._wp ) THEN
2011               IF( zhbat(ji,jj) == 0._wp )   hbatv(ji,jj) = rn_sbot_min
2012               IF( zhbat(ji,jj) /= 0._wp )   hbatv(ji,jj) = zhbat(ji,jj)
2013            ENDIF
2014         END DO
2015      END DO
2016      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk( hbatf, 'F', 1._wp )
2017      DO jj = 1, jpj
2018         DO ji = 1, jpi
2019            IF( hbatf(ji,jj) == 0._wp ) THEN
2020               IF( zhbat(ji,jj) == 0._wp )   hbatf(ji,jj) = rn_sbot_min
2021               IF( zhbat(ji,jj) /= 0._wp )   hbatf(ji,jj) = zhbat(ji,jj)
2022            ENDIF
2023         END DO
2024      END DO
2025
2026!!bug:  key_helsinki a verifer
2027      hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
2028      hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
2029      hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
2030      hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
2031
2032      IF( nprint == 1 .AND. lwp )   THEN
2033         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
2034            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
2035         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
2036            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
2037         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
2038            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
2039         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
2040            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
2041      ENDIF
2042!! helsinki
2043
2044      !                                            ! =======================
2045      !                                            !   s-ccordinate fields     (gdep., e3.)
2046      !                                            ! =======================
2047      !
2048      ! non-dimensional "sigma" for model level depth at w- and t-levels
2049
2050
2051!========================================================================
2052! Song and Haidvogel  1994 (ln_s_sh94=T)
2053! Siddorn and Furner 2012 (ln_sf12=T)
2054! or  tanh function       (both false)                   
2055!========================================================================
2056      IF      ( ln_s_sh94 ) THEN
2057                           CALL s_sh94()
2058      ELSE IF ( ln_s_sf12 ) THEN
2059                           CALL s_sf12()
2060      ELSE                 
2061                           CALL s_tanh()
2062      ENDIF
2063
2064      CALL lbc_lnk( e3t_0 , 'T', 1._wp )
2065      CALL lbc_lnk( e3u_0 , 'U', 1._wp )
2066      CALL lbc_lnk( e3v_0 , 'V', 1._wp )
2067      CALL lbc_lnk( e3f_0 , 'F', 1._wp )
2068      CALL lbc_lnk( e3w_0 , 'W', 1._wp )
2069      CALL lbc_lnk( e3uw_0, 'U', 1._wp )
2070      CALL lbc_lnk( e3vw_0, 'V', 1._wp )
2071
2072      fsdepw(:,:,:) = gdepw_0 (:,:,:)
2073      fsde3w(:,:,:) = gdep3w_0(:,:,:)
2074      !
2075      where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0
2076      where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0
2077      where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0
2078      where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0
2079      where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0
2080      where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0
2081      where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0
2082
2083#if defined key_agrif
2084      ! Ensure meaningful vertical scale factors in ghost lines/columns
2085      IF( .NOT. Agrif_Root() ) THEN
2086         
2087         IF((nbondi == -1).OR.(nbondi == 2)) THEN
2088            e3u_0(1,:,:) = e3u_0(2,:,:)
2089         ENDIF
2090         !
2091         IF((nbondi ==  1).OR.(nbondi == 2)) THEN
2092            e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:)
2093         ENDIF
2094         !
2095         IF((nbondj == -1).OR.(nbondj == 2)) THEN
2096            e3v_0(:,1,:) = e3v_0(:,2,:)
2097         ENDIF
2098         !
2099         IF((nbondj ==  1).OR.(nbondj == 2)) THEN
2100            e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:)
2101         ENDIF
2102         !
2103      ENDIF
2104#endif
2105
2106      fsdept(:,:,:) = gdept_0 (:,:,:)
2107      fsdepw(:,:,:) = gdepw_0 (:,:,:)
2108      fsde3w(:,:,:) = gdep3w_0(:,:,:)
2109      fse3t (:,:,:) = e3t_0   (:,:,:)
2110      fse3u (:,:,:) = e3u_0   (:,:,:)
2111      fse3v (:,:,:) = e3v_0   (:,:,:)
2112      fse3f (:,:,:) = e3f_0   (:,:,:)
2113      fse3w (:,:,:) = e3w_0   (:,:,:)
2114      fse3uw(:,:,:) = e3uw_0  (:,:,:)
2115      fse3vw(:,:,:) = e3vw_0  (:,:,:)
2116!!
2117      ! HYBRID :
2118      DO jj = 1, jpj
2119         DO ji = 1, jpi
2120            DO jk = 1, jpkm1
2121               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
2122            END DO
2123            IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0
2124         END DO
2125      END DO
2126      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
2127         &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
2128
2129      IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain
2130         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) )
2131         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
2132            &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w '  , MINVAL( gdep3w_0(:,:,:) )
2133         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0   (:,:,:) ),   &
2134            &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0   (:,:,:) ),   &
2135            &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0  (:,:,:) ),   &
2136            &                          ' w ', MINVAL( e3w_0  (:,:,:) )
2137
2138         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
2139            &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w '  , MAXVAL( gdep3w_0(:,:,:) )
2140         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0   (:,:,:) ),   &
2141            &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0   (:,:,:) ),   &
2142            &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0  (:,:,:) ),   &
2143            &                          ' w ', MAXVAL( e3w_0  (:,:,:) )
2144      ENDIF
2145      !  END DO
2146      IF(lwp) THEN                                  ! selected vertical profiles
2147         WRITE(numout,*)
2148         WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
2149         WRITE(numout,*) ' ~~~~~~  --------------------'
2150         WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2151         WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     &
2152            &                                 e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
2153         DO jj = mj0(20), mj1(20)
2154            DO ji = mi0(20), mi1(20)
2155               WRITE(numout,*)
2156               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
2157               WRITE(numout,*) ' ~~~~~~  --------------------'
2158               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2159               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
2160                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2161            END DO
2162         END DO
2163         DO jj = mj0(74), mj1(74)
2164            DO ji = mi0(100), mi1(100)
2165               WRITE(numout,*)
2166               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
2167               WRITE(numout,*) ' ~~~~~~  --------------------'
2168               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2169               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
2170                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2171            END DO
2172         END DO
2173      ENDIF
2174
2175!================================================================================
2176! check the coordinate makes sense
2177!================================================================================
2178      DO ji = 1, jpi
2179         DO jj = 1, jpj
2180
2181            IF( hbatt(ji,jj) > 0._wp) THEN
2182               DO jk = 1, mbathy(ji,jj)
2183                 ! check coordinate is monotonically increasing
2184                 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
2185                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
2186                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
2187                    WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
2188                    WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
2189                    CALL ctl_stop( ctmp1 )
2190                 ENDIF
2191                 ! and check it has never gone negative
2192                 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
2193                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
2194                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk
2195                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2196                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2197                    CALL ctl_stop( ctmp1 )
2198                 ENDIF
2199                 ! and check it never exceeds the total depth
2200                 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
2201                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
2202                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
2203                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2204                    CALL ctl_stop( ctmp1 )
2205                 ENDIF
2206               END DO
2207
2208               DO jk = 1, mbathy(ji,jj)-1
2209                 ! and check it never exceeds the total depth
2210                IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
2211                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
2212                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
2213                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2214                    CALL ctl_stop( ctmp1 )
2215                 ENDIF
2216               END DO
2217
2218            ENDIF
2219
2220         END DO
2221      END DO
2222      !
2223      CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
2224      !
2225      IF( nn_timing == 1 )  CALL timing_stop('zgr_sco')
2226      !
2227   END SUBROUTINE zgr_sco
2228
2229!!======================================================================
2230   SUBROUTINE s_sh94()
2231
2232      !!----------------------------------------------------------------------
2233      !!                  ***  ROUTINE s_sh94  ***
2234      !!                     
2235      !! ** Purpose :   stretch the s-coordinate system
2236      !!
2237      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994
2238      !!                mixed S/sigma coordinate
2239      !!
2240      !! Reference : Song and Haidvogel 1994.
2241      !!----------------------------------------------------------------------
2242      !
2243      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2244      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
2245      !
2246      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2247      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
2248
2249      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2250      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2251
2252      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
2253      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
2254      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
2255      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
2256
2257      DO ji = 1, jpi
2258         DO jj = 1, jpj
2259
2260            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
2261               DO jk = 1, jpk
2262                  z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
2263                  z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
2264               END DO
2265            ELSE ! shallow water, uniform sigma
2266               DO jk = 1, jpk
2267                  z_gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
2268                  z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
2269                  END DO
2270            ENDIF
2271            !
2272            DO jk = 1, jpkm1
2273               z_esigt3(ji,jj,jk  ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2274               z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2275            END DO
2276            z_esigw3(ji,jj,1  ) = 2._wp * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) )
2277            z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
2278            !
2279            ! Coefficients for vertical depth as the sum of e3w scale factors
2280            z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
2281            DO jk = 2, jpk
2282               z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2283            END DO
2284            !
2285            DO jk = 1, jpk
2286               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2287               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2288               gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
2289               gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
2290               gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
2291            END DO
2292           !
2293         END DO   ! for all jj's
2294      END DO    ! for all ji's
2295
2296      DO ji = 1, jpim1
2297         DO jj = 1, jpjm1
2298            DO jk = 1, jpk
2299               z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
2300                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2301               z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
2302                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2303               z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     &
2304                  &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   &
2305                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2306               z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
2307                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2308               z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
2309                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2310               !
2311               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2312               e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2313               e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2314               e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2315               !
2316               e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2317               e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2318               e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2319            END DO
2320        END DO
2321      END DO
2322
2323      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2324      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2325
2326   END SUBROUTINE s_sh94
2327
2328   SUBROUTINE s_sf12
2329
2330      !!----------------------------------------------------------------------
2331      !!                  ***  ROUTINE s_sf12 ***
2332      !!                     
2333      !! ** Purpose :   stretch the s-coordinate system
2334      !!
2335      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
2336      !!                mixed S/sigma/Z coordinate
2337      !!
2338      !!                This method allows the maintenance of fixed surface and or
2339      !!                bottom cell resolutions (cf. geopotential coordinates)
2340      !!                within an analytically derived stretched S-coordinate framework.
2341      !!
2342      !!
2343      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
2344      !!----------------------------------------------------------------------
2345      !
2346      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2347      REAL(wp) ::   zsmth               ! smoothing around critical depth
2348      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
2349      !
2350      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2351      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
2352
2353      !
2354      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2355      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2356
2357      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
2358      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
2359      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
2360      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
2361
2362      DO ji = 1, jpi
2363         DO jj = 1, jpj
2364
2365          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
2366             
2367              zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
2368                                                     ! could be changed by users but care must be taken to do so carefully
2369              zzb = 1.0_wp-(zzb/hbatt(ji,jj))
2370           
2371              zzs = rn_zs / hbatt(ji,jj) 
2372             
2373              IF (rn_efold /= 0.0_wp) THEN
2374                zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
2375              ELSE
2376                zsmth = 1.0_wp 
2377              ENDIF
2378               
2379              DO jk = 1, jpk
2380                z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
2381                z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
2382              ENDDO
2383              z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth  )
2384              z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth  )
2385 
2386          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
2387
2388            DO jk = 1, jpk
2389              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
2390              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
2391            END DO
2392
2393          ELSE  ! shallow water, z coordinates
2394
2395            DO jk = 1, jpk
2396              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2397              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2398            END DO
2399
2400          ENDIF
2401
2402          DO jk = 1, jpkm1
2403             z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2404             z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2405          END DO
2406          z_esigw3(ji,jj,1  ) = 2.0_wp * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ))
2407          z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
2408
2409          ! Coefficients for vertical depth as the sum of e3w scale factors
2410          z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
2411          DO jk = 2, jpk
2412             z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2413          END DO
2414
2415          DO jk = 1, jpk
2416             gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
2417             gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
2418             gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
2419          END DO
2420
2421        ENDDO   ! for all jj's
2422      ENDDO    ! for all ji's
2423
2424      DO ji=1,jpi-1
2425        DO jj=1,jpj-1
2426
2427          DO jk = 1, jpk
2428                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / &
2429                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2430                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / &
2431                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2432                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  &
2433                                      hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / &
2434                                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2435                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / &
2436                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2437                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / &
2438                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2439
2440             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
2441             e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
2442             e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
2443             e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
2444             !
2445             e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
2446             e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
2447             e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
2448          END DO
2449
2450        ENDDO
2451      ENDDO
2452      !
2453      CALL lbc_lnk(e3t_0 ,'T',1.) ; CALL lbc_lnk(e3u_0 ,'T',1.)
2454      CALL lbc_lnk(e3v_0 ,'T',1.) ; CALL lbc_lnk(e3f_0 ,'T',1.)
2455      CALL lbc_lnk(e3w_0 ,'T',1.)
2456      CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.)
2457      !
2458      !                                               ! =============
2459
2460      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2461      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2462
2463   END SUBROUTINE s_sf12
2464
2465   SUBROUTINE s_tanh()
2466
2467      !!----------------------------------------------------------------------
2468      !!                  ***  ROUTINE s_tanh***
2469      !!                     
2470      !! ** Purpose :   stretch the s-coordinate system
2471      !!
2472      !! ** Method  :   s-coordinate stretch
2473      !!
2474      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
2475      !!----------------------------------------------------------------------
2476
2477      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2478      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
2479
2480      REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
2481      REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw
2482
2483      CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      )
2484      CALL wrk_alloc( jpk, z_esigt, z_esigw                                               )
2485
2486      z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp   ;   z_gsi3w  = 0._wp
2487      z_esigt  = 0._wp   ;   z_esigw  = 0._wp 
2488
2489      DO jk = 1, jpk
2490        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
2491        z_gsigt(jk) = -fssig( REAL(jk,wp)        )
2492      END DO
2493      IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk)
2494      !
2495      ! Coefficients for vertical scale factors at w-, t- levels
2496!!gm bug :  define it from analytical function, not like juste bellow....
2497!!gm        or betteroffer the 2 possibilities....
2498      DO jk = 1, jpkm1
2499         z_esigt(jk  ) = z_gsigw(jk+1) - z_gsigw(jk)
2500         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
2501      END DO
2502      z_esigw( 1 ) = 2._wp * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
2503      z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
2504      !
2505      ! Coefficients for vertical depth as the sum of e3w scale factors
2506      z_gsi3w(1) = 0.5_wp * z_esigw(1)
2507      DO jk = 2, jpk
2508         z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
2509      END DO
2510!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
2511      DO jk = 1, jpk
2512         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2513         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2514         gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
2515         gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
2516         gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )
2517      END DO
2518!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
2519      DO jj = 1, jpj
2520         DO ji = 1, jpi
2521            DO jk = 1, jpk
2522              e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2523              e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2524              e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2525              e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
2526              !
2527              e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2528              e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2529              e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2530            END DO
2531         END DO
2532      END DO
2533
2534      CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      )
2535      CALL wrk_dealloc( jpk, z_esigt, z_esigw                                               )
2536
2537   END SUBROUTINE s_tanh
2538
2539   FUNCTION fssig( pk ) RESULT( pf )
2540      !!----------------------------------------------------------------------
2541      !!                 ***  ROUTINE fssig ***
2542      !!       
2543      !! ** Purpose :   provide the analytical function in s-coordinate
2544      !!         
2545      !! ** Method  :   the function provide the non-dimensional position of
2546      !!                T and W (i.e. between 0 and 1)
2547      !!                T-points at integer values (between 1 and jpk)
2548      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2549      !!----------------------------------------------------------------------
2550      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
2551      REAL(wp)             ::   pf   ! sigma value
2552      !!----------------------------------------------------------------------
2553      !
2554      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   &
2555         &     - TANH( rn_thetb * rn_theta                                )  )   &
2556         & * (   COSH( rn_theta                           )                      &
2557         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              &
2558         & / ( 2._wp * SINH( rn_theta ) )
2559      !
2560   END FUNCTION fssig
2561
2562
2563   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
2564      !!----------------------------------------------------------------------
2565      !!                 ***  ROUTINE fssig1 ***
2566      !!
2567      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
2568      !!
2569      !! ** Method  :   the function provides the non-dimensional position of
2570      !!                T and W (i.e. between 0 and 1)
2571      !!                T-points at integer values (between 1 and jpk)
2572      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2573      !!----------------------------------------------------------------------
2574      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
2575      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
2576      REAL(wp)             ::   pf1   ! sigma value
2577      !!----------------------------------------------------------------------
2578      !
2579      IF ( rn_theta == 0 ) then      ! uniform sigma
2580         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
2581      ELSE                        ! stretched sigma
2582         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
2583            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
2584            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
2585      ENDIF
2586      !
2587   END FUNCTION fssig1
2588
2589
2590   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
2591      !!----------------------------------------------------------------------
2592      !!                 ***  ROUTINE fgamma  ***
2593      !!
2594      !! ** Purpose :   provide analytical function for the s-coordinate
2595      !!
2596      !! ** Method  :   the function provides the non-dimensional position of
2597      !!                T and W (i.e. between 0 and 1)
2598      !!                T-points at integer values (between 1 and jpk)
2599      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2600      !!
2601      !!                This method allows the maintenance of fixed surface and or
2602      !!                bottom cell resolutions (cf. geopotential coordinates)
2603      !!                within an analytically derived stretched S-coordinate framework.
2604      !!
2605      !! Reference  :   Siddorn and Furner, in prep
2606      !!----------------------------------------------------------------------
2607      REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate
2608      REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate
2609      REAL(wp), INTENT(in   ) ::   pzb           ! Bottom box depth
2610      REAL(wp), INTENT(in   ) ::   pzs           ! surface box depth
2611      REAL(wp), INTENT(in   ) ::   psmth       ! Smoothing parameter
2612      REAL(wp)                ::   za1,za2,za3    ! local variables
2613      REAL(wp)                ::   zn1,zn2        ! local variables
2614      REAL(wp)                ::   za,zb,zx       ! local variables
2615      integer                 ::   jk
2616      !!----------------------------------------------------------------------
2617      !
2618
2619      zn1  =  1./(jpk-1.)
2620      zn2  =  1. -  zn1
2621
2622      za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 
2623      za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp)
2624      za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
2625     
2626      za = pzb - za3*(pzs-za1)-za2
2627      za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) )
2628      zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
2629      zx = 1.0_wp-za/2.0_wp-zb
2630 
2631      DO jk = 1, jpk
2632        p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp +  &
2633                    & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &
2634                    &      (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
2635        p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth)
2636      ENDDO 
2637
2638      !
2639   END FUNCTION fgamma
2640
2641   !!======================================================================
2642END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.