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

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

source: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 4747

Last change on this file since 4747 was 4747, checked in by mathiot, 10 years ago

ISF branch: change to deal with non mask bathymetry (land processor definition, building bathy and ice shelf draft variable), update of hpg (definition of ze3wu in case of zps and vvl) and bfr (in case of 2 cell water column thickness, each cell feels top and bottom friction).

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