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

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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5015

Last change on this file since 5015 was 5015, checked in by timgraham, 9 years ago

Fix for unitialised arrays in domzgr.F90 (#1431)

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