source: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 4724

Last change on this file since 4724 was 4724, checked in by mathiot, 6 years ago

ISF branch: add comments, fix mpp and restar issues, add test to stop if incompatible options and fix mask issue in sbcice and sbcblk.

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