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

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

source: branches/UKMO/r5518_amm15_test/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 7144

Last change on this file since 7144 was 7144, checked in by jcastill, 7 years ago

Remove svn keywords

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