source: branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 9513

Last change on this file since 9513 was 9513, checked in by mathiot, 3 years ago

Add option to detect and remove subglacial lake (do not affect closed sea option)

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