source: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5944

Last change on this file since 5944 was 5944, checked in by mathiot, 5 years ago

ISF: change related to reviewers comments

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