source: branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5098

Last change on this file since 5098 was 5098, checked in by mathiot, 6 years ago

add wmask, wumask, wvmask and restore loop order and add flag to ignore specific isf code if no isf

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