source: branches/UKMO/2015_CO6_CO5_zenv_wr_direct_dwl_temp/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5428

Last change on this file since 5428 was 5428, checked in by deazer, 6 years ago

Added incode changes to allow compoarisoin with CO5 after svn keywords removed.
Extracts, merges builds and runs as expected from working copy.

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