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

Last change on this file since 9316 was 9316, checked in by mathiot, 2 years ago

update domzgr and namelist (rewrite of zgr_isf based on discussion with P. Holland, R. Smith, A. Siaahan and J. Gregory)

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