source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5443

Last change on this file since 5443 was 5443, checked in by davestorkey, 5 years ago

Update 2015/dev_r5021_UKMO1_CICE_coupling branch to revision 5442 of the trunk.

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