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

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

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

Last change on this file since 9732 was 9732, checked in by mathiot, 4 years ago

replace the specific treatment of shallow ice shelves by a more general case

File size: 125.5 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
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_glhw_min   = 1.0e-3                             ! threshold on hw to define grounding line water
1301      REAL(wp) :: rn_isfhw_min  = 1.0e-3                            ! threshold on hw to define isf draft into the cavity
1302      REAL(wp) :: rn_isfshallow = 0.0_wp                            ! threshold to define shallow ice shelf cavity
1303      REAL(wp) :: rn_zisfmax    = 6000.0_wp                         ! maximun meter of ice we are allowed to dig to assure connectivity
1304      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)
1305      !!---------------------------------------------------------------------
1306      NAMELIST/namzgr_isf/nn_kisfmax, rn_zisfmax, &
1307              &           rn_isfdep_min, rn_glhw_min, rn_isfhw_min, &
1308              &           ln_isfcheminey, ln_isfconnect, ln_isfchannel, ln_isfsubgl, rn_isfsubgllon, rn_isfsubgllat
1309      !
1310      IF( nn_timing == 1 )  CALL timing_start('zgr_isf')
1311      !
1312      CALL wrk_alloc( jpi, jpj, zdummy, zmask, zrisfdep)
1313      CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy )
1314      !
1315      !
1316      REWIND( numnam_ref )              ! Namelist namzgr_isf in reference namelist : ice shelf geometry definition
1317      READ  ( numnam_ref, namzgr_isf, IOSTAT = ios, ERR = 901)
1318901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_isf in reference namelist', lwp )
1319
1320      REWIND( numnam_cfg )              ! Namelist namzgr_sco in configuration namelist : ice shelf geometry definition
1321      READ  ( numnam_cfg, namzgr_isf, IOSTAT = ios, ERR = 902 )
1322902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_isf in configuration namelist', lwp )
1323      IF(lwm) WRITE ( numond, namzgr_isf )
1324      !
1325      IF (lwp) THEN
1326         WRITE(numout,*) ' nn_kisfmax       = ',nn_kisfmax           !
1327         WRITE(numout,*) ' rn_zisfmax       = ',rn_zisfmax           !
1328         WRITE(numout,*) ' rn_isfdep_min    = ',rn_isfdep_min        !
1329         WRITE(numout,*) ' rn_isfhw_min     = ',rn_isfhw_min         !
1330         WRITE(numout,*) ' rn_glhw_min      = ',rn_glhw_min          !
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 => set the grounding
1342      ! line position
1343      WHERE( bathy(:,:) - risfdep(:,:) < rn_glhw_min ) 
1344         risfdep(:,:)=0.0 ; misfdep(:,:) = 1
1345         bathy  (:,:)=0.0 ; mbathy (:,:) = 0
1346      END WHERE
1347      !
1348      ! 1.2 ensure a minimum thickness for iceshelf cavity => avoid to negative
1349      ! e3t if ssh + sum(e3t*tmask) < 0
1350      DO jj = 1, jpj
1351         DO ji = 1, jpi
1352            IF (bathy(ji,jj)-risfdep(ji,jj) < rn_isfhw_min) THEN
1353               IF ((bathy(ji,jj) - rn_isfhw_min) > MAX(rn_isfdep_min,e3t_1d(1))) THEN
1354                  risfdep(ji,jj) = bathy(ji,jj) - rn_isfhw_min
1355               ELSE
1356                  risfdep(ji,jj)=0.0 ; misfdep(ji,jj) = 1
1357                  bathy  (ji,jj)=0.0 ; mbathy (ji,jj) = 0
1358               END IF
1359            END IF
1360         END DO
1361      END DO
1362      !
1363      ! 1.3 Remove channels and single point 'bay'.
1364      ! channel could be created if connectivity is enforced later.
1365      IF (ln_isfchannel) THEN
1366         zmask(:,:)=1._wp
1367         WHERE (mbathy(:,:)==0.0)
1368            zmask(:,:)=0._wp
1369         END WHERE
1370         DO jj = 2, jpjm1
1371            DO ji = 2, jpim1
1372               IF(  (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN
1373                  umskip1=zmask(ji,jj)*zmask(ji+1,jj  )  ! 1 = connexion
1374                  vmskjp1=zmask(ji,jj)*zmask(ji  ,jj+1)  ! 1 = connexion
1375                  umskim1=zmask(ji,jj)*zmask(ji-1,jj  )  ! 1 = connexion
1376                  vmskjm1=zmask(ji,jj)*zmask(ji  ,jj-1)  ! 1 = connexion
1377               ! zonal channel and single bay
1378                  IF ((umskip1+umskim1>=1) .AND. (vmskjp1+vmskjm1==0)) THEN
1379                     risfdep(ji,jj)=0.0 ; misfdep(ji,jj) = 1
1380                     bathy  (ji,jj)=0.0 ; mbathy (ji,jj) = 0
1381                  END IF
1382               ! meridionnal channel and single bay
1383                  IF ((vmskjp1+vmskjm1>=1) .AND. (umskip1+umskim1==0)) THEN
1384                     risfdep(ji,jj)=0.0 ; misfdep(ji,jj) = 1
1385                     bathy  (ji,jj)=0.0 ; mbathy (ji,jj) = 0
1386                  END IF
1387               END IF
1388            END DO
1389         END DO
1390         ! ensure halo correct
1391         IF( lk_mpp ) THEN
1392            zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1393            zdummy(:,:) = FLOAT( mbathy(:,:)  ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:)  = INT( zdummy(:,:) )
1394            CALL lbc_lnk( risfdep, 'T', 1. )
1395            CALL lbc_lnk( bathy  , 'T', 1. )
1396         ENDIF       
1397      END IF
1398      !
1399      ! 1.4 Closed single isolated point
1400      zmask(:,:)=1._wp
1401      WHERE (mbathy(:,:)==0.0)
1402         zmask(:,:)=0._wp
1403      END WHERE
1404      DO jj = 2, jpjm1
1405         DO ji = 2, jpim1
1406            IF(  (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN
1407               umskip1=zmask(ji,jj)*zmask(ji+1,jj  )  ! 1 = connexion
1408               vmskjp1=zmask(ji,jj)*zmask(ji  ,jj+1)  ! 1 = connexion
1409               umskim1=zmask(ji,jj)*zmask(ji-1,jj  )  ! 1 = connexion
1410               vmskjm1=zmask(ji,jj)*zmask(ji  ,jj-1)  ! 1 = connexion
1411               ! remove single point
1412               IF (umskip1+umskim1+vmskjp1+vmskjm1==0.0_wp) THEN
1413                  risfdep(ji,jj)=0.0 ; misfdep(ji,jj) = 1
1414                  bathy  (ji,jj)=0.0 ; mbathy (ji,jj) = 0
1415               END IF
1416            END IF
1417         END DO
1418      END DO
1419      ! ensure halo correct
1420      IF( lk_mpp ) THEN
1421         zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1422         zdummy(:,:) = FLOAT( mbathy(:,:)  ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:)  = INT( zdummy(:,:) )
1423         CALL lbc_lnk( risfdep, 'T', 1. )
1424         CALL lbc_lnk( bathy  , 'T', 1. )
1425      ENDIF
1426      !
1427      ! 2 Compute misfdep for ocean points (i.e. first wet level)
1428      ! find the first ocean level such that the first level thickness
1429      ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where
1430      ! e3t_0 is the reference level thickness
1431      !
1432      ! 2.0 set misfdep to 1 on open water and initialisation beneath ice shelf
1433      WHERE( risfdep(:,:) == 0.0 ) ;   misfdep(:,:) = 1   ! open water or grounded : set misfdep to 1 
1434      ELSEWHERE                    ;   misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level
1435      END WHERE 
1436      !
1437      DO jk = 2, jpkm1 
1438         zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
1439         WHERE( risfdep(:,:) > 0.0 .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1 
1440      END DO 
1441      !
1442      ! 2.1 be sure misfdep not below mbathy
1443      ! warning because of condition 4 we could have 'wet cell with misfdep below mbathy
1444      ! risfdep of these cells will be fix later on (see 5)
1445      WHERE( misfdep > mbathy .AND. mbathy > 0 ) misfdep=mbathy
1446      !
1447      ! 3.0 Assure 2 wet cells in the column at T point and along the edge.
1448      ! first do T, then loop 4 times on U-1, U, V-1, V and then at T point
1449      ! If un-luky, digging cell U-1 can trigger case for U+1, then V-1, then V+1
1450      !
1451      ! find the deepest isfdep level that fit the 2 wet cell on the water column
1452      ! on all the sides (still need 4 pass)
1453      DO jn = 1, 4
1454         zmisfdep = misfdep
1455         DO jj = 2, jpjm1
1456            DO ji = 2, jpim1
1457               ! ISF cell only
1458               IF(  (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN
1459                  zmbathyij  = mbathy(ji  ,jj)
1460                  ! set ground edge value to jpk to skip it later on
1461                  zmbathyip1 = mbathy(ji+1,jj) ; IF (zmbathyip1 < misfdep(ji,jj)) zmbathyip1 = jpk ! not wet cell in common on this edge
1462                  zmbathyim1 = mbathy(ji-1,jj) ; IF (zmbathyim1 < misfdep(ji,jj)) zmbathyim1 = jpk ! not wet cell in common on this edge
1463                  zmbathyjp1 = mbathy(ji,jj+1) ; IF (zmbathyjp1 < misfdep(ji,jj)) zmbathyjp1 = jpk ! not wet cell in common on this edge
1464                  zmbathyjm1 = mbathy(ji,jj-1) ; IF (zmbathyjm1 < misfdep(ji,jj)) zmbathyjm1 = jpk ! not wet cell in common on this edge
1465                  ! shallowest bathy level
1466                  zmbathyij = MIN(zmbathyij,zmbathyip1,zmbathyim1,zmbathyjp1,zmbathyjm1)
1467                  ! misfdep need to be <= zmbathyij-1 to fit 2 wet cell on the
1468                  ! water column
1469                  jk = MIN(misfdep(ji,jj),zmbathyij-1)
1470                  ! update isf draft if needed
1471                  IF ( jk < misfdep(ji,jj) ) THEN
1472                     zmisfdep(ji,jj) = jk
1473                     risfdep(ji,jj)  = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1474                  ELSE
1475                     ! isf depth not modify, nothing to do
1476                  END IF
1477               ENDIF
1478            END DO
1479         END DO
1480         misfdep=zmisfdep
1481         ! ensure halo correct before new pass
1482         IF( lk_mpp ) THEN
1483            zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1484            zdummy(:,:) = FLOAT( mbathy(:,:)  ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:)  = INT( zdummy(:,:) )
1485            CALL lbc_lnk( risfdep, 'T', 1. )
1486            CALL lbc_lnk( bathy  , 'T', 1. )
1487         ENDIF
1488      END DO ! jn
1489      !
1490      ! 3.1 condition block to assure connectivity every where beneath an ice shelf
1491      IF (ln_isfconnect) THEN
1492         zmask(:,:)=1.0
1493         WHERE (mbathy==0)
1494            zmask(:,:)=jpk
1495         END WHERE
1496
1497         zmbathy  = mbathy
1498         zmisfdep = misfdep
1499         zrisfdep = risfdep
1500         WHERE (mbathy == 0) zmbathy=jpk
1501         DO jj = 2, jpjm1
1502            DO ji = 2, jpim1
1503               IF(  (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN
1504! what it should be
1505                  umskip1=zmask(ji,jj)*zmask(ji+1,jj  )  ! 1 = should be connected
1506                  umskim1=zmask(ji,jj)*zmask(ji-1,jj  )  ! 1 = should be connected
1507                  vmskjp1=zmask(ji,jj)*zmask(ji  ,jj+1)  ! 1 = should be connected
1508                  vmskjm1=zmask(ji,jj)*zmask(ji  ,jj-1)  ! 1 = should be connected
1509! what it is
1510                  umskip1_r=jpk ; umskim1_r=jpk; vmskjp1_r=jpk; vmskjm1_r=jpk
1511                  IF (misfdep(ji,jj) > zmbathy(ji+1,jj  )) umskip1_r=1.0 ! 1 = no effective connection
1512                  IF (misfdep(ji,jj) > zmbathy(ji-1,jj  )) umskim1_r=1.0
1513                  IF (misfdep(ji,jj) > zmbathy(ji  ,jj+1)) vmskjp1_r=1.0
1514                  IF (misfdep(ji,jj) > zmbathy(ji  ,jj-1)) vmskjm1_r=1.0
1515                  ! defining level need for connectivity
1516                  jk=NINT(MIN(zmbathy(ji+1,jj  ) * umskip1_r * umskip1, &
1517                     &        zmbathy(ji-1,jj  ) * umskim1_r * umskim1, &
1518                     &        zmbathy(ji  ,jj+1) * vmskjp1_r * vmskjp1, &
1519                     &        zmbathy(ji  ,jj-1) * vmskjm1_r * vmskjm1, &
1520                     &        jpk+1.0 )) ! add jpk in the MIN to avoid out of boundary later on
1521! if connectivity is OK or no connection needed (grounding line) or grounded, zmisfdep=misfdep
1522                  zmisfdep(ji,jj)=MIN(misfdep(ji,jj),jk-1)
1523                  !
1524                  ! update isf draft if needed (need to be done here because next test is in meter and level)
1525                  IF ( zmisfdep(ji,jj) < misfdep(ji,jj) ) THEN
1526                     jk = zmisfdep(ji,jj)
1527                     zrisfdep(ji,jj)  = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1528                  END IF
1529   
1530                  ! sanity check
1531                  ! check if we dig more than nn_kisfmax level or reach the surface
1532                  ! check if we dig more than rn_zisfmax meter
1533                  ! => if this is the case, undo what has been done before
1534                  IF (      (misfdep(ji,jj)-zmisfdep(ji,jj) > MIN(nn_kisfmax,misfdep(ji,jj)-2)) &
1535                     & .OR. (risfdep(ji,jj)-zrisfdep(ji,jj) > MIN(rn_zisfmax,risfdep(ji,jj)  )) ) THEN
1536                     zmisfdep(ji,jj)=misfdep(ji,jj) 
1537                     zrisfdep(ji,jj)=risfdep(ji,jj)
1538                  END IF
1539               END IF
1540            END DO
1541         END DO
1542         misfdep=zmisfdep
1543         risfdep=zrisfdep
1544         ! ensure halo correct
1545         IF( lk_mpp ) THEN
1546            zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1547            zdummy(:,:) = FLOAT( mbathy(:,:)  ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:)  = INT( zdummy(:,:) )
1548            CALL lbc_lnk( risfdep, 'T', 1. )
1549            CALL lbc_lnk( bathy  , 'T', 1. )
1550         ENDIF
1551      END IF
1552      !
1553      IF (ln_isfcheminey) THEN
1554      ! 3.2 fill hole in ice shelf (ie cell with no velocity point)
1555      !      => misfdep = MIN(misfdep at U, U-1, V, V-1)
1556      !         risfdep = gdepw(misfdep) (ie full cell)
1557         zmisfdep = misfdep
1558         WHERE (mbathy == 0) zmisfdep=jpk   ! grounded
1559         DO jj = 2, jpjm1
1560            DO ji = 2, jpim1
1561               ibtest    = zmisfdep(ji  ,jj  )
1562               ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  )
1563               ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1)
1564               ! case unconected wet cell (T point) where isfdraft(ii,jj) deeper
1565               ! than bathy U-1/U+1/V-1/V+1
1566               IF( zmisfdep(ji,jj) > mbathy(ji-1,jj  ) ) ibtestim1 = jpk ! mask at U-1
1567               IF( zmisfdep(ji,jj) > mbathy(ji+1,jj  ) ) ibtestip1 = jpk ! mask at U+1
1568               IF( zmisfdep(ji,jj) > mbathy(ji  ,jj-1) ) ibtestjm1 = jpk ! mask at V-1
1569               IF( zmisfdep(ji,jj) > mbathy(ji  ,jj+1) ) ibtestjp1 = jpk ! mask at V+1
1570               ! case unconected wet cell (T point) where bathy(ii,jj) shallower
1571               ! than isfdraft U-1/U+1/V-1/V+1
1572               IF( mbathy(ji,jj) < zmisfdep(ji-1,jj  ) ) ibtestim1 = jpk ! mask at U-1
1573               IF( mbathy(ji,jj) < zmisfdep(ji+1,jj  ) ) ibtestip1 = jpk ! mask at U+1
1574               IF( mbathy(ji,jj) < zmisfdep(ji  ,jj-1) ) ibtestjm1 = jpk ! mask at V-1
1575               IF( mbathy(ji,jj) < zmisfdep(ji  ,jj+1) ) ibtestjp1 = jpk ! mask at V+1
1576               ! if surround in fully mask => mask cell  (1)
1577               ! if hole in the ice shelf, set to the min of surrounding (the MIN is doing the job)
1578               ! if misfdep is not changed, nothing is done (the MAX is doing the
1579               ! job)
1580               ibtest = MAX(ibtest,MIN(ibtestim1, ibtestip1, ibtestjm1,ibtestjp1))
1581               IF (misfdep(ji,jj) < ibtest) THEN
1582                  misfdep(ji,jj) = ibtest
1583                  risfdep(ji,jj) = gdepw_1d(ibtest)
1584               END IF
1585            ENDDO
1586         ENDDO
1587         WHERE( misfdep == jpk)   ! ground case (1)
1588            mbathy = 0 ; bathy = 0.0_wp ; misfdep = 1 ; risfdep = 0.0_wp
1589         END WHERE
1590         !
1591         ! ensure halo correct
1592         IF( lk_mpp ) THEN
1593            zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1594            zdummy(:,:) = FLOAT( mbathy(:,:)  ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:)  = INT( zdummy(:,:) )
1595            CALL lbc_lnk( risfdep, 'T', 1. )
1596            CALL lbc_lnk( bathy  , 'T', 1. )
1597         END IF
1598      END IF
1599      !
1600      CALL wrk_dealloc( jpi, jpj, zmask, zdummy, zrisfdep)
1601      CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy)
1602      !
1603      IF( nn_timing == 1 )  CALL timing_stop('zgr_isf')
1604     
1605   END SUBROUTINE
1606
1607   SUBROUTINE zgr_sco
1608      !!----------------------------------------------------------------------
1609      !!                  ***  ROUTINE zgr_sco  ***
1610      !!                     
1611      !! ** Purpose :   define the s-coordinate system
1612      !!
1613      !! ** Method  :   s-coordinate
1614      !!         The depth of model levels is defined as the product of an
1615      !!      analytical function by the local bathymetry, while the vertical
1616      !!      scale factors are defined as the product of the first derivative
1617      !!      of the analytical function by the bathymetry.
1618      !!      (this solution save memory as depth and scale factors are not
1619      !!      3d fields)
1620      !!          - Read bathymetry (in meters) at t-point and compute the
1621      !!         bathymetry at u-, v-, and f-points.
1622      !!            hbatu = mi( hbatt )
1623      !!            hbatv = mj( hbatt )
1624      !!            hbatf = mi( mj( hbatt ) )
1625      !!          - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
1626      !!         function and its derivative given as function.
1627      !!            z_gsigt(k) = fssig (k    )
1628      !!            z_gsigw(k) = fssig (k-0.5)
1629      !!            z_esigt(k) = fsdsig(k    )
1630      !!            z_esigw(k) = fsdsig(k-0.5)
1631      !!      Three options for stretching are give, and they can be modified
1632      !!      following the users requirements. Nevertheless, the output as
1633      !!      well as the way to compute the model levels and scale factors
1634      !!      must be respected in order to insure second order accuracy
1635      !!      schemes.
1636      !!
1637      !!      The three methods for stretching available are:
1638      !!
1639      !!           s_sh94 (Song and Haidvogel 1994)
1640      !!                a sinh/tanh function that allows sigma and stretched sigma
1641      !!
1642      !!           s_sf12 (Siddorn and Furner 2012?)
1643      !!                allows the maintenance of fixed surface and or
1644      !!                bottom cell resolutions (cf. geopotential coordinates)
1645      !!                within an analytically derived stretched S-coordinate framework.
1646      !!
1647      !!          s_tanh  (Madec et al 1996)
1648      !!                a cosh/tanh function that gives stretched coordinates       
1649      !!
1650      !!----------------------------------------------------------------------
1651      !
1652      INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
1653      INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
1654      INTEGER  ::   ios                      ! Local integer output status for namelist read
1655      REAL(wp) ::   zrmax, ztaper   ! temporary scalars
1656      REAL(wp) ::   zrfact
1657      !
1658      REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
1659      REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
1660
1661      NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
1662                           rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
1663     !!----------------------------------------------------------------------
1664      !
1665      IF( nn_timing == 1 )  CALL timing_start('zgr_sco')
1666      !
1667      CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
1668      !
1669      REWIND( numnam_ref )              ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
1670      READ  ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901)
1671901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist', lwp )
1672
1673      REWIND( numnam_cfg )              ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters
1674      READ  ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 )
1675902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist', lwp )
1676      IF(lwm) WRITE ( numond, namzgr_sco )
1677
1678      IF(lwp) THEN                           ! control print
1679         WRITE(numout,*)
1680         WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
1681         WRITE(numout,*) '~~~~~~~~~~~'
1682         WRITE(numout,*) '   Namelist namzgr_sco'
1683         WRITE(numout,*) '     stretching coeffs '
1684         WRITE(numout,*) '        maximum depth of s-bottom surface (>0)       rn_sbot_max   = ',rn_sbot_max
1685         WRITE(numout,*) '        minimum depth of s-bottom surface (>0)       rn_sbot_min   = ',rn_sbot_min
1686         WRITE(numout,*) '        Critical depth                               rn_hc         = ',rn_hc
1687         WRITE(numout,*) '        maximum cut-off r-value allowed              rn_rmax       = ',rn_rmax
1688         WRITE(numout,*) '     Song and Haidvogel 1994 stretching              ln_s_sh94     = ',ln_s_sh94
1689         WRITE(numout,*) '        Song and Haidvogel 1994 stretching coefficients'
1690         WRITE(numout,*) '        surface control parameter (0<=rn_theta<=20)  rn_theta      = ',rn_theta
1691         WRITE(numout,*) '        bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ',rn_thetb
1692         WRITE(numout,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb
1693         WRITE(numout,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12
1694         WRITE(numout,*) '        switching to sigma (T) or Z (F) at H<Hc      ln_sigcrit    = ',ln_sigcrit
1695         WRITE(numout,*) '        Siddorn and Furner 2012 stretching coefficients'
1696         WRITE(numout,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha
1697         WRITE(numout,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold
1698         WRITE(numout,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs
1699         WRITE(numout,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a
1700         WRITE(numout,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b
1701         WRITE(numout,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b'
1702      ENDIF
1703
1704      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate
1705      hifu(:,:) = rn_sbot_min
1706      hifv(:,:) = rn_sbot_min
1707      hiff(:,:) = rn_sbot_min
1708
1709      !                                        ! set maximum ocean depth
1710      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1711
1712      DO jj = 1, jpj
1713         DO ji = 1, jpi
1714           IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1715         END DO
1716      END DO
1717      !                                        ! =============================
1718      !                                        ! Define the envelop bathymetry   (hbatt)
1719      !                                        ! =============================
1720      ! use r-value to create hybrid coordinates
1721      zenv(:,:) = bathy(:,:)
1722      !
1723      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
1724      DO jj = 1, jpj
1725         DO ji = 1, jpi
1726           IF( bathy(ji,jj) == 0._wp ) THEN
1727             iip1 = MIN( ji+1, jpi )
1728             ijp1 = MIN( jj+1, jpj )
1729             iim1 = MAX( ji-1, 1 )
1730             ijm1 = MAX( jj-1, 1 )
1731             IF( ( + bathy(iim1,ijp1) + bathy(ji,ijp1) + bathy(iip1,ijp1)  &
1732                &  + bathy(iim1,jj  )                  + bathy(iip1,jj  )  &
1733                &  + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1)  ) > 0._wp ) THEN
1734                zenv(ji,jj) = rn_sbot_min
1735             ENDIF
1736           ENDIF
1737         END DO
1738      END DO
1739      ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1740      CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1741      !
1742      ! smooth the bathymetry (if required)
1743      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
1744      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
1745      !
1746      jl = 0
1747      zrmax = 1._wp
1748      !   
1749      !     
1750      ! set scaling factor used in reducing vertical gradients
1751      zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
1752      !
1753      ! initialise temporary evelope depth arrays
1754      ztmpi1(:,:) = zenv(:,:)
1755      ztmpi2(:,:) = zenv(:,:)
1756      ztmpj1(:,:) = zenv(:,:)
1757      ztmpj2(:,:) = zenv(:,:)
1758      !
1759      ! initialise temporary r-value arrays
1760      zri(:,:) = 1._wp
1761      zrj(:,:) = 1._wp
1762      !                                                            ! ================ !
1763      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) !  Iterative loop  !
1764         !                                                         ! ================ !
1765         jl = jl + 1
1766         zrmax = 0._wp
1767         ! we set zrmax from previous r-values (zri and zrj) first
1768         ! if set after current r-value calculation (as previously)
1769         ! we could exit DO WHILE prematurely before checking r-value
1770         ! of current zenv
1771         DO jj = 1, nlcj
1772            DO ji = 1, nlci
1773               zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
1774            END DO
1775         END DO
1776         zri(:,:) = 0._wp
1777         zrj(:,:) = 0._wp
1778         DO jj = 1, nlcj
1779            DO ji = 1, nlci
1780               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
1781               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
1782               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
1783                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
1784               END IF
1785               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
1786                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
1787               END IF
1788               IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact
1789               IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact
1790               IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact
1791               IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact
1792            END DO
1793         END DO
1794         IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
1795         !
1796         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax
1797         !
1798         DO jj = 1, nlcj
1799            DO ji = 1, nlci
1800               zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
1801            END DO
1802         END DO
1803         ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1804         CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1805         !                                                  ! ================ !
1806      END DO                                                !     End loop     !
1807      !                                                     ! ================ !
1808      DO jj = 1, jpj
1809         DO ji = 1, jpi
1810            zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
1811         END DO
1812      END DO
1813      !
1814      ! Envelope bathymetry saved in hbatt
1815      hbatt(:,:) = zenv(:,:) 
1816      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
1817         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1818         DO jj = 1, jpj
1819            DO ji = 1, jpi
1820               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
1821               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
1822            END DO
1823         END DO
1824      ENDIF
1825      !
1826      IF(lwp) THEN                             ! Control print
1827         WRITE(numout,*)
1828         WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
1829         WRITE(numout,*)
1830         CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )
1831         IF( nprint == 1 )   THEN       
1832            WRITE(numout,*) ' bathy  MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
1833            WRITE(numout,*) ' hbatt  MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
1834         ENDIF
1835      ENDIF
1836
1837      !                                        ! ==============================
1838      !                                        !   hbatu, hbatv, hbatf fields
1839      !                                        ! ==============================
1840      IF(lwp) THEN
1841         WRITE(numout,*)
1842         WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
1843      ENDIF
1844      hbatu(:,:) = rn_sbot_min
1845      hbatv(:,:) = rn_sbot_min
1846      hbatf(:,:) = rn_sbot_min
1847      DO jj = 1, jpjm1
1848        DO ji = 1, jpim1   ! NO vector opt.
1849           hbatu(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
1850           hbatv(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
1851           hbatf(ji,jj) = 0.25_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
1852              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
1853        END DO
1854      END DO
1855      !
1856      ! Apply lateral boundary condition
1857!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
1858      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk( hbatu, 'U', 1._wp )
1859      DO jj = 1, jpj
1860         DO ji = 1, jpi
1861            IF( hbatu(ji,jj) == 0._wp ) THEN
1862               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
1863               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
1864            ENDIF
1865         END DO
1866      END DO
1867      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk( hbatv, 'V', 1._wp )
1868      DO jj = 1, jpj
1869         DO ji = 1, jpi
1870            IF( hbatv(ji,jj) == 0._wp ) THEN
1871               IF( zhbat(ji,jj) == 0._wp )   hbatv(ji,jj) = rn_sbot_min
1872               IF( zhbat(ji,jj) /= 0._wp )   hbatv(ji,jj) = zhbat(ji,jj)
1873            ENDIF
1874         END DO
1875      END DO
1876      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk( hbatf, 'F', 1._wp )
1877      DO jj = 1, jpj
1878         DO ji = 1, jpi
1879            IF( hbatf(ji,jj) == 0._wp ) THEN
1880               IF( zhbat(ji,jj) == 0._wp )   hbatf(ji,jj) = rn_sbot_min
1881               IF( zhbat(ji,jj) /= 0._wp )   hbatf(ji,jj) = zhbat(ji,jj)
1882            ENDIF
1883         END DO
1884      END DO
1885
1886!!bug:  key_helsinki a verifer
1887      hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1888      hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1889      hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1890      hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1891
1892      IF( nprint == 1 .AND. lwp )   THEN
1893         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
1894            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
1895         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
1896            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
1897         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
1898            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1899         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
1900            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1901      ENDIF
1902!! helsinki
1903
1904      !                                            ! =======================
1905      !                                            !   s-ccordinate fields     (gdep., e3.)
1906      !                                            ! =======================
1907      !
1908      ! non-dimensional "sigma" for model level depth at w- and t-levels
1909
1910
1911!========================================================================
1912! Song and Haidvogel  1994 (ln_s_sh94=T)
1913! Siddorn and Furner 2012 (ln_sf12=T)
1914! or  tanh function       (both false)                   
1915!========================================================================
1916      IF      ( ln_s_sh94 ) THEN
1917                           CALL s_sh94()
1918      ELSE IF ( ln_s_sf12 ) THEN
1919                           CALL s_sf12()
1920      ELSE                 
1921                           CALL s_tanh()
1922      ENDIF
1923
1924      CALL lbc_lnk( e3t_0 , 'T', 1._wp )
1925      CALL lbc_lnk( e3u_0 , 'U', 1._wp )
1926      CALL lbc_lnk( e3v_0 , 'V', 1._wp )
1927      CALL lbc_lnk( e3f_0 , 'F', 1._wp )
1928      CALL lbc_lnk( e3w_0 , 'W', 1._wp )
1929      CALL lbc_lnk( e3uw_0, 'U', 1._wp )
1930      CALL lbc_lnk( e3vw_0, 'V', 1._wp )
1931
1932      fsdepw(:,:,:) = gdepw_0 (:,:,:)
1933      fsde3w(:,:,:) = gdep3w_0(:,:,:)
1934      !
1935      where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0
1936      where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0
1937      where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0
1938      where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0
1939      where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0
1940      where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0
1941      where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0
1942
1943#if defined key_agrif
1944      ! Ensure meaningful vertical scale factors in ghost lines/columns
1945      IF( .NOT. Agrif_Root() ) THEN
1946         
1947         IF((nbondi == -1).OR.(nbondi == 2)) THEN
1948            e3u_0(1,:,:) = e3u_0(2,:,:)
1949         ENDIF
1950         !
1951         IF((nbondi ==  1).OR.(nbondi == 2)) THEN
1952            e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:)
1953         ENDIF
1954         !
1955         IF((nbondj == -1).OR.(nbondj == 2)) THEN
1956            e3v_0(:,1,:) = e3v_0(:,2,:)
1957         ENDIF
1958         !
1959         IF((nbondj ==  1).OR.(nbondj == 2)) THEN
1960            e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:)
1961         ENDIF
1962         !
1963      ENDIF
1964#endif
1965
1966      fsdept(:,:,:) = gdept_0 (:,:,:)
1967      fsdepw(:,:,:) = gdepw_0 (:,:,:)
1968      fsde3w(:,:,:) = gdep3w_0(:,:,:)
1969      fse3t (:,:,:) = e3t_0   (:,:,:)
1970      fse3u (:,:,:) = e3u_0   (:,:,:)
1971      fse3v (:,:,:) = e3v_0   (:,:,:)
1972      fse3f (:,:,:) = e3f_0   (:,:,:)
1973      fse3w (:,:,:) = e3w_0   (:,:,:)
1974      fse3uw(:,:,:) = e3uw_0  (:,:,:)
1975      fse3vw(:,:,:) = e3vw_0  (:,:,:)
1976!!
1977      ! HYBRID :
1978      DO jj = 1, jpj
1979         DO ji = 1, jpi
1980            DO jk = 1, jpkm1
1981               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
1982            END DO
1983            IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0
1984         END DO
1985      END DO
1986      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
1987         &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
1988
1989      IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain
1990         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) )
1991         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
1992            &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w '  , MINVAL( gdep3w_0(:,:,:) )
1993         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0   (:,:,:) ),   &
1994            &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0   (:,:,:) ),   &
1995            &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0  (:,:,:) ),   &
1996            &                          ' w ', MINVAL( e3w_0  (:,:,:) )
1997
1998         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
1999            &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w '  , MAXVAL( gdep3w_0(:,:,:) )
2000         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0   (:,:,:) ),   &
2001            &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0   (:,:,:) ),   &
2002            &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0  (:,:,:) ),   &
2003            &                          ' w ', MAXVAL( e3w_0  (:,:,:) )
2004      ENDIF
2005      !  END DO
2006      IF(lwp) THEN                                  ! selected vertical profiles
2007         WRITE(numout,*)
2008         WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
2009         WRITE(numout,*) ' ~~~~~~  --------------------'
2010         WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2011         WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     &
2012            &                                 e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
2013         DO jj = mj0(20), mj1(20)
2014            DO ji = mi0(20), mi1(20)
2015               WRITE(numout,*)
2016               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,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         DO jj = mj0(74), mj1(74)
2024            DO ji = mi0(100), mi1(100)
2025               WRITE(numout,*)
2026               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
2027               WRITE(numout,*) ' ~~~~~~  --------------------'
2028               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2029               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
2030                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2031            END DO
2032         END DO
2033      ENDIF
2034
2035!================================================================================
2036! check the coordinate makes sense
2037!================================================================================
2038      DO ji = 1, jpi
2039         DO jj = 1, jpj
2040
2041            IF( hbatt(ji,jj) > 0._wp) THEN
2042               DO jk = 1, mbathy(ji,jj)
2043                 ! check coordinate is monotonically increasing
2044                 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
2045                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
2046                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
2047                    WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
2048                    WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
2049                    CALL ctl_stop( ctmp1 )
2050                 ENDIF
2051                 ! and check it has never gone negative
2052                 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
2053                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
2054                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk
2055                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2056                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2057                    CALL ctl_stop( ctmp1 )
2058                 ENDIF
2059                 ! and check it never exceeds the total depth
2060                 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
2061                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
2062                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
2063                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2064                    CALL ctl_stop( ctmp1 )
2065                 ENDIF
2066               END DO
2067
2068               DO jk = 1, mbathy(ji,jj)-1
2069                 ! and check it never exceeds the total depth
2070                IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
2071                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
2072                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
2073                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2074                    CALL ctl_stop( ctmp1 )
2075                 ENDIF
2076               END DO
2077
2078            ENDIF
2079
2080         END DO
2081      END DO
2082      !
2083      CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
2084      !
2085      IF( nn_timing == 1 )  CALL timing_stop('zgr_sco')
2086      !
2087   END SUBROUTINE zgr_sco
2088
2089!!======================================================================
2090   SUBROUTINE s_sh94()
2091
2092      !!----------------------------------------------------------------------
2093      !!                  ***  ROUTINE s_sh94  ***
2094      !!                     
2095      !! ** Purpose :   stretch the s-coordinate system
2096      !!
2097      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994
2098      !!                mixed S/sigma coordinate
2099      !!
2100      !! Reference : Song and Haidvogel 1994.
2101      !!----------------------------------------------------------------------
2102      !
2103      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2104      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
2105      !
2106      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2107      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
2108
2109      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2110      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2111
2112      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
2113      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
2114      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
2115      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
2116
2117      DO ji = 1, jpi
2118         DO jj = 1, jpj
2119
2120            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
2121               DO jk = 1, jpk
2122                  z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
2123                  z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
2124               END DO
2125            ELSE ! shallow water, uniform sigma
2126               DO jk = 1, jpk
2127                  z_gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
2128                  z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
2129                  END DO
2130            ENDIF
2131            !
2132            DO jk = 1, jpkm1
2133               z_esigt3(ji,jj,jk  ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2134               z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2135            END DO
2136            z_esigw3(ji,jj,1  ) = 2._wp * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) )
2137            z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
2138            !
2139            ! Coefficients for vertical depth as the sum of e3w scale factors
2140            z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
2141            DO jk = 2, jpk
2142               z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2143            END DO
2144            !
2145            DO jk = 1, jpk
2146               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2147               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2148               gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
2149               gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
2150               gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
2151            END DO
2152           !
2153         END DO   ! for all jj's
2154      END DO    ! for all ji's
2155
2156      DO ji = 1, jpim1
2157         DO jj = 1, jpjm1
2158            DO jk = 1, jpk
2159               z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
2160                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2161               z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
2162                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2163               z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     &
2164                  &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   &
2165                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2166               z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
2167                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2168               z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
2169                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2170               !
2171               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2172               e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2173               e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2174               e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2175               !
2176               e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2177               e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2178               e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2179            END DO
2180        END DO
2181      END DO
2182
2183      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2184      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2185
2186   END SUBROUTINE s_sh94
2187
2188   SUBROUTINE s_sf12
2189
2190      !!----------------------------------------------------------------------
2191      !!                  ***  ROUTINE s_sf12 ***
2192      !!                     
2193      !! ** Purpose :   stretch the s-coordinate system
2194      !!
2195      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
2196      !!                mixed S/sigma/Z coordinate
2197      !!
2198      !!                This method allows the maintenance of fixed surface and or
2199      !!                bottom cell resolutions (cf. geopotential coordinates)
2200      !!                within an analytically derived stretched S-coordinate framework.
2201      !!
2202      !!
2203      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
2204      !!----------------------------------------------------------------------
2205      !
2206      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2207      REAL(wp) ::   zsmth               ! smoothing around critical depth
2208      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
2209      !
2210      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2211      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
2212
2213      !
2214      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2215      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2216
2217      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
2218      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
2219      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
2220      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
2221
2222      DO ji = 1, jpi
2223         DO jj = 1, jpj
2224
2225          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
2226             
2227              zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
2228                                                     ! could be changed by users but care must be taken to do so carefully
2229              zzb = 1.0_wp-(zzb/hbatt(ji,jj))
2230           
2231              zzs = rn_zs / hbatt(ji,jj) 
2232             
2233              IF (rn_efold /= 0.0_wp) THEN
2234                zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
2235              ELSE
2236                zsmth = 1.0_wp 
2237              ENDIF
2238               
2239              DO jk = 1, jpk
2240                z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
2241                z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
2242              ENDDO
2243              z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth  )
2244              z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth  )
2245 
2246          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
2247
2248            DO jk = 1, jpk
2249              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
2250              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
2251            END DO
2252
2253          ELSE  ! shallow water, z coordinates
2254
2255            DO jk = 1, jpk
2256              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2257              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2258            END DO
2259
2260          ENDIF
2261
2262          DO jk = 1, jpkm1
2263             z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2264             z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2265          END DO
2266          z_esigw3(ji,jj,1  ) = 2.0_wp * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ))
2267          z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
2268
2269          ! Coefficients for vertical depth as the sum of e3w scale factors
2270          z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
2271          DO jk = 2, jpk
2272             z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2273          END DO
2274
2275          DO jk = 1, jpk
2276             gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
2277             gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
2278             gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
2279          END DO
2280
2281        ENDDO   ! for all jj's
2282      ENDDO    ! for all ji's
2283
2284      DO ji=1,jpi-1
2285        DO jj=1,jpj-1
2286
2287          DO jk = 1, jpk
2288                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / &
2289                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2290                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / &
2291                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2292                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  &
2293                                      hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / &
2294                                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2295                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / &
2296                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2297                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / &
2298                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2299
2300             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
2301             e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
2302             e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
2303             e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
2304             !
2305             e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
2306             e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
2307             e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
2308          END DO
2309
2310        ENDDO
2311      ENDDO
2312      !
2313      CALL lbc_lnk(e3t_0 ,'T',1.) ; CALL lbc_lnk(e3u_0 ,'T',1.)
2314      CALL lbc_lnk(e3v_0 ,'T',1.) ; CALL lbc_lnk(e3f_0 ,'T',1.)
2315      CALL lbc_lnk(e3w_0 ,'T',1.)
2316      CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.)
2317      !
2318      !                                               ! =============
2319
2320      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2321      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2322
2323   END SUBROUTINE s_sf12
2324
2325   SUBROUTINE s_tanh()
2326
2327      !!----------------------------------------------------------------------
2328      !!                  ***  ROUTINE s_tanh***
2329      !!                     
2330      !! ** Purpose :   stretch the s-coordinate system
2331      !!
2332      !! ** Method  :   s-coordinate stretch
2333      !!
2334      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
2335      !!----------------------------------------------------------------------
2336
2337      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2338      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
2339
2340      REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
2341      REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw
2342
2343      CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      )
2344      CALL wrk_alloc( jpk, z_esigt, z_esigw                                               )
2345
2346      z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp   ;   z_gsi3w  = 0._wp
2347      z_esigt  = 0._wp   ;   z_esigw  = 0._wp 
2348
2349      DO jk = 1, jpk
2350        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
2351        z_gsigt(jk) = -fssig( REAL(jk,wp)        )
2352      END DO
2353      IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk)
2354      !
2355      ! Coefficients for vertical scale factors at w-, t- levels
2356!!gm bug :  define it from analytical function, not like juste bellow....
2357!!gm        or betteroffer the 2 possibilities....
2358      DO jk = 1, jpkm1
2359         z_esigt(jk  ) = z_gsigw(jk+1) - z_gsigw(jk)
2360         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
2361      END DO
2362      z_esigw( 1 ) = 2._wp * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
2363      z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
2364      !
2365      ! Coefficients for vertical depth as the sum of e3w scale factors
2366      z_gsi3w(1) = 0.5_wp * z_esigw(1)
2367      DO jk = 2, jpk
2368         z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
2369      END DO
2370!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
2371      DO jk = 1, jpk
2372         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2373         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2374         gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
2375         gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
2376         gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )
2377      END DO
2378!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
2379      DO jj = 1, jpj
2380         DO ji = 1, jpi
2381            DO jk = 1, jpk
2382              e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2383              e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2384              e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2385              e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
2386              !
2387              e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2388              e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2389              e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2390            END DO
2391         END DO
2392      END DO
2393
2394      CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      )
2395      CALL wrk_dealloc( jpk, z_esigt, z_esigw                                               )
2396
2397   END SUBROUTINE s_tanh
2398
2399   FUNCTION fssig( pk ) RESULT( pf )
2400      !!----------------------------------------------------------------------
2401      !!                 ***  ROUTINE fssig ***
2402      !!       
2403      !! ** Purpose :   provide the analytical function in s-coordinate
2404      !!         
2405      !! ** Method  :   the function provide the non-dimensional position of
2406      !!                T and W (i.e. between 0 and 1)
2407      !!                T-points at integer values (between 1 and jpk)
2408      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2409      !!----------------------------------------------------------------------
2410      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
2411      REAL(wp)             ::   pf   ! sigma value
2412      !!----------------------------------------------------------------------
2413      !
2414      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   &
2415         &     - TANH( rn_thetb * rn_theta                                )  )   &
2416         & * (   COSH( rn_theta                           )                      &
2417         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              &
2418         & / ( 2._wp * SINH( rn_theta ) )
2419      !
2420   END FUNCTION fssig
2421
2422
2423   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
2424      !!----------------------------------------------------------------------
2425      !!                 ***  ROUTINE fssig1 ***
2426      !!
2427      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
2428      !!
2429      !! ** Method  :   the function provides the non-dimensional position of
2430      !!                T and W (i.e. between 0 and 1)
2431      !!                T-points at integer values (between 1 and jpk)
2432      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2433      !!----------------------------------------------------------------------
2434      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
2435      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
2436      REAL(wp)             ::   pf1   ! sigma value
2437      !!----------------------------------------------------------------------
2438      !
2439      IF ( rn_theta == 0 ) then      ! uniform sigma
2440         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
2441      ELSE                        ! stretched sigma
2442         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
2443            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
2444            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
2445      ENDIF
2446      !
2447   END FUNCTION fssig1
2448
2449
2450   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
2451      !!----------------------------------------------------------------------
2452      !!                 ***  ROUTINE fgamma  ***
2453      !!
2454      !! ** Purpose :   provide analytical function for the s-coordinate
2455      !!
2456      !! ** Method  :   the function provides the non-dimensional position of
2457      !!                T and W (i.e. between 0 and 1)
2458      !!                T-points at integer values (between 1 and jpk)
2459      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2460      !!
2461      !!                This method allows the maintenance of fixed surface and or
2462      !!                bottom cell resolutions (cf. geopotential coordinates)
2463      !!                within an analytically derived stretched S-coordinate framework.
2464      !!
2465      !! Reference  :   Siddorn and Furner, in prep
2466      !!----------------------------------------------------------------------
2467      REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate
2468      REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate
2469      REAL(wp), INTENT(in   ) ::   pzb           ! Bottom box depth
2470      REAL(wp), INTENT(in   ) ::   pzs           ! surface box depth
2471      REAL(wp), INTENT(in   ) ::   psmth       ! Smoothing parameter
2472      REAL(wp)                ::   za1,za2,za3    ! local variables
2473      REAL(wp)                ::   zn1,zn2        ! local variables
2474      REAL(wp)                ::   za,zb,zx       ! local variables
2475      integer                 ::   jk
2476      !!----------------------------------------------------------------------
2477      !
2478
2479      zn1  =  1./(jpk-1.)
2480      zn2  =  1. -  zn1
2481
2482      za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 
2483      za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp)
2484      za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
2485     
2486      za = pzb - za3*(pzs-za1)-za2
2487      za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) )
2488      zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
2489      zx = 1.0_wp-za/2.0_wp-zb
2490 
2491      DO jk = 1, jpk
2492        p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp +  &
2493                    & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &
2494                    &      (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
2495        p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth)
2496      ENDDO 
2497
2498      !
2499   END FUNCTION fgamma
2500
2501   !!======================================================================
2502END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.