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 NEMO/branches/UKMO/r8395_cpl-pressure/NEMOGCM/TOOLS/DOMAINcfg/src – NEMO

source: NEMO/branches/UKMO/r8395_cpl-pressure/NEMOGCM/TOOLS/DOMAINcfg/src/domzgr.f90 @ 10797

Last change on this file since 10797 was 10797, checked in by jcastill, 5 years ago

Remove svn keywords

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