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

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

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

Last change on this file since 5120 was 5120, checked in by acc, 9 years ago

#1473. Re-organisation and optimisation of ice shelf cavity option. This commit merges changes from the dev_r5094_UKMO_ISFCLEAN branch onto the trunk. Results will change, even with ln_isfcav=F, due to a return to original definitions of the vertical metrics. All changes have been reviewed and SETTE tested.

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