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

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

source: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5952

Last change on this file since 5952 was 5952, checked in by mathiot, 8 years ago

ice sheet coupling: minor changes before merge

  • Property svn:keywords set to Id
File size: 134.8 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
531               ! set grounded point to 0
532               ! (a treshold could be set here if needed, or set it offline based on the grounded fraction)
533               WHERE ( bathy(:,:) <= risfdep(:,:) + rn_isfhmin )
534                  misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp
535                  mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp
536               END WHERE
537            END IF
538            !       
539            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
540               !
541              IF( nn_cla == 0 ) THEN
542                 ii0 = 140   ;   ii1 = 140                   ! Gibraltar Strait open
543                 ij0 = 102   ;   ij1 = 102                   ! (Thomson, Ocean Modelling, 1995)
544                 DO ji = mi0(ii0), mi1(ii1)
545                    DO jj = mj0(ij0), mj1(ij1)
546                       bathy(ji,jj) = 284._wp
547                    END DO
548                 END DO
549                 IF(lwp) WRITE(numout,*)     
550                 IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
551                 !
552                 ii0 = 160   ;   ii1 = 160                   ! Bab el mandeb Strait open
553                 ij0 = 88    ;   ij1 = 88                    ! (Thomson, Ocean Modelling, 1995)
554                 DO ji = mi0(ii0), mi1(ii1)
555                    DO jj = mj0(ij0), mj1(ij1)
556                       bathy(ji,jj) = 137._wp
557                    END DO
558                 END DO
559                 IF(lwp) WRITE(numout,*)
560                 IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
561              ENDIF
562              !
563           ENDIF
564            !
565        ENDIF
566         !                                            ! =============== !
567      ELSE                                            !      error      !
568         !                                            ! =============== !
569         WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
570         CALL ctl_stop( '    zgr_bat : '//trim(ctmp1) )
571      ENDIF
572      !
573      IF( nn_closea == 0 )   CALL clo_bat( bathy, mbathy )    !==  NO closed seas or lakes  ==!
574      !                       
575      IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==!
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      !! ** Action  : - test compatibility between isfdraft and bathy
1208      !!              - bathy and isfdraft are modified
1209      !!----------------------------------------------------------------------
1210      !!   
1211      INTEGER  ::   ji, jj, jl, jk       ! dummy loop indices
1212      INTEGER  ::   ik, it               ! temporary integers
1213      INTEGER  ::   icompt, ibtest       ! (ISF)
1214      INTEGER  ::   ibtestim1, ibtestip1 ! (ISF)
1215      INTEGER  ::   ibtestjm1, ibtestjp1 ! (ISF)
1216      REAL(wp) ::   zdepth           ! Ajusted ocean depth to avoid too small e3t
1217      REAL(wp) ::   zmax             ! Maximum and minimum depth
1218      REAL(wp) ::   zbathydiff       ! isf temporary scalar
1219      REAL(wp) ::   zrisfdepdiff     ! isf temporary scalar
1220      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
1221      REAL(wp) ::   zdepwp           ! Ajusted ocean depth to avoid too small e3t
1222      REAL(wp) ::   zdiff            ! temporary scalar
1223      REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 2D workspace (ISH)
1224      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep         ! 2D workspace (ISH)
1225      !!---------------------------------------------------------------------
1226      !
1227      IF( nn_timing == 1 )  CALL timing_start('zgr_isf')
1228      !
1229      CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep)
1230      CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy )
1231
1232      ! (ISF) compute misfdep
1233      WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) /= 0 ) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1 
1234      ELSEWHERE                      ;                         misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level
1235      END WHERE 
1236
1237      ! Compute misfdep for ocean points (i.e. first wet level)
1238      ! find the first ocean level such that the first level thickness
1239      ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where
1240      ! e3t_0 is the reference level thickness
1241      DO jk = 2, jpkm1 
1242         zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
1243         WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1 
1244      END DO
1245      WHERE ( 0._wp < risfdep(:,:) .AND. risfdep(:,:) <= e3t_1d(1) )
1246         risfdep(:,:) = 0. ; misfdep(:,:) = 1
1247      END WHERE
1248
1249! 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
1250      icompt = 0 
1251! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together
1252      DO jl = 1, 10     
1253         ! check at each iteration if isf is grounded or not (1cm treshold have to be update after first coupling experiments)
1254         WHERE (bathy(:,:) <= risfdep(:,:) + rn_isfhmin)
1255            misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp
1256            mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp
1257         END WHERE
1258         WHERE (mbathy(:,:) <= 0) 
1259            misfdep(:,:) = 0; risfdep(:,:) = 0._wp 
1260            mbathy (:,:) = 0; bathy  (:,:) = 0._wp
1261         ENDWHERE
1262         IF( lk_mpp ) THEN
1263            zbathy(:,:)  = FLOAT( misfdep(:,:) )
1264            CALL lbc_lnk( zbathy, 'T', 1. )
1265            misfdep(:,:) = INT( zbathy(:,:) )
1266
1267            CALL lbc_lnk( risfdep,'T', 1. )
1268            CALL lbc_lnk( bathy,  'T', 1. )
1269
1270            zbathy(:,:)  = FLOAT( mbathy(:,:) )
1271            CALL lbc_lnk( zbathy, 'T', 1. )
1272            mbathy(:,:)  = INT( zbathy(:,:) )
1273         ENDIF
1274         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
1275            misfdep( 1 ,:) = misfdep(jpim1,:)            ! local domain is cyclic east-west
1276            misfdep(jpi,:) = misfdep(  2  ,:) 
1277            mbathy( 1 ,:)  = mbathy(jpim1,:)             ! local domain is cyclic east-west
1278            mbathy(jpi,:)  = mbathy(  2  ,:)
1279         ENDIF
1280
1281         ! split last cell if possible (only where water column is 2 cell or less)
1282         ! if coupled to ice sheet, we do not modify the bathymetry (can be discuss).
1283         IF ( .NOT. ln_iscpl) THEN
1284            DO jk = jpkm1, 1, -1
1285               zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1286               WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)
1287                  mbathy(:,:) = jk
1288                  bathy(:,:)  = zmax
1289               END WHERE
1290            END DO
1291         END IF
1292 
1293         ! split top cell if possible (only where water column is 2 cell or less)
1294         DO jk = 2, jpkm1
1295            zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1296            WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy)
1297               misfdep(:,:) = jk
1298               risfdep(:,:) = zmax
1299            END WHERE
1300         END DO
1301
1302 
1303 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition
1304         DO jj = 1, jpj
1305            DO ji = 1, jpi
1306               ! find the minimum change option:
1307               ! test bathy
1308               IF (risfdep(ji,jj) > 1) THEN
1309                  IF ( .NOT. ln_iscpl ) THEN
1310                     zbathydiff  =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) &
1311                         &            + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
1312                     zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) &
1313                         &            - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
1314 
1315                     IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) <  misfdep(ji,jj)) THEN
1316                        IF (zbathydiff <= zrisfdepdiff) THEN
1317                           bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )
1318                           mbathy(ji,jj)= mbathy(ji,jj) + 1
1319                        ELSE
1320                           risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )
1321                           misfdep(ji,jj) = misfdep(ji,jj) - 1
1322                        END IF
1323                     ENDIF
1324                  ELSE
1325                     IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) <  misfdep(ji,jj)) THEN
1326                        risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )
1327                        misfdep(ji,jj) = misfdep(ji,jj) - 1
1328                     END IF
1329                  END IF
1330               END IF
1331            END DO
1332         END DO
1333 
1334         ! At least 2 levels for water thickness at T, U, and V point.
1335         DO jj = 1, jpj
1336            DO ji = 1, jpi
1337               ! find the minimum change option:
1338               ! test bathy
1339               IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN
1340                  IF ( .NOT. ln_iscpl ) 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                  ELSE
1353                     misfdep(ji,jj)= misfdep(ji,jj) - 1
1354                     risfdep(ji,jj)= gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )
1355                  END IF
1356               ENDIF
1357            END DO
1358         END DO
1359 
1360 ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1)
1361         DO jj = 1, jpjm1
1362            DO ji = 1, jpim1
1363               IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN
1364                  IF ( .NOT. ln_iscpl ) THEN
1365                     zbathydiff  =ABS(bathy(ji,jj    ) - ( gdepw_1d(mbathy (ji,jj)+1) &
1366                          &                              + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat )))
1367                     zrisfdepdiff=ABS(risfdep(ji,jj+1) - ( gdepw_1d(misfdep(ji,jj+1)) &
1368                          &                              - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))
1369                     IF (zbathydiff <= zrisfdepdiff) THEN
1370                        mbathy(ji,jj) = mbathy(ji,jj) + 1
1371                        bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj   )+1)*e3zps_rat )
1372                     ELSE
1373                        misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1
1374                        risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )
1375                     END IF
1376                  ELSE
1377                     misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1
1378                     risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )
1379                  END IF
1380               ENDIF
1381            END DO
1382         END DO
1383 
1384         IF( lk_mpp ) THEN
1385            zbathy(:,:)  = FLOAT( misfdep(:,:) )
1386            CALL lbc_lnk( zbathy, 'T', 1. )
1387            misfdep(:,:) = INT( zbathy(:,:) )
1388
1389            CALL lbc_lnk( risfdep,'T', 1. )
1390            CALL lbc_lnk( bathy,  'T', 1. )
1391
1392            zbathy(:,:)  = FLOAT( mbathy(:,:) )
1393            CALL lbc_lnk( zbathy, 'T', 1. )
1394            mbathy(:,:)  = INT( zbathy(:,:) )
1395         ENDIF
1396 ! point V misdep(ji,jj) EQ mbathy(ji,jj+1)
1397         DO jj = 1, jpjm1
1398            DO ji = 1, jpim1
1399               IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) > 1) THEN
1400                  IF ( .NOT. ln_iscpl ) THEN
1401                     zbathydiff  =ABS(  bathy(ji,jj+1) - ( gdepw_1d(mbathy (ji,jj+1)+1) &
1402                           &                             + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))
1403                     zrisfdepdiff=ABS(risfdep(ji,jj  ) - ( gdepw_1d(misfdep(ji,jj  )  ) &
1404                           &                             - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat )))
1405                     IF (zbathydiff <= zrisfdepdiff) THEN
1406                        mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1
1407                        bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )
1408                     ELSE
1409                        misfdep(ji,jj)   = misfdep(ji,jj) - 1
1410                        risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat )
1411                     END IF
1412                  ELSE
1413                     misfdep(ji,jj)   = misfdep(ji,jj) - 1
1414                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat )
1415                  END IF
1416               ENDIF
1417            END DO
1418         END DO
1419 
1420 
1421         IF( lk_mpp ) THEN
1422            zbathy(:,:)  = FLOAT( misfdep(:,:) )
1423            CALL lbc_lnk( zbathy, 'T', 1. )
1424            misfdep(:,:) = INT( zbathy(:,:) )
1425
1426            CALL lbc_lnk( risfdep,'T', 1. )
1427            CALL lbc_lnk( bathy,  'T', 1. )
1428
1429            zbathy(:,:)  = FLOAT( mbathy(:,:) )
1430            CALL lbc_lnk( zbathy, 'T', 1. )
1431            mbathy(:,:)  = INT( zbathy(:,:) )
1432         ENDIF 
1433 
1434 ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1)
1435         DO jj = 1, jpjm1
1436            DO ji = 1, jpim1
1437               IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN
1438                  IF ( .NOT. ln_iscpl ) THEN
1439                  zbathydiff  =ABS(  bathy(ji  ,jj) - ( gdepw_1d(mbathy (ji,jj)+1) &
1440                       &                              + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat )))
1441                  zrisfdepdiff=ABS(risfdep(ji+1,jj) - ( gdepw_1d(misfdep(ji+1,jj)) &
1442                       &                              - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))
1443                  IF (zbathydiff <= zrisfdepdiff) THEN
1444                     mbathy(ji,jj) = mbathy(ji,jj) + 1
1445                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
1446                  ELSE
1447                     misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1
1448                     risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )
1449                  END IF
1450                  ELSE
1451                     misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1
1452                     risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )
1453                  ENDIF
1454               ENDIF
1455            ENDDO
1456         ENDDO
1457 
1458         IF( lk_mpp ) THEN
1459            zbathy(:,:)  = FLOAT( misfdep(:,:) )
1460            CALL lbc_lnk( zbathy, 'T', 1. )
1461            misfdep(:,:) = INT( zbathy(:,:) )
1462
1463            CALL lbc_lnk( risfdep,'T', 1. )
1464            CALL lbc_lnk( bathy,  'T', 1. )
1465
1466            zbathy(:,:)  = FLOAT( mbathy(:,:) )
1467            CALL lbc_lnk( zbathy, 'T', 1. )
1468            mbathy(:,:)  = INT( zbathy(:,:) )
1469         ENDIF 
1470 
1471 ! point U misfdep(ji,jj) EQ bathy(ji,jj+1)
1472         DO jj = 1, jpjm1
1473            DO ji = 1, jpim1
1474               IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) > 1) THEN
1475                  IF ( .NOT. ln_iscpl ) THEN
1476                     zbathydiff  =ABS(  bathy(ji+1,jj) - ( gdepw_1d(mbathy (ji+1,jj)+1) &
1477                          &                              + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))
1478                     zrisfdepdiff=ABS(risfdep(ji  ,jj) - ( gdepw_1d(misfdep(ji  ,jj)  ) &
1479                          &                              - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat )))
1480                     IF (zbathydiff <= zrisfdepdiff) THEN
1481                        mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1
1482                        bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat )
1483                     ELSE
1484                        misfdep(ji,jj)   = misfdep(ji  ,jj) - 1
1485                        risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat )
1486                     END IF
1487                  ELSE
1488                     misfdep(ji,jj)   = misfdep(ji  ,jj) - 1
1489                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat )
1490                  ENDIF
1491               ENDIF
1492            ENDDO
1493         ENDDO
1494 
1495         IF( lk_mpp ) THEN
1496            zbathy(:,:)  = FLOAT( misfdep(:,:) )
1497            CALL lbc_lnk( zbathy, 'T', 1. )
1498            misfdep(:,:) = INT( zbathy(:,:) )
1499
1500            CALL lbc_lnk( risfdep,'T', 1. )
1501            CALL lbc_lnk( bathy,  'T', 1. )
1502
1503            zbathy(:,:)  = FLOAT( mbathy(:,:) )
1504            CALL lbc_lnk( zbathy, 'T', 1. )
1505            mbathy(:,:)  = INT( zbathy(:,:) )
1506         ENDIF
1507      END DO
1508      ! end dig bathy/ice shelf to be compatible
1509      ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness
1510      DO jl = 1,20
1511 
1512 ! remove single point "bay" on isf coast line in the ice shelf draft'
1513         DO jk = 2, jpk
1514            WHERE (misfdep==0) misfdep=jpk
1515            zmask=0._wp
1516            WHERE (misfdep <= jk) zmask=1
1517            DO jj = 2, jpjm1
1518               DO ji = 2, jpim1
1519                  IF (misfdep(ji,jj) == jk) THEN
1520                     ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
1521                     IF (ibtest <= 1) THEN
1522                        risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1
1523                        IF (misfdep(ji,jj) > mbathy(ji,jj)) misfdep(ji,jj) = jpk
1524                     END IF
1525                  END IF
1526               END DO
1527            END DO
1528         END DO
1529         WHERE (misfdep==jpk)
1530             misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp
1531         END WHERE
1532         IF( lk_mpp ) THEN
1533            zbathy(:,:)  = FLOAT( misfdep(:,:) )
1534            CALL lbc_lnk( zbathy, 'T', 1. )
1535            misfdep(:,:) = INT( zbathy(:,:) )
1536
1537            CALL lbc_lnk( risfdep,'T', 1. )
1538            CALL lbc_lnk( bathy,  'T', 1. )
1539
1540            zbathy(:,:)  = FLOAT( mbathy(:,:) )
1541            CALL lbc_lnk( zbathy, 'T', 1. )
1542            mbathy(:,:)  = INT( zbathy(:,:) )
1543         ENDIF
1544 
1545 ! remove single point "bay" on bathy coast line beneath an ice shelf'
1546         DO jk = jpk,1,-1
1547            zmask=0._wp
1548            WHERE (mbathy >= jk ) zmask=1
1549            DO jj = 2, jpjm1
1550               DO ji = 2, jpim1
1551                  IF (mbathy(ji,jj) == jk .AND. misfdep(ji,jj) >= 2) THEN
1552                     ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
1553                     IF (ibtest <= 1) THEN
1554                        bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1
1555                        IF (misfdep(ji,jj) > mbathy(ji,jj)) mbathy(ji,jj) = 0
1556                     END IF
1557                  END IF
1558               END DO
1559            END DO
1560         END DO
1561         WHERE (mbathy==0)
1562             misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp
1563         END WHERE
1564         IF( lk_mpp ) THEN
1565            zbathy(:,:)  = FLOAT( misfdep(:,:) )
1566            CALL lbc_lnk( zbathy, 'T', 1. )
1567            misfdep(:,:) = INT( zbathy(:,:) )
1568
1569            CALL lbc_lnk( risfdep,'T', 1. )
1570            CALL lbc_lnk( bathy,  'T', 1. )
1571
1572            zbathy(:,:)  = FLOAT( mbathy(:,:) )
1573            CALL lbc_lnk( zbathy, 'T', 1. )
1574            mbathy(:,:)  = INT( zbathy(:,:) )
1575         ENDIF
1576 
1577 ! fill hole in ice shelf
1578         zmisfdep = misfdep
1579         zrisfdep = risfdep
1580         WHERE (zmisfdep <= 1._wp) zmisfdep=jpk
1581         DO jj = 2, jpjm1
1582            DO ji = 2, jpim1
1583               ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  )
1584               ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1)
1585               IF( zmisfdep(ji,jj) >= mbathy(ji-1,jj  ) ) ibtestim1 = jpk
1586               IF( zmisfdep(ji,jj) >= mbathy(ji+1,jj  ) ) ibtestip1 = jpk
1587               IF( zmisfdep(ji,jj) >= mbathy(ji  ,jj-1) ) ibtestjm1 = jpk
1588               IF( zmisfdep(ji,jj) >= mbathy(ji  ,jj+1) ) ibtestjp1 = jpk
1589               ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
1590               IF( ibtest == jpk .AND. misfdep(ji,jj) >= 2) THEN
1591                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp
1592               END IF
1593               IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) >= 2) THEN
1594                  misfdep(ji,jj) = ibtest
1595                  risfdep(ji,jj) = gdepw_1d(ibtest)
1596               ENDIF
1597            ENDDO
1598         ENDDO
1599 
1600         IF( lk_mpp ) THEN
1601            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
1602            CALL lbc_lnk( zbathy,  'T', 1. ) 
1603            misfdep(:,:) = INT( zbathy(:,:) ) 
1604
1605            CALL lbc_lnk( risfdep, 'T', 1. ) 
1606            CALL lbc_lnk( bathy,   'T', 1. )
1607
1608            zbathy(:,:) = FLOAT( mbathy(:,:) )
1609            CALL lbc_lnk( zbathy,  'T', 1. )
1610            mbathy(:,:) = INT( zbathy(:,:) )
1611         ENDIF 
1612 !
1613 !! fill hole in bathymetry
1614         zmbathy (:,:)=mbathy (:,:)
1615         DO jj = 2, jpjm1
1616            DO ji = 2, jpim1
1617               ibtestim1 = zmbathy(ji-1,jj  ) ; ibtestip1 = zmbathy(ji+1,jj  )
1618               ibtestjm1 = zmbathy(ji  ,jj-1) ; ibtestjp1 = zmbathy(ji  ,jj+1)
1619               IF( zmbathy(ji,jj) <  misfdep(ji-1,jj  ) ) ibtestim1 = 0
1620               IF( zmbathy(ji,jj) <  misfdep(ji+1,jj  ) ) ibtestip1 = 0
1621               IF( zmbathy(ji,jj) <  misfdep(ji  ,jj-1) ) ibtestjm1 = 0
1622               IF( zmbathy(ji,jj) <  misfdep(ji  ,jj+1) ) ibtestjp1 = 0
1623               ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
1624               IF( ibtest == 0 .AND. misfdep(ji,jj) >= 2) THEN
1625                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ;
1626               END IF
1627               IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) >= 2) THEN
1628                  mbathy(ji,jj) = ibtest
1629                  bathy(ji,jj)  = gdepw_1d(ibtest+1) 
1630               ENDIF
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 (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN
1649                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,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 U point water column less than 2 cells), mask U
1666         DO jj = 1, jpjm1
1667            DO ji = 1, jpim1
1668               IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN
1669                  mbathy(ji+1,jj)  = mbathy(ji+1,jj) - 1;   bathy(ji+1,jj)   = gdepw_1d(mbathy(ji+1,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 (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN
1689                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+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 (ie V point water column less than 2 cells), mask V
1706         DO jj = 1, jpjm1
1707            DO ji = 1, jpi
1708               IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN
1709                  mbathy(ji,jj+1)  = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ;
1710               END IF
1711            END DO
1712         END DO
1713         IF( lk_mpp ) THEN
1714            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
1715            CALL lbc_lnk( zbathy, 'T', 1. ) 
1716            misfdep(:,:) = INT( zbathy(:,:) ) 
1717
1718            CALL lbc_lnk( risfdep,'T', 1. ) 
1719            CALL lbc_lnk( bathy,  'T', 1. )
1720
1721            zbathy(:,:) = FLOAT( mbathy(:,:) )
1722            CALL lbc_lnk( zbathy, 'T', 1. )
1723            mbathy(:,:) = INT( zbathy(:,:) )
1724         ENDIF 
1725 ! if not compatible after all check, mask T
1726         DO jj = 1, jpj
1727            DO ji = 1, jpi
1728               IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN
1729                  misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj)  = 0 ; bathy(ji,jj)   = 0._wp ;
1730               END IF
1731            END DO
1732         END DO
1733 
1734         WHERE (mbathy(:,:) == 1)
1735            mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp
1736         END WHERE
1737      END DO 
1738! end check compatibility ice shelf/bathy
1739      ! remove very shallow ice shelf (less than ~ 10m if 75L)
1740      IF( icompt == 0 ) THEN
1741         IF(lwp) WRITE(numout,*)'     no points with ice shelf too close to bathymetry' 
1742      ELSE
1743         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry' 
1744      ENDIF 
1745
1746      ! compute scale factor and depth at T- and W- points
1747      DO jj = 1, jpj
1748         DO ji = 1, jpi
1749            ik = mbathy(ji,jj)
1750            IF( ik > 0 ) THEN               ! ocean point only
1751               ! max ocean level case
1752               IF( ik == jpkm1 ) THEN
1753                  zdepwp = bathy(ji,jj)
1754                  ze3tp  = bathy(ji,jj) - gdepw_1d(ik)
1755                  ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
1756                  e3t_0(ji,jj,ik  ) = ze3tp
1757                  e3t_0(ji,jj,ik+1) = ze3tp
1758                  e3w_0(ji,jj,ik  ) = ze3wp
1759                  e3w_0(ji,jj,ik+1) = ze3tp
1760                  gdepw_0(ji,jj,ik+1) = zdepwp
1761                  gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp
1762                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
1763                  !
1764               ELSE                         ! standard case
1765                  IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
1766                  ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
1767                  ENDIF
1768      !            gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
1769!gm Bug?  check the gdepw_1d
1770                  !       ... on ik
1771                  gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   &
1772                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   &
1773                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) ))
1774                  e3t_0  (ji,jj,ik  ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik  )
1775                  e3w_0  (ji,jj,ik  ) = gdept_0(ji,jj,ik  ) - gdept_1d(ik-1)
1776                  !       ... on ik+1
1777                  e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
1778                  e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
1779               ENDIF
1780            ENDIF
1781         END DO
1782      END DO
1783      !
1784      it = 0
1785      DO jj = 1, jpj
1786         DO ji = 1, jpi
1787            ik = mbathy(ji,jj)
1788            IF( ik > 0 ) THEN               ! ocean point only
1789               e3tp (ji,jj) = e3t_0(ji,jj,ik)
1790               e3wp (ji,jj) = e3w_0(ji,jj,ik)
1791               ! test
1792               zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  )
1793               IF( zdiff <= 0._wp .AND. lwp ) THEN
1794                  it = it + 1
1795                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj
1796                  WRITE(numout,*) ' bathy = ', bathy(ji,jj)
1797                  WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
1798                  WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  )
1799               ENDIF
1800            ENDIF
1801         END DO
1802      END DO
1803      !
1804      ! (ISF) Definition of e3t, u, v, w for ISF case
1805      DO jj = 1, jpj 
1806         DO ji = 1, jpi 
1807            ik = misfdep(ji,jj) 
1808            IF( ik > 1 ) THEN               ! ice shelf point only
1809               IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik) 
1810               gdepw_0(ji,jj,ik) = risfdep(ji,jj) 
1811!gm Bug?  check the gdepw_0
1812            !       ... on ik
1813               gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   & 
1814                  &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   & 
1815                  &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      ) 
1816               e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 
1817               e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik)
1818
1819               IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)
1820                  e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 
1821               ENDIF 
1822            !       ... on ik / ik-1
1823               e3w_0  (ji,jj,ik  ) = e3t_0  (ji,jj,ik) 
1824               e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
1825! The next line isn't required and doesn't affect results - included for consistency with bathymetry code
1826               gdept_0(ji,jj,ik-1) = gdept_1d(ik-1)
1827            ENDIF
1828         END DO
1829      END DO
1830   
1831      it = 0 
1832      DO jj = 1, jpj 
1833         DO ji = 1, jpi 
1834            ik = misfdep(ji,jj) 
1835            IF( ik > 1 ) THEN               ! ice shelf point only
1836               e3tp (ji,jj) = e3t_0(ji,jj,ik  ) 
1837               e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 
1838            ! test
1839               zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  ) 
1840               IF( zdiff <= 0. .AND. lwp ) THEN 
1841                  it = it + 1 
1842                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
1843                  WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 
1844                  WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
1845                  WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj) 
1846               ENDIF
1847            ENDIF
1848         END DO
1849      END DO
1850
1851      CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep )
1852      CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy )
1853
1854      IF( nn_timing == 1 )  CALL timing_stop('zgr_isf')
1855     
1856   END SUBROUTINE zgr_isf
1857
1858   SUBROUTINE zgr_sco
1859      !!----------------------------------------------------------------------
1860      !!                  ***  ROUTINE zgr_sco  ***
1861      !!                     
1862      !! ** Purpose :   define the s-coordinate system
1863      !!
1864      !! ** Method  :   s-coordinate
1865      !!         The depth of model levels is defined as the product of an
1866      !!      analytical function by the local bathymetry, while the vertical
1867      !!      scale factors are defined as the product of the first derivative
1868      !!      of the analytical function by the bathymetry.
1869      !!      (this solution save memory as depth and scale factors are not
1870      !!      3d fields)
1871      !!          - Read bathymetry (in meters) at t-point and compute the
1872      !!         bathymetry at u-, v-, and f-points.
1873      !!            hbatu = mi( hbatt )
1874      !!            hbatv = mj( hbatt )
1875      !!            hbatf = mi( mj( hbatt ) )
1876      !!          - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
1877      !!         function and its derivative given as function.
1878      !!            z_gsigt(k) = fssig (k    )
1879      !!            z_gsigw(k) = fssig (k-0.5)
1880      !!            z_esigt(k) = fsdsig(k    )
1881      !!            z_esigw(k) = fsdsig(k-0.5)
1882      !!      Three options for stretching are give, and they can be modified
1883      !!      following the users requirements. Nevertheless, the output as
1884      !!      well as the way to compute the model levels and scale factors
1885      !!      must be respected in order to insure second order accuracy
1886      !!      schemes.
1887      !!
1888      !!      The three methods for stretching available are:
1889      !!
1890      !!           s_sh94 (Song and Haidvogel 1994)
1891      !!                a sinh/tanh function that allows sigma and stretched sigma
1892      !!
1893      !!           s_sf12 (Siddorn and Furner 2012?)
1894      !!                allows the maintenance of fixed surface and or
1895      !!                bottom cell resolutions (cf. geopotential coordinates)
1896      !!                within an analytically derived stretched S-coordinate framework.
1897      !!
1898      !!          s_tanh  (Madec et al 1996)
1899      !!                a cosh/tanh function that gives stretched coordinates       
1900      !!
1901      !!----------------------------------------------------------------------
1902      !
1903      INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
1904      INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
1905      INTEGER  ::   ios                      ! Local integer output status for namelist read
1906      REAL(wp) ::   zrmax, ztaper   ! temporary scalars
1907      REAL(wp) ::   zrfact
1908      !
1909      REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
1910      REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
1911
1912      NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
1913                           rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
1914     !!----------------------------------------------------------------------
1915      !
1916      IF( nn_timing == 1 )  CALL timing_start('zgr_sco')
1917      !
1918      CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
1919      !
1920      REWIND( numnam_ref )              ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
1921      READ  ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901)
1922901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist', lwp )
1923
1924      REWIND( numnam_cfg )              ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters
1925      READ  ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 )
1926902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist', lwp )
1927      IF(lwm) WRITE ( numond, namzgr_sco )
1928
1929      IF(lwp) THEN                           ! control print
1930         WRITE(numout,*)
1931         WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
1932         WRITE(numout,*) '~~~~~~~~~~~'
1933         WRITE(numout,*) '   Namelist namzgr_sco'
1934         WRITE(numout,*) '     stretching coeffs '
1935         WRITE(numout,*) '        maximum depth of s-bottom surface (>0)       rn_sbot_max   = ',rn_sbot_max
1936         WRITE(numout,*) '        minimum depth of s-bottom surface (>0)       rn_sbot_min   = ',rn_sbot_min
1937         WRITE(numout,*) '        Critical depth                               rn_hc         = ',rn_hc
1938         WRITE(numout,*) '        maximum cut-off r-value allowed              rn_rmax       = ',rn_rmax
1939         WRITE(numout,*) '     Song and Haidvogel 1994 stretching              ln_s_sh94     = ',ln_s_sh94
1940         WRITE(numout,*) '        Song and Haidvogel 1994 stretching coefficients'
1941         WRITE(numout,*) '        surface control parameter (0<=rn_theta<=20)  rn_theta      = ',rn_theta
1942         WRITE(numout,*) '        bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ',rn_thetb
1943         WRITE(numout,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb
1944         WRITE(numout,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12
1945         WRITE(numout,*) '        switching to sigma (T) or Z (F) at H<Hc      ln_sigcrit    = ',ln_sigcrit
1946         WRITE(numout,*) '        Siddorn and Furner 2012 stretching coefficients'
1947         WRITE(numout,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha
1948         WRITE(numout,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold
1949         WRITE(numout,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs
1950         WRITE(numout,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a
1951         WRITE(numout,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b
1952         WRITE(numout,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b'
1953      ENDIF
1954
1955      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate
1956      hifu(:,:) = rn_sbot_min
1957      hifv(:,:) = rn_sbot_min
1958      hiff(:,:) = rn_sbot_min
1959
1960      !                                        ! set maximum ocean depth
1961      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1962
1963      DO jj = 1, jpj
1964         DO ji = 1, jpi
1965           IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1966         END DO
1967      END DO
1968      !                                        ! =============================
1969      !                                        ! Define the envelop bathymetry   (hbatt)
1970      !                                        ! =============================
1971      ! use r-value to create hybrid coordinates
1972      zenv(:,:) = bathy(:,:)
1973      !
1974      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
1975      DO jj = 1, jpj
1976         DO ji = 1, jpi
1977           IF( bathy(ji,jj) == 0._wp ) THEN
1978             iip1 = MIN( ji+1, jpi )
1979             ijp1 = MIN( jj+1, jpj )
1980             iim1 = MAX( ji-1, 1 )
1981             ijm1 = MAX( jj-1, 1 )
1982             IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) +              &
1983        &         bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
1984               zenv(ji,jj) = rn_sbot_min
1985             ENDIF
1986           ENDIF
1987         END DO
1988      END DO
1989      ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1990      CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1991      !
1992      ! smooth the bathymetry (if required)
1993      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
1994      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
1995      !
1996      jl = 0
1997      zrmax = 1._wp
1998      !   
1999      !     
2000      ! set scaling factor used in reducing vertical gradients
2001      zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
2002      !
2003      ! initialise temporary evelope depth arrays
2004      ztmpi1(:,:) = zenv(:,:)
2005      ztmpi2(:,:) = zenv(:,:)
2006      ztmpj1(:,:) = zenv(:,:)
2007      ztmpj2(:,:) = zenv(:,:)
2008      !
2009      ! initialise temporary r-value arrays
2010      zri(:,:) = 1._wp
2011      zrj(:,:) = 1._wp
2012      !                                                            ! ================ !
2013      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) !  Iterative loop  !
2014         !                                                         ! ================ !
2015         jl = jl + 1
2016         zrmax = 0._wp
2017         ! we set zrmax from previous r-values (zri and zrj) first
2018         ! if set after current r-value calculation (as previously)
2019         ! we could exit DO WHILE prematurely before checking r-value
2020         ! of current zenv
2021         DO jj = 1, nlcj
2022            DO ji = 1, nlci
2023               zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
2024            END DO
2025         END DO
2026         zri(:,:) = 0._wp
2027         zrj(:,:) = 0._wp
2028         DO jj = 1, nlcj
2029            DO ji = 1, nlci
2030               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
2031               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
2032               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
2033                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
2034               END IF
2035               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
2036                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
2037               END IF
2038               IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact
2039               IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact
2040               IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact
2041               IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact
2042            END DO
2043         END DO
2044         IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
2045         !
2046         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax
2047         !
2048         DO jj = 1, nlcj
2049            DO ji = 1, nlci
2050               zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
2051            END DO
2052         END DO
2053         ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
2054         CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
2055         !                                                  ! ================ !
2056      END DO                                                !     End loop     !
2057      !                                                     ! ================ !
2058      DO jj = 1, jpj
2059         DO ji = 1, jpi
2060            zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
2061         END DO
2062      END DO
2063      !
2064      ! Envelope bathymetry saved in hbatt
2065      hbatt(:,:) = zenv(:,:) 
2066      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
2067         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
2068         DO jj = 1, jpj
2069            DO ji = 1, jpi
2070               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
2071               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
2072            END DO
2073         END DO
2074      ENDIF
2075      !
2076      IF(lwp) THEN                             ! Control print
2077         WRITE(numout,*)
2078         WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
2079         WRITE(numout,*)
2080         CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )
2081         IF( nprint == 1 )   THEN       
2082            WRITE(numout,*) ' bathy  MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
2083            WRITE(numout,*) ' hbatt  MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
2084         ENDIF
2085      ENDIF
2086
2087      !                                        ! ==============================
2088      !                                        !   hbatu, hbatv, hbatf fields
2089      !                                        ! ==============================
2090      IF(lwp) THEN
2091         WRITE(numout,*)
2092         WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
2093      ENDIF
2094      hbatu(:,:) = rn_sbot_min
2095      hbatv(:,:) = rn_sbot_min
2096      hbatf(:,:) = rn_sbot_min
2097      DO jj = 1, jpjm1
2098        DO ji = 1, jpim1   ! NO vector opt.
2099           hbatu(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
2100           hbatv(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
2101           hbatf(ji,jj) = 0.25_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
2102              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
2103        END DO
2104      END DO
2105      !
2106      ! Apply lateral boundary condition
2107!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
2108      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk( hbatu, 'U', 1._wp )
2109      DO jj = 1, jpj
2110         DO ji = 1, jpi
2111            IF( hbatu(ji,jj) == 0._wp ) THEN
2112               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
2113               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
2114            ENDIF
2115         END DO
2116      END DO
2117      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk( hbatv, 'V', 1._wp )
2118      DO jj = 1, jpj
2119         DO ji = 1, jpi
2120            IF( hbatv(ji,jj) == 0._wp ) THEN
2121               IF( zhbat(ji,jj) == 0._wp )   hbatv(ji,jj) = rn_sbot_min
2122               IF( zhbat(ji,jj) /= 0._wp )   hbatv(ji,jj) = zhbat(ji,jj)
2123            ENDIF
2124         END DO
2125      END DO
2126      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk( hbatf, 'F', 1._wp )
2127      DO jj = 1, jpj
2128         DO ji = 1, jpi
2129            IF( hbatf(ji,jj) == 0._wp ) THEN
2130               IF( zhbat(ji,jj) == 0._wp )   hbatf(ji,jj) = rn_sbot_min
2131               IF( zhbat(ji,jj) /= 0._wp )   hbatf(ji,jj) = zhbat(ji,jj)
2132            ENDIF
2133         END DO
2134      END DO
2135
2136!!bug:  key_helsinki a verifer
2137      hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
2138      hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
2139      hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
2140      hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
2141
2142      IF( nprint == 1 .AND. lwp )   THEN
2143         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
2144            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
2145         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
2146            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
2147         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
2148            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
2149         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
2150            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
2151      ENDIF
2152!! helsinki
2153
2154      !                                            ! =======================
2155      !                                            !   s-ccordinate fields     (gdep., e3.)
2156      !                                            ! =======================
2157      !
2158      ! non-dimensional "sigma" for model level depth at w- and t-levels
2159
2160
2161!========================================================================
2162! Song and Haidvogel  1994 (ln_s_sh94=T)
2163! Siddorn and Furner 2012 (ln_sf12=T)
2164! or  tanh function       (both false)                   
2165!========================================================================
2166      IF      ( ln_s_sh94 ) THEN
2167                           CALL s_sh94()
2168      ELSE IF ( ln_s_sf12 ) THEN
2169                           CALL s_sf12()
2170      ELSE                 
2171                           CALL s_tanh()
2172      ENDIF
2173
2174      CALL lbc_lnk( e3t_0 , 'T', 1._wp )
2175      CALL lbc_lnk( e3u_0 , 'U', 1._wp )
2176      CALL lbc_lnk( e3v_0 , 'V', 1._wp )
2177      CALL lbc_lnk( e3f_0 , 'F', 1._wp )
2178      CALL lbc_lnk( e3w_0 , 'W', 1._wp )
2179      CALL lbc_lnk( e3uw_0, 'U', 1._wp )
2180      CALL lbc_lnk( e3vw_0, 'V', 1._wp )
2181
2182      fsdepw(:,:,:) = gdepw_0 (:,:,:)
2183      fsde3w(:,:,:) = gdep3w_0(:,:,:)
2184      !
2185      where (e3t_0   (:,:,:) == 0.0)  e3t_0(:,:,:)  = 1.0_wp
2186      where (e3u_0   (:,:,:) == 0.0)  e3u_0(:,:,:)  = 1.0_wp
2187      where (e3v_0   (:,:,:) == 0.0)  e3v_0(:,:,:)  = 1.0_wp
2188      where (e3f_0   (:,:,:) == 0.0)  e3f_0(:,:,:)  = 1.0_wp
2189      where (e3w_0   (:,:,:) == 0.0)  e3w_0(:,:,:)  = 1.0_wp
2190      where (e3uw_0  (:,:,:) == 0.0)  e3uw_0(:,:,:) = 1.0_wp
2191      where (e3vw_0  (:,:,:) == 0.0)  e3vw_0(:,:,:) = 1.0_wp
2192
2193#if defined key_agrif
2194      ! Ensure meaningful vertical scale factors in ghost lines/columns
2195      IF( .NOT. Agrif_Root() ) THEN
2196         
2197         IF((nbondi == -1).OR.(nbondi == 2)) THEN
2198            e3u_0(1,:,:) = e3u_0(2,:,:)
2199         ENDIF
2200         !
2201         IF((nbondi ==  1).OR.(nbondi == 2)) THEN
2202            e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:)
2203         ENDIF
2204         !
2205         IF((nbondj == -1).OR.(nbondj == 2)) THEN
2206            e3v_0(:,1,:) = e3v_0(:,2,:)
2207         ENDIF
2208         !
2209         IF((nbondj ==  1).OR.(nbondj == 2)) THEN
2210            e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:)
2211         ENDIF
2212         !
2213      ENDIF
2214#endif
2215
2216      fsdept(:,:,:) = gdept_0 (:,:,:)
2217      fsdepw(:,:,:) = gdepw_0 (:,:,:)
2218      fsde3w(:,:,:) = gdep3w_0(:,:,:)
2219      fse3t (:,:,:) = e3t_0   (:,:,:)
2220      fse3u (:,:,:) = e3u_0   (:,:,:)
2221      fse3v (:,:,:) = e3v_0   (:,:,:)
2222      fse3f (:,:,:) = e3f_0   (:,:,:)
2223      fse3w (:,:,:) = e3w_0   (:,:,:)
2224      fse3uw(:,:,:) = e3uw_0  (:,:,:)
2225      fse3vw(:,:,:) = e3vw_0  (:,:,:)
2226!!
2227      ! HYBRID :
2228      DO jj = 1, jpj
2229         DO ji = 1, jpi
2230            DO jk = 1, jpkm1
2231               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
2232            END DO
2233            IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0
2234         END DO
2235      END DO
2236      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
2237         &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
2238
2239      IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain
2240         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) )
2241         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
2242            &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w '  , MINVAL( gdep3w_0(:,:,:) )
2243         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0   (:,:,:) ),   &
2244            &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0   (:,:,:) ),   &
2245            &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0  (:,:,:) ),   &
2246            &                          ' w ', MINVAL( e3w_0  (:,:,:) )
2247
2248         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
2249            &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w '  , MAXVAL( gdep3w_0(:,:,:) )
2250         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0   (:,:,:) ),   &
2251            &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0   (:,:,:) ),   &
2252            &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0  (:,:,:) ),   &
2253            &                          ' w ', MAXVAL( e3w_0  (:,:,:) )
2254      ENDIF
2255      !  END DO
2256      IF(lwp) THEN                                  ! selected vertical profiles
2257         WRITE(numout,*)
2258         WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
2259         WRITE(numout,*) ' ~~~~~~  --------------------'
2260         WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2261         WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     &
2262            &                                 e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
2263         DO jj = mj0(20), mj1(20)
2264            DO ji = mi0(20), mi1(20)
2265               WRITE(numout,*)
2266               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
2267               WRITE(numout,*) ' ~~~~~~  --------------------'
2268               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2269               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
2270                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2271            END DO
2272         END DO
2273         DO jj = mj0(74), mj1(74)
2274            DO ji = mi0(100), mi1(100)
2275               WRITE(numout,*)
2276               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
2277               WRITE(numout,*) ' ~~~~~~  --------------------'
2278               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2279               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
2280                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2281            END DO
2282         END DO
2283      ENDIF
2284
2285!================================================================================
2286! check the coordinate makes sense
2287!================================================================================
2288      DO ji = 1, jpi
2289         DO jj = 1, jpj
2290
2291            IF( hbatt(ji,jj) > 0._wp) THEN
2292               DO jk = 1, mbathy(ji,jj)
2293                 ! check coordinate is monotonically increasing
2294                 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
2295                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
2296                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
2297                    WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
2298                    WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
2299                    CALL ctl_stop( ctmp1 )
2300                 ENDIF
2301                 ! and check it has never gone negative
2302                 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
2303                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
2304                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk
2305                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2306                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2307                    CALL ctl_stop( ctmp1 )
2308                 ENDIF
2309                 ! and check it never exceeds the total depth
2310                 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
2311                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
2312                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
2313                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2314                    CALL ctl_stop( ctmp1 )
2315                 ENDIF
2316               END DO
2317
2318               DO jk = 1, mbathy(ji,jj)-1
2319                 ! and check it never exceeds the total depth
2320                IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
2321                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
2322                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
2323                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2324                    CALL ctl_stop( ctmp1 )
2325                 ENDIF
2326               END DO
2327
2328            ENDIF
2329
2330         END DO
2331      END DO
2332      !
2333      CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
2334      !
2335      IF( nn_timing == 1 )  CALL timing_stop('zgr_sco')
2336      !
2337   END SUBROUTINE zgr_sco
2338
2339!!======================================================================
2340   SUBROUTINE s_sh94()
2341
2342      !!----------------------------------------------------------------------
2343      !!                  ***  ROUTINE s_sh94  ***
2344      !!                     
2345      !! ** Purpose :   stretch the s-coordinate system
2346      !!
2347      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994
2348      !!                mixed S/sigma coordinate
2349      !!
2350      !! Reference : Song and Haidvogel 1994.
2351      !!----------------------------------------------------------------------
2352      !
2353      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2354      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
2355      !
2356      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2357      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
2358
2359      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2360      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2361
2362      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
2363      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
2364      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
2365      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
2366
2367      DO ji = 1, jpi
2368         DO jj = 1, jpj
2369
2370            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
2371               DO jk = 1, jpk
2372                  z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
2373                  z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
2374               END DO
2375            ELSE ! shallow water, uniform sigma
2376               DO jk = 1, jpk
2377                  z_gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
2378                  z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
2379                  END DO
2380            ENDIF
2381            !
2382            DO jk = 1, jpkm1
2383               z_esigt3(ji,jj,jk  ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2384               z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2385            END DO
2386            z_esigw3(ji,jj,1  ) = 2._wp * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) )
2387            z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
2388            !
2389            ! Coefficients for vertical depth as the sum of e3w scale factors
2390            z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
2391            DO jk = 2, jpk
2392               z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2393            END DO
2394            !
2395            DO jk = 1, jpk
2396               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2397               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2398               gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
2399               gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
2400               gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
2401            END DO
2402           !
2403         END DO   ! for all jj's
2404      END DO    ! for all ji's
2405
2406      DO ji = 1, jpim1
2407         DO jj = 1, jpjm1
2408            DO jk = 1, jpk
2409               z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
2410                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2411               z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
2412                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2413               z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     &
2414                  &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   &
2415                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2416               z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
2417                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2418               z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
2419                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2420               !
2421               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2422               e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2423               e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2424               e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2425               !
2426               e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2427               e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2428               e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2429            END DO
2430        END DO
2431      END DO
2432
2433      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2434      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2435
2436   END SUBROUTINE s_sh94
2437
2438   SUBROUTINE s_sf12
2439
2440      !!----------------------------------------------------------------------
2441      !!                  ***  ROUTINE s_sf12 ***
2442      !!                     
2443      !! ** Purpose :   stretch the s-coordinate system
2444      !!
2445      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
2446      !!                mixed S/sigma/Z coordinate
2447      !!
2448      !!                This method allows the maintenance of fixed surface and or
2449      !!                bottom cell resolutions (cf. geopotential coordinates)
2450      !!                within an analytically derived stretched S-coordinate framework.
2451      !!
2452      !!
2453      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
2454      !!----------------------------------------------------------------------
2455      !
2456      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2457      REAL(wp) ::   zsmth               ! smoothing around critical depth
2458      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
2459      !
2460      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2461      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
2462
2463      !
2464      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2465      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2466
2467      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
2468      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
2469      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
2470      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
2471
2472      DO ji = 1, jpi
2473         DO jj = 1, jpj
2474
2475          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
2476             
2477              zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
2478                                                     ! could be changed by users but care must be taken to do so carefully
2479              zzb = 1.0_wp-(zzb/hbatt(ji,jj))
2480           
2481              zzs = rn_zs / hbatt(ji,jj) 
2482             
2483              IF (rn_efold /= 0.0_wp) THEN
2484                zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
2485              ELSE
2486                zsmth = 1.0_wp 
2487              ENDIF
2488               
2489              DO jk = 1, jpk
2490                z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
2491                z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
2492              ENDDO
2493              z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth  )
2494              z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth  )
2495 
2496          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
2497
2498            DO jk = 1, jpk
2499              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
2500              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
2501            END DO
2502
2503          ELSE  ! shallow water, z coordinates
2504
2505            DO jk = 1, jpk
2506              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2507              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2508            END DO
2509
2510          ENDIF
2511
2512          DO jk = 1, jpkm1
2513             z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2514             z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2515          END DO
2516          z_esigw3(ji,jj,1  ) = 2.0_wp * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ))
2517          z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
2518
2519          ! Coefficients for vertical depth as the sum of e3w scale factors
2520          z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
2521          DO jk = 2, jpk
2522             z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2523          END DO
2524
2525          DO jk = 1, jpk
2526             gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
2527             gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
2528             gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
2529          END DO
2530
2531        ENDDO   ! for all jj's
2532      ENDDO    ! for all ji's
2533
2534      DO ji=1,jpi-1
2535        DO jj=1,jpj-1
2536
2537          DO jk = 1, jpk
2538                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / &
2539                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2540                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / &
2541                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2542                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  &
2543                                      hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / &
2544                                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2545                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / &
2546                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2547                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / &
2548                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2549
2550             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
2551             e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
2552             e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
2553             e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
2554             !
2555             e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
2556             e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
2557             e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
2558          END DO
2559
2560        ENDDO
2561      ENDDO
2562      !
2563      CALL lbc_lnk(e3t_0 ,'T',1.) ; CALL lbc_lnk(e3u_0 ,'T',1.)
2564      CALL lbc_lnk(e3v_0 ,'T',1.) ; CALL lbc_lnk(e3f_0 ,'T',1.)
2565      CALL lbc_lnk(e3w_0 ,'T',1.)
2566      CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.)
2567      !
2568      !                                               ! =============
2569
2570      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2571      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2572
2573   END SUBROUTINE s_sf12
2574
2575   SUBROUTINE s_tanh()
2576
2577      !!----------------------------------------------------------------------
2578      !!                  ***  ROUTINE s_tanh***
2579      !!                     
2580      !! ** Purpose :   stretch the s-coordinate system
2581      !!
2582      !! ** Method  :   s-coordinate stretch
2583      !!
2584      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
2585      !!----------------------------------------------------------------------
2586
2587      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2588      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
2589
2590      REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
2591      REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw
2592
2593      CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      )
2594      CALL wrk_alloc( jpk, z_esigt, z_esigw                                               )
2595
2596      z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp   ;   z_gsi3w  = 0._wp
2597      z_esigt  = 0._wp   ;   z_esigw  = 0._wp 
2598
2599      DO jk = 1, jpk
2600        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
2601        z_gsigt(jk) = -fssig( REAL(jk,wp)        )
2602      END DO
2603      IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk)
2604      !
2605      ! Coefficients for vertical scale factors at w-, t- levels
2606!!gm bug :  define it from analytical function, not like juste bellow....
2607!!gm        or betteroffer the 2 possibilities....
2608      DO jk = 1, jpkm1
2609         z_esigt(jk  ) = z_gsigw(jk+1) - z_gsigw(jk)
2610         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
2611      END DO
2612      z_esigw( 1 ) = 2._wp * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
2613      z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
2614      !
2615      ! Coefficients for vertical depth as the sum of e3w scale factors
2616      z_gsi3w(1) = 0.5_wp * z_esigw(1)
2617      DO jk = 2, jpk
2618         z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
2619      END DO
2620!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
2621      DO jk = 1, jpk
2622         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2623         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2624         gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
2625         gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
2626         gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )
2627      END DO
2628!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
2629      DO jj = 1, jpj
2630         DO ji = 1, jpi
2631            DO jk = 1, jpk
2632              e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2633              e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2634              e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2635              e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
2636              !
2637              e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2638              e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2639              e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2640            END DO
2641         END DO
2642      END DO
2643
2644      CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      )
2645      CALL wrk_dealloc( jpk, z_esigt, z_esigw                                               )
2646
2647   END SUBROUTINE s_tanh
2648
2649   FUNCTION fssig( pk ) RESULT( pf )
2650      !!----------------------------------------------------------------------
2651      !!                 ***  ROUTINE fssig ***
2652      !!       
2653      !! ** Purpose :   provide the analytical function in s-coordinate
2654      !!         
2655      !! ** Method  :   the function provide the non-dimensional position of
2656      !!                T and W (i.e. between 0 and 1)
2657      !!                T-points at integer values (between 1 and jpk)
2658      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2659      !!----------------------------------------------------------------------
2660      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
2661      REAL(wp)             ::   pf   ! sigma value
2662      !!----------------------------------------------------------------------
2663      !
2664      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   &
2665         &     - TANH( rn_thetb * rn_theta                                )  )   &
2666         & * (   COSH( rn_theta                           )                      &
2667         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              &
2668         & / ( 2._wp * SINH( rn_theta ) )
2669      !
2670   END FUNCTION fssig
2671
2672
2673   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
2674      !!----------------------------------------------------------------------
2675      !!                 ***  ROUTINE fssig1 ***
2676      !!
2677      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
2678      !!
2679      !! ** Method  :   the function provides the non-dimensional position of
2680      !!                T and W (i.e. between 0 and 1)
2681      !!                T-points at integer values (between 1 and jpk)
2682      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2683      !!----------------------------------------------------------------------
2684      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
2685      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
2686      REAL(wp)             ::   pf1   ! sigma value
2687      !!----------------------------------------------------------------------
2688      !
2689      IF ( rn_theta == 0 ) then      ! uniform sigma
2690         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
2691      ELSE                        ! stretched sigma
2692         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
2693            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
2694            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
2695      ENDIF
2696      !
2697   END FUNCTION fssig1
2698
2699
2700   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
2701      !!----------------------------------------------------------------------
2702      !!                 ***  ROUTINE fgamma  ***
2703      !!
2704      !! ** Purpose :   provide analytical function for the s-coordinate
2705      !!
2706      !! ** Method  :   the function provides the non-dimensional position of
2707      !!                T and W (i.e. between 0 and 1)
2708      !!                T-points at integer values (between 1 and jpk)
2709      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2710      !!
2711      !!                This method allows the maintenance of fixed surface and or
2712      !!                bottom cell resolutions (cf. geopotential coordinates)
2713      !!                within an analytically derived stretched S-coordinate framework.
2714      !!
2715      !! Reference  :   Siddorn and Furner, in prep
2716      !!----------------------------------------------------------------------
2717      REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate
2718      REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate
2719      REAL(wp), INTENT(in   ) ::   pzb           ! Bottom box depth
2720      REAL(wp), INTENT(in   ) ::   pzs           ! surface box depth
2721      REAL(wp), INTENT(in   ) ::   psmth       ! Smoothing parameter
2722      REAL(wp)                ::   za1,za2,za3    ! local variables
2723      REAL(wp)                ::   zn1,zn2        ! local variables
2724      REAL(wp)                ::   za,zb,zx       ! local variables
2725      integer                 ::   jk
2726      !!----------------------------------------------------------------------
2727      !
2728
2729      zn1  =  1./(jpk-1.)
2730      zn2  =  1. -  zn1
2731
2732      za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 
2733      za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp)
2734      za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
2735     
2736      za = pzb - za3*(pzs-za1)-za2
2737      za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) )
2738      zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
2739      zx = 1.0_wp-za/2.0_wp-zb
2740 
2741      DO jk = 1, jpk
2742        p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp +  &
2743                    & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &
2744                    &      (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
2745        p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth)
2746      ENDDO 
2747
2748      !
2749   END FUNCTION fgamma
2750
2751   !!======================================================================
2752END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.