source: utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/domzgr.F90 @ 13056

Last change on this file since 13056 was 13056, checked in by rblod, 4 months ago

ticket #2129 : cleaning domcfg

File size: 108.8 KB
Line 
1MODULE domzgr
2   !!==============================================================================
3   !!                       ***  MODULE domzgr   ***
4   !! Ocean domain : definition of the vertical coordinate system
5   !!==============================================================================
6   !! History :  OPA  ! 1995-12  (G. Madec)  Original code : s vertical coordinate
7   !!                 ! 1997-07  (G. Madec)  lbc_lnk call
8   !!                 ! 1997-04  (J.-O. Beismann)
9   !!            8.5  ! 2002-09  (A. Bozec, G. Madec)  F90: Free form and module
10   !!             -   ! 2002-09  (A. de Miranda)  rigid-lid + islands
11   !!  NEMO      1.0  ! 2003-08  (G. Madec)  F90: Free form and module
12   !!             -   ! 2005-10  (A. Beckmann)  modifications for hybrid s-ccordinates & new stretching function
13   !!            2.0  ! 2006-04  (R. Benshila, G. Madec)  add zgr_zco
14   !!            3.0  ! 2008-06  (G. Madec)  insertion of domzgr_zps.h90 & conding style
15   !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option
16   !!            3.3  ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level
17   !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn and Furner stretching function
18   !!            3.4  ! 2012-12  (R. Bourdalle-Badie and G. Reffray)  modify C1D case 
19   !!            3.6  ! 2014-11  (P. Mathiot and C. Harris) add ice shelf capabilitye
20   !!            3.?  ! 2015-11  (H. Liu) Modifications for Wetting/Drying
21   !!----------------------------------------------------------------------
22
23   !!----------------------------------------------------------------------
24   !!   dom_zgr          : defined the ocean vertical coordinate system
25   !!       zgr_bat      : bathymetry fields (levels and meters)
26   !!       zgr_bat_zoom : modify the bathymetry field if zoom domain
27   !!       zgr_bat_ctl  : check the bathymetry files
28   !!       zgr_bot_level: deepest ocean level for t-, u, and v-points
29   !!       zgr_z        : reference z-coordinate
30   !!       zgr_zco      : z-coordinate
31   !!       zgr_zps      : z-coordinate with partial steps
32   !!       zgr_sco      : s-coordinate
33   !!       fssig        : tanh stretch function
34   !!       fssig1       : Song and Haidvogel 1994 stretch function
35   !!       fgamma       : Siddorn and Furner 2012 stretching function
36   !!---------------------------------------------------------------------
37   USE dom_oce           ! ocean domain
38   USE depth_e3       ! depth <=> e3
39   !
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 lib_fortran
45   USE dombat
46   USE domisf
47   USE agrif_domzgr
48
49   IMPLICIT NONE
50   PRIVATE
51
52   PUBLIC   dom_zgr        ! called by dom_init.F90
53   !
54   !                              !!* Namelist namzgr_sco *
55   LOGICAL  ::   ln_s_sh94         ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T)
56   LOGICAL  ::   ln_s_sf12         ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T)
57   !
58   REAL(wp) ::   rn_sbot_min       ! minimum depth of s-bottom surface (>0) (m)
59   REAL(wp) ::   rn_sbot_max       ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
60   REAL(wp) ::   rn_rmax           ! maximum cut-off r-value allowed (0<rn_rmax<1)
61   REAL(wp) ::   rn_hc             ! Critical depth for transition from sigma to stretched coordinates
62
63   ! Song and Haidvogel 1994 stretching parameters
64   REAL(wp) ::   rn_theta          ! surface control parameter (0<=rn_theta<=20)
65   REAL(wp) ::   rn_thetb          ! bottom control parameter  (0<=rn_thetb<= 1)
66   REAL(wp) ::   rn_bb             ! stretching parameter
67   !                                        ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
68   ! Siddorn and Furner stretching parameters
69   LOGICAL  ::   ln_sigcrit        ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
70   REAL(wp) ::   rn_alpha          ! control parameter ( > 1 stretch towards surface, < 1 towards seabed)
71   REAL(wp) ::   rn_efold          !  efold length scale for transition to stretched coord
72   REAL(wp) ::   rn_zs             !  depth of surface grid box
73                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
74   REAL(wp) ::   rn_zb_a           !  bathymetry scaling factor for calculating Zb
75   REAL(wp) ::   rn_zb_b           !  offset for calculating Zb
76
77  !! * Substitutions
78   !!----------------------------------------------------------------------
79   !!                   ***  vectopt_loop_substitute  ***
80   !!----------------------------------------------------------------------
81   !! ** purpose :   substitute the inner loop start/end indices with CPP macro
82   !!                allow unrolling of do-loop (useful with vector processors)
83   !!----------------------------------------------------------------------
84   !!----------------------------------------------------------------------
85   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
86   !! $Id: vectopt_loop_substitute.h90 4990 2014-12-15 16:42:49Z timgraham $
87   !! Software governed by the CeCILL licence (./LICENSE)
88   !!----------------------------------------------------------------------
89   !!----------------------------------------------------------------------
90   !! NEMO/OPA 3.3.1 , NEMO Consortium (2011)
91   !! $Id: domzgr.F90 6827 2016-08-01 13:37:15Z flavoni $
92   !! Software governed by the CeCILL licence     (./LICENSE)
93   !!----------------------------------------------------------------------
94CONTAINS       
95
96   SUBROUTINE dom_zgr( k_top, k_bot )
97      !!----------------------------------------------------------------------
98      !!                ***  ROUTINE dom_zgr  ***
99      !!                   
100      !! ** Purpose :   set the depth of model levels and the resulting
101      !!              vertical scale factors.
102      !!
103      !! ** Method  : - reference 1D vertical coordinate (gdep._1d, e3._1d)
104      !!              - read/set ocean depth and ocean levels (bathy, mbathy)
105      !!              - vertical coordinate (gdep., e3.) depending on the
106      !!                coordinate chosen :
107      !!                   ln_zco=T   z-coordinate   
108      !!                   ln_zps=T   z-coordinate with partial steps
109      !!                   ln_zco=T   s-coordinate
110      !!
111      !! ** Action  :   define gdep., e3., mbathy and bathy
112      !!----------------------------------------------------------------------
113      INTEGER, DIMENSION(:,:), INTENT(out) ::   k_top, k_bot   ! ocean first and last level indices
114      !
115      INTEGER ::   ioptio, ibat   ! local integer
116      INTEGER ::   ios
117      !
118      INTEGER :: jk
119      REAL(wp) ::   zrefdep             ! depth of the reference level (~10m)
120
121
122      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav
123      !!----------------------------------------------------------------------
124      !
125      !
126      REWIND( numnam_ref )              ! Namelist namzgr in reference namelist : Vertical coordinate
127      READ  ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901 )
128901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist', lwp )
129
130      REWIND( numnam_cfg )              ! Namelist namzgr in configuration namelist : Vertical coordinate
131      READ  ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 )
132902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp )
133      IF(lwm) WRITE ( numond, namzgr )
134
135      IF(ln_read_cfg) THEN
136         IF(lwp) WRITE(numout,*)
137         IF(lwp) WRITE(numout,*) '   ==>>>   Read vertical mesh in ', TRIM( cn_domcfg ), ' file'
138         !
139         CALL zgr_read   ( ln_zco  , ln_zps  , ln_sco, ln_isfcav,   &
140            &              gdept_1d, gdepw_1d, e3t_1d, e3w_1d   ,   &    ! 1D gridpoints depth
141            &              gdept_0 , gdepw_0                    ,   &    ! gridpoints depth
142            &              e3t_0   , e3u_0   , e3v_0 , e3f_0    ,   &    ! vertical scale factors
143            &              e3w_0   , e3uw_0  , e3vw_0           ,   &    ! vertical scale factors
144            &              k_top   , k_bot            )                  ! 1st & last ocean level
145            !
146!!gm to be remove when removing the OLD definition of e3 scale factors so that gde3w disappears
147!      ! Compute gde3w_0 (vertical sum of e3w)
148!      gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
149!      DO jk = 2, jpk
150!         gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)
151!      END DO
152      !
153
154      !                                ! top/bottom ocean level indices for t-, u- and v-points (f-point also for top)
155      CALL zgr_top_bot( k_top, k_bot )      ! with a minimum value set to 1
156
157      !                                ! deepest/shallowest W level Above/Below ~10m
158!!gm BUG in s-coordinate this does not work!
159      zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d )                   ! ref. depth with tolerance (10% of minimum layer thickness)
160      nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
161      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m
162!!gm end bug
163
164      ENDIF       
165
166      IF(lwp) THEN                     ! Control print
167         WRITE(numout,*)
168         WRITE(numout,*) 'dom_zgr : vertical coordinate'
169         WRITE(numout,*) '~~~~~~~'
170         WRITE(numout,*) '   Namelist namzgr : set vertical coordinate'
171         WRITE(numout,*) '      z-coordinate - full steps      ln_zco    = ', ln_zco
172         WRITE(numout,*) '      z-coordinate - partial steps   ln_zps    = ', ln_zps
173         WRITE(numout,*) '      s- or hybrid z-s-coordinate    ln_sco    = ', ln_sco
174         WRITE(numout,*) '      ice shelf cavities             ln_isfcav = ', ln_isfcav
175      ENDIF
176
177      ioptio = 0                       ! Check Vertical coordinate options
178      IF( ln_zco      )   ioptio = ioptio + 1
179      IF( ln_zps      )   ioptio = ioptio + 1
180      IF( ln_sco      )   ioptio = ioptio + 1
181      IF( ioptio /= 1 )   CALL ctl_stop( ' none or several vertical coordinate options used' )
182      !
183      IF ( ln_isfcav ) CALL zgr_isf_nam
184      ioptio = 0
185      IF ( ln_zco .AND. ln_isfcav ) ioptio = ioptio + 1
186      IF ( ln_sco .AND. ln_isfcav ) ioptio = ioptio + 1
187      IF( ioptio > 0 )   CALL ctl_stop( ' Cavity not tested/compatible with full step (zco) and sigma (ln_sco) ' )
188
189      IF(.NOT.ln_read_cfg) THEN
190      !
191      ! Build the vertical coordinate system
192      ! ------------------------------------
193                          CALL zgr_z            ! Reference z-coordinate system (always called)
194                          CALL zgr_bat          ! Bathymetry fields (levels and meters)
195      IF( ln_zco      )   CALL zgr_zco          ! z-coordinate
196      IF( ln_zps      )   CALL zgr_zps          ! Partial step z-coordinate
197      IF( ln_sco      )   CALL zgr_sco          ! s-coordinate or hybrid z-s coordinate
198                          CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isolated ocean points
199      !
200      ! final adjustment of mbathy & check
201      ! -----------------------------------
202                          CALL zgr_bot_level    ! deepest ocean level for t-, u- and v-points
203                          k_bot = mbkt
204                          CALL zgr_top_level    ! shallowest ocean level for T-, U-, V- points
205                          k_top = mikt
206
207      ENDIF
208      !
209      IF( lwp )   THEN
210         WRITE(numout,*) ' MIN val k_top   ', MINVAL(   k_top(:,:) ), ' MAX ', MAXVAL( k_top(:,:) )
211         WRITE(numout,*) ' MIN val k_bot   ', MINVAL(   k_bot(:,:) ), ' MAX ', MAXVAL( k_bot(:,:) )
212         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
213            &                          ' w ', MINVAL( gdepw_0(:,:,:) )
214         WRITE(numout,*) ' MIN val e3    t ', MINVAL(   e3t_0(:,:,:) ), ' f ', MINVAL(   e3f_0(:,:,:) ),  &
215            &                          ' u ', MINVAL(   e3u_0(:,:,:) ), ' u ', MINVAL(   e3v_0(:,:,:) ),  &
216            &                          ' uw', MINVAL(  e3uw_0(:,:,:) ), ' vw', MINVAL(  e3vw_0(:,:,:)),   &
217            &                          ' w ', MINVAL(   e3w_0(:,:,:) )
218
219         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
220            &                          ' w ', MAXVAL( gdepw_0(:,:,:) )
221         WRITE(numout,*) ' MAX val e3    t ', MAXVAL(   e3t_0(:,:,:) ), ' f ', MAXVAL(   e3f_0(:,:,:) ),  &
222            &                          ' u ', MAXVAL(   e3u_0(:,:,:) ), ' u ', MAXVAL(   e3v_0(:,:,:) ),  &
223            &                          ' uw', MAXVAL(  e3uw_0(:,:,:) ), ' vw', MAXVAL(  e3vw_0(:,:,:) ),  &
224            &                          ' w ', MAXVAL(   e3w_0(:,:,:) )
225      ENDIF
226
227      !
228   END SUBROUTINE dom_zgr
229
230   SUBROUTINE zgr_read( ld_zco  , ld_zps  , ld_sco  , ld_isfcav,   &   ! type of vertical coordinate
231      &                 pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d  ,   &   ! 1D reference vertical coordinate
232      &                 pdept , pdepw ,                            &   ! 3D t & w-points depth
233      &                 pe3t  , pe3u  , pe3v   , pe3f ,            &   ! vertical scale factors
234      &                 pe3w  , pe3uw , pe3vw         ,            &   !     -      -      -
235      &                 k_top  , k_bot    )                            ! top & bottom ocean level
236      !!---------------------------------------------------------------------
237      !!              ***  ROUTINE zgr_read  ***
238      !!
239      !! ** Purpose :   Read the vertical information in the domain configuration file
240      !!
241      !!----------------------------------------------------------------------
242      LOGICAL                   , INTENT(out) ::   ld_zco, ld_zps, ld_sco      ! vertical coordinate flags
243      LOGICAL                   , INTENT(out) ::   ld_isfcav                   ! under iceshelf cavity flag
244      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth       [m]
245      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d           ! 1D vertical scale factors [m]
246      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pdept, pdepw                ! grid-point depth          [m]
247      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors    [m]
248      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3w , pe3uw, pe3vw         !    -       -      -
249      INTEGER , DIMENSION(:,:)  , INTENT(out) ::   k_top , k_bot               ! first & last ocean level
250      !
251      INTEGER  ::   jk     ! dummy loop index
252      INTEGER  ::   inum   ! local logical unit
253      REAL(WP) ::   z_zco, z_zps, z_sco, z_cav
254      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D workspace
255      !!----------------------------------------------------------------------
256      !
257      IF(lwp) THEN
258         WRITE(numout,*)
259         WRITE(numout,*) '   zgr_read : read the vertical coordinates in ', TRIM( cn_domcfg ), ' file'
260         WRITE(numout,*) '   ~~~~~~~~'
261      ENDIF
262      !
263      CALL iom_open( cn_domcfg, inum )
264      !
265      !                          !* type of vertical coordinate
266      CALL iom_get( inum, 'ln_zco'   , z_zco )
267      CALL iom_get( inum, 'ln_zps'   , z_zps )
268      CALL iom_get( inum, 'ln_sco'   , z_sco )
269      IF( z_zco == 0._wp ) THEN   ;   ld_zco = .false.   ;   ELSE   ;   ld_zco = .true.   ;   ENDIF
270      IF( z_zps == 0._wp ) THEN   ;   ld_zps = .false.   ;   ELSE   ;   ld_zps = .true.   ;   ENDIF
271      IF( z_sco == 0._wp ) THEN   ;   ld_sco = .false.   ;   ELSE   ;   ld_sco = .true.   ;   ENDIF
272      !
273      !                          !* ocean cavities under iceshelves
274      CALL iom_get( inum, 'ln_isfcav', z_cav )
275      IF( z_cav == 0._wp ) THEN   ;   ld_isfcav = .false.   ;   ELSE   ;   ld_isfcav = .true.   ;   ENDIF
276      !
277      !                          !* vertical scale factors
278      CALL iom_get( inum, jpdom_unknown, 'e3t_1d'  , pe3t_1d  )                     ! 1D reference coordinate
279      CALL iom_get( inum, jpdom_unknown, 'e3w_1d'  , pe3w_1d  )
280      !
281      CALL iom_get( inum, jpdom_data, 'e3t_0'  , pe3t  , lrowattr=ln_use_jattr )    ! 3D coordinate
282      CALL iom_get( inum, jpdom_data, 'e3u_0'  , pe3u  , lrowattr=ln_use_jattr )
283      CALL iom_get( inum, jpdom_data, 'e3v_0'  , pe3v  , lrowattr=ln_use_jattr )
284      CALL iom_get( inum, jpdom_data, 'e3f_0'  , pe3f  , lrowattr=ln_use_jattr )
285      CALL iom_get( inum, jpdom_data, 'e3w_0'  , pe3w  , lrowattr=ln_use_jattr )
286      CALL iom_get( inum, jpdom_data, 'e3uw_0' , pe3uw , lrowattr=ln_use_jattr )
287      CALL iom_get( inum, jpdom_data, 'e3vw_0' , pe3vw , lrowattr=ln_use_jattr )
288      !
289      !                          !* depths
290      !                                   !- old depth definition (obsolescent feature)
291      IF(  iom_varid( inum, 'gdept_1d', ldstop = .FALSE. ) > 0  .AND.  &
292         & iom_varid( inum, 'gdepw_1d', ldstop = .FALSE. ) > 0  .AND.  &
293         & iom_varid( inum, 'gdept_0' , ldstop = .FALSE. ) > 0  .AND.  &
294         & iom_varid( inum, 'gdepw_0' , ldstop = .FALSE. ) > 0    ) THEN
295         CALL ctl_warn( 'zgr_read : old definition of depths and scale factors used ', & 
296            &           '           depths at t- and w-points read in the domain configuration file')
297         CALL iom_get( inum, jpdom_unknown, 'gdept_1d', pdept_1d )   
298         CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', pdepw_1d )
299         CALL iom_get( inum, jpdom_data   , 'gdept_0' , pdept , lrowattr=ln_use_jattr )
300         CALL iom_get( inum, jpdom_data   , 'gdepw_0' , pdepw , lrowattr=ln_use_jattr )
301         !
302      ELSE                                !- depths computed from e3. scale factors
303         CALL e3_to_depth( pe3t_1d, pe3w_1d, pdept_1d, pdepw_1d )    ! 1D reference depth
304         CALL e3_to_depth( pe3t   , pe3w   , pdept   , pdepw    )    ! 3D depths
305         IF(lwp) THEN
306            WRITE(numout,*)
307            WRITE(numout,*) '              Reference 1D z-coordinate depth and scale factors:'
308            WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
309            WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
310         ENDIF
311      ENDIF
312      !
313      !                          !* ocean top and bottom level
314      CALL iom_get( inum, jpdom_data, 'top_level'    , z2d  , lrowattr=ln_use_jattr )   ! 1st wet T-points (ISF)
315      k_top(:,:) = NINT( z2d(:,:) )
316      CALL iom_get( inum, jpdom_data, 'bottom_level' , z2d  , lrowattr=ln_use_jattr )   ! last wet T-points
317      k_bot(:,:) = NINT( z2d(:,:) )
318      !
319      ! reference depth for negative bathy (wetting and drying only)
320      ! IF( ll_wd )  CALL iom_get( inum,  'rn_wd_ref_depth' , ssh_ref   )
321      !
322      CALL iom_close( inum )
323      !
324   END SUBROUTINE zgr_read
325
326   SUBROUTINE zgr_top_bot( k_top, k_bot )
327      !!----------------------------------------------------------------------
328      !!                    ***  ROUTINE zgr_top_bot  ***
329      !!
330      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
331      !!
332      !! ** Method  :   computes from k_top and k_bot with a minimum value of 1 over land
333      !!
334      !! ** Action  :   mikt, miku, mikv :   vertical indices of the shallowest
335      !!                                     ocean level at t-, u- & v-points
336      !!                                     (min value = 1)
337      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest
338      !!                                     ocean level at t-, u- & v-points
339      !!                                     (min value = 1 over land)
340      !!----------------------------------------------------------------------
341      INTEGER , DIMENSION(:,:), INTENT(in) ::   k_top, k_bot   ! top & bottom ocean level indices
342      !
343      INTEGER ::   ji, jj   ! dummy loop indices
344      REAL(wp), DIMENSION(jpi,jpj) ::   zk   ! workspace
345      !!----------------------------------------------------------------------
346      !
347      IF(lwp) WRITE(numout,*)
348      IF(lwp) WRITE(numout,*) '    zgr_top_bot : ocean top and bottom k-index of T-, U-, V- and W-levels '
349      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
350      !
351      mikt(:,:) = MAX( k_top(:,:) , 1 )    ! top    ocean k-index of T-level (=1 over land)
352      !
353      mbkt(:,:) = MAX( k_bot(:,:) , 1 )    ! bottom ocean k-index of T-level (=1 over land)
354 
355      !                                    ! N.B.  top     k-index of W-level = mikt
356      !                                    !       bottom  k-index of W-level = mbkt+1
357      DO jj = 1, jpjm1
358         DO ji = 1, jpim1
359            miku(ji,jj) = MAX(  mikt(ji+1,jj  ) , mikt(ji,jj)  )
360            mikv(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj)  )
361            mikf(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj), mikt(ji+1,jj  ), mikt(ji+1,jj+1)  )
362            !
363            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
364            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
365         END DO
366      END DO
367      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
368      zk(:,:) = REAL( miku(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'U', 1. )   ;   miku(:,:) = MAX( NINT( zk(:,:) ), 1 )
369      zk(:,:) = REAL( mikv(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'V', 1. )   ;   mikv(:,:) = MAX( NINT( zk(:,:) ), 1 )
370      zk(:,:) = REAL( mikf(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'F', 1. )   ;   mikf(:,:) = MAX( NINT( zk(:,:) ), 1 )
371      !
372      zk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'U', 1. )   ;   mbku(:,:) = MAX( NINT( zk(:,:) ), 1 )
373      zk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'V', 1. )   ;   mbkv(:,:) = MAX( NINT( zk(:,:) ), 1 )
374      !
375   END SUBROUTINE zgr_top_bot
376
377   SUBROUTINE zgr_z
378      !!----------------------------------------------------------------------
379      !!                   ***  ROUTINE zgr_z  ***
380      !!                   
381      !! ** Purpose :   set the depth of model levels and the resulting
382      !!      vertical scale factors.
383      !!
384      !! ** Method  :   z-coordinate system (use in all type of coordinate)
385      !!        The depth of model levels is defined from an analytical
386      !!      function the derivative of which gives the scale factors.
387      !!        both depth and scale factors only depend on k (1d arrays).
388      !!              w-level: gdepw_1d  = gdep(k)
389      !!                       e3w_1d(k) = dk(gdep)(k)     = e3(k)
390      !!              t-level: gdept_1d  = gdep(k+0.5)
391      !!                       e3t_1d(k) = dk(gdep)(k+0.5) = e3(k+0.5)
392      !!
393      !! ** Action  : - gdept_1d, gdepw_1d : depth of T- and W-point (m)
394      !!              - e3t_1d  , e3w_1d   : scale factors at T- and W-levels (m)
395      !!
396      !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
397      !!----------------------------------------------------------------------
398      INTEGER  ::   jk                     ! dummy loop indices
399      REAL(wp) ::   zt, zw                 ! temporary scalars
400      REAL(wp) ::   zsur, za0, za1, zkth   ! Values set from parameters in
401      REAL(wp) ::   zacr, zdzmin, zhmax    ! par_CONFIG_Rxx.h90
402      REAL(wp) ::   zrefdep                ! depth of the reference level (~10m)
403      REAL(wp) ::   za2, zkth2, zacr2      ! Values for optional double tanh function set from parameters
404      !!----------------------------------------------------------------------
405      !
406      ! Set variables from parameters
407      ! ------------------------------
408       zkth = ppkth       ;   zacr = ppacr
409       zdzmin = ppdzmin   ;   zhmax = pphmax
410       zkth2 = ppkth2     ;   zacr2 = ppacr2   ! optional (ldbletanh=T) double tanh parameters
411
412      ! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
413      !  za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
414      IF(   ppa1  == pp_to_be_computed  .AND.  &
415         &  ppa0  == pp_to_be_computed  .AND.  &
416         &  ppsur == pp_to_be_computed           ) THEN
417         !
418         za1  = (  ppdzmin - pphmax / FLOAT(jpkm1)  )                                                      &
419            & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * (  LOG( COSH( (jpk - ppkth) / ppacr) )      &
420            &                                                   - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
421         za0  = ppdzmin - za1 *              TANH( (1-ppkth) / ppacr )
422         zsur =   - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
423      ELSE
424         za1 = ppa1 ;       za0 = ppa0 ;          zsur = ppsur
425         za2 = ppa2                            ! optional (ldbletanh=T) double tanh parameter
426      ENDIF
427
428      IF(lwp) THEN                         ! Parameter print
429         WRITE(numout,*)
430         WRITE(numout,*) '    zgr_z   : Reference vertical z-coordinates'
431         WRITE(numout,*) '    ~~~~~~~'
432         IF(  ppkth == 0._wp ) THEN             
433              WRITE(numout,*) '            Uniform grid with ',jpk-1,' layers'
434              WRITE(numout,*) '            Total depth    :', zhmax
435              WRITE(numout,*) '            Layer thickness:', zhmax/(jpk-1)
436         ELSE
437            IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN
438               WRITE(numout,*) '         zsur, za0, za1 computed from '
439               WRITE(numout,*) '                 zdzmin = ', zdzmin
440               WRITE(numout,*) '                 zhmax  = ', zhmax
441            ENDIF
442            WRITE(numout,*) '           Value of coefficients for vertical mesh:'
443            WRITE(numout,*) '                 zsur = ', zsur
444            WRITE(numout,*) '                 za0  = ', za0
445            WRITE(numout,*) '                 za1  = ', za1
446            WRITE(numout,*) '                 zkth = ', zkth
447            WRITE(numout,*) '                 zacr = ', zacr
448            IF( ldbletanh ) THEN
449               WRITE(numout,*) ' (Double tanh    za2  = ', za2
450               WRITE(numout,*) '  parameters)    zkth2= ', zkth2
451               WRITE(numout,*) '                 zacr2= ', zacr2
452            ENDIF
453         ENDIF
454      ENDIF
455
456
457      ! Reference z-coordinate (depth - scale factor at T- and W-points)
458      ! ======================
459      IF( ppkth == 0._wp ) THEN            !  uniform vertical grid
460
461
462
463         za1 = zhmax / FLOAT(jpk-1) 
464
465         DO jk = 1, jpk
466            zw = FLOAT( jk )
467            zt = FLOAT( jk ) + 0.5_wp
468            gdepw_1d(jk) = ( zw - 1 ) * za1
469            gdept_1d(jk) = ( zt - 1 ) * za1
470            e3w_1d  (jk) =  za1
471            e3t_1d  (jk) =  za1
472         END DO
473      ELSE                                ! Madec & Imbard 1996 function
474         IF( .NOT. ldbletanh ) THEN
475            DO jk = 1, jpk
476               zw = REAL( jk , wp )
477               zt = REAL( jk , wp ) + 0.5_wp
478               gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) )  )
479               gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) )  )
480               e3w_1d  (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   )
481               e3t_1d  (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   )
482            END DO
483         ELSE
484            DO jk = 1, jpk
485               zw = FLOAT( jk )
486               zt = FLOAT( jk ) + 0.5_wp
487               ! Double tanh function
488               gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr  ) )    &
489                  &                             + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) )  )
490               gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr  ) )    &
491                  &                             + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) )  )
492               e3w_1d  (jk) =          za0      + za1        * TANH(       (zw-zkth ) / zacr  )      &
493                  &                             + za2        * TANH(       (zw-zkth2) / zacr2 )
494               e3t_1d  (jk) =          za0      + za1        * TANH(       (zt-zkth ) / zacr  )      &
495                  &                             + za2        * TANH(       (zt-zkth2) / zacr2 )
496            END DO
497         ENDIF
498         gdepw_1d(1) = 0._wp                    ! force first w-level to be exactly at zero
499      ENDIF
500
501      IF ( ln_isfcav .OR. ln_e3_dep ) THEN      ! e3. = dk[gdep]   
502         !
503         DO jk = 1, jpkm1
504            e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) 
505         END DO
506         e3t_1d(jpk) = e3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO
507
508         DO jk = 2, jpk
509            e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1) 
510         END DO
511         e3w_1d(1  ) = 2._wp * (gdept_1d(1) - gdepw_1d(1)) 
512      END IF
513
514!!gm BUG in s-coordinate this does not work!
515      ! deepest/shallowest W level Above/Below ~10m
516      zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d )                   ! ref. depth with tolerance (10% of minimum layer thickness)
517      nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
518      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m
519!!gm end bug
520
521      IF(lwp) THEN                        ! control print
522         WRITE(numout,*)
523         WRITE(numout,*) '              Reference z-coordinate depth and scale factors:'
524         WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
525         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk )
526      ENDIF
527      DO jk = 1, jpk                      ! control positivity
528         IF( e3w_1d  (jk) <= 0._wp .OR. e3t_1d  (jk) <= 0._wp )   CALL ctl_stop( 'dom:zgr_z: e3w_1d or e3t_1d =< 0 '    )
529         IF( gdepw_1d(jk) <  0._wp .OR. gdept_1d(jk) <  0._wp )   CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' )
530      END DO
531      !
532   END SUBROUTINE zgr_z
533
534
535   SUBROUTINE zgr_bat
536      !!----------------------------------------------------------------------
537      !!                    ***  ROUTINE zgr_bat  ***
538      !!
539      !! ** Purpose :   set bathymetry both in levels and meters
540      !!
541      !! ** Method  :   read or define mbathy and bathy arrays
542      !!       * level bathymetry:
543      !!      The ocean basin geometry is given by a two-dimensional array,
544      !!      mbathy, which is defined as follow :
545      !!            mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
546      !!                              at t-point (ji,jj).
547      !!                            = 0  over the continental t-point.
548      !!      The array mbathy is checked to verified its consistency with
549      !!      model option. in particular:
550      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
551      !!                  along closed boundary.
552      !!            mbathy must be cyclic IF jperio=1.
553      !!            mbathy must be lower or equal to jpk-1.
554      !!            isolated ocean grid points are suppressed from mbathy
555      !!                  since they are only connected to remaining
556      !!                  ocean through vertical diffusion.
557      !!      ntopo=-1 :   rectangular channel or bassin with a bump
558      !!      ntopo= 0 :   flat rectangular channel or basin
559      !!      ntopo= 1 :   mbathy is read in 'bathy_level.nc' NetCDF file
560      !!                   bathy  is read in 'bathy_meter.nc' NetCDF file
561      !!
562      !! ** Action  : - mbathy: level bathymetry (in level index)
563      !!              - bathy : meter bathymetry (in meters)
564      !!----------------------------------------------------------------------
565      INTEGER  ::   ji, jj, jk            ! dummy loop indices
566      INTEGER  ::   inum                      ! temporary logical unit
567      INTEGER  ::   ierror                    ! error flag
568      INTEGER  ::   ii_bump, ij_bump, ih      ! bump center position
569      INTEGER  ::   ii0, ii1, ij0, ij1, ik    ! local indices
570      REAL(wp) ::   r_bump , h_bump , h_oce   ! bump characteristics
571      REAL(wp) ::   zi, zj, zh, zhmin         ! local scalars
572      INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   idta   ! global domain integer data
573      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdta   ! global domain scalar data
574      !!----------------------------------------------------------------------
575      !
576      IF(lwp) WRITE(numout,*)
577      IF(lwp) WRITE(numout,*) '    zgr_bat : defines level and meter bathymetry'
578      IF(lwp) WRITE(numout,*) '    ~~~~~~~'
579      !                                               ! ================== !
580      IF( ntopo == 0 .OR. ntopo == -1 ) THEN          !   defined by hand  !
581         !                                            ! ================== !
582         !                                            ! global domain level and meter bathymetry (idta,zdta)
583         !
584         ALLOCATE( idta(jpiglo,jpjglo), STAT=ierror )
585         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate idta array' )
586         ALLOCATE( zdta(jpiglo,jpjglo), STAT=ierror )
587         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zdta array' )
588         !
589         IF( ntopo == 0 ) THEN                        ! flat basin
590            IF(lwp) WRITE(numout,*)
591            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin'
592            IF( rn_bathy > 0.01 ) THEN
593               IF(lwp) WRITE(numout,*) '         Depth = rn_bathy read in namelist'
594               zdta(:,:) = rn_bathy
595               IF( ln_sco ) THEN                                   ! s-coordinate (zsc       ): idta()=jpk
596                  idta(:,:) = jpkm1
597               ELSE                                                ! z-coordinate (zco or zps): step-like topography
598                  idta(:,:) = jpkm1
599                  DO jk = 1, jpkm1
600                     WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) )   idta(:,:) = jk
601                  END DO
602               ENDIF
603            ELSE
604               IF(lwp) WRITE(numout,*) '         Depth = depthw(jpkm1)'
605               idta(:,:) = jpkm1                            ! before last level
606               zdta(:,:) = gdepw_1d(jpk)                     ! last w-point depth
607               h_oce     = gdepw_1d(jpk)
608            ENDIF
609         ELSE                                         ! bump centered in the basin
610            IF(lwp) WRITE(numout,*)
611            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin with a bump'
612            ii_bump = jpiglo / 2                           ! i-index of the bump center
613            ij_bump = jpjglo / 2                           ! j-index of the bump center
614            r_bump  = 50000._wp                            ! bump radius (meters)       
615            h_bump  =  2700._wp                            ! bump height (meters)
616            h_oce   = gdepw_1d(jpk)                        ! background ocean depth (meters)
617            IF(lwp) WRITE(numout,*) '            bump characteristics: '
618            IF(lwp) WRITE(numout,*) '               bump center (i,j)   = ', ii_bump, ii_bump
619            IF(lwp) WRITE(numout,*) '               bump height         = ', h_bump , ' meters'
620            IF(lwp) WRITE(numout,*) '               bump radius         = ', r_bump , ' index'
621            IF(lwp) WRITE(numout,*) '            background ocean depth = ', h_oce  , ' meters'
622            !                                       
623            DO jj = 1, jpjglo                              ! zdta :
624               DO ji = 1, jpiglo
625                  zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
626                  zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
627                  zdta(ji,jj) = h_oce - h_bump * EXP( -( zi*zi + zj*zj ) )
628               END DO
629            END DO
630            !                                              ! idta :
631            IF( ln_sco ) THEN                                   ! s-coordinate (zsc       ): idta()=jpk
632               idta(:,:) = jpkm1
633            ELSE                                                ! z-coordinate (zco or zps): step-like topography
634               idta(:,:) = jpkm1
635               DO jk = 1, jpkm1
636                  WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) )   idta(:,:) = jk
637               END DO
638            ENDIF
639         ENDIF
640         !
641         !                                            ! set GLOBAL boundary conditions
642         IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
643            idta( :    , 1    ) = -1                ;      zdta( :    , 1    ) = -1._wp
644            idta( :    ,jpjglo) =  0                ;      zdta( :    ,jpjglo) =  0._wp
645         ELSEIF( jperio == 2 ) THEN
646            idta( :    , 1    ) = idta( : ,  3  )   ;      zdta( :    , 1    ) = zdta( : ,  3  )
647            idta( :    ,jpjglo) = 0                 ;      zdta( :    ,jpjglo) =  0._wp
648            idta( 1    , :    ) = 0                 ;      zdta( 1    , :    ) =  0._wp
649            idta(jpiglo, :    ) = 0                 ;      zdta(jpiglo, :    ) =  0._wp
650         ELSE
651            ih = 0                                  ;      zh = 0._wp
652            IF( ln_sco )   ih = jpkm1               ;      IF( ln_sco )   zh = h_oce
653            idta( :    , 1    ) = ih                ;      zdta( :    , 1    ) =  zh
654            idta( :    ,jpjglo) = ih                ;      zdta( :    ,jpjglo) =  zh
655            idta( 1    , :    ) = ih                ;      zdta( 1    , :    ) =  zh
656            idta(jpiglo, :    ) = ih                ;      zdta(jpiglo, :    ) =  zh
657         ENDIF
658
659         !                                            ! local domain level and meter bathymetries (mbathy,bathy)
660         mbathy(:,:) = 0                                   ! set to zero extra halo points
661         bathy (:,:) = 0._wp                               ! (require for mpp case)
662         DO jj = 1, nlcj                                   ! interior values
663            DO ji = 1, nlci
664               mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
665               bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
666            END DO
667         END DO
668         risfdep(:,:)=0.e0
669         misfdep(:,:)=1
670         !
671         DEALLOCATE( idta, zdta )
672         !
673         !                                            ! ================ !
674      ELSEIF( ntopo == 1 .OR. ntopo ==2 ) THEN                       !   read in file   ! (over the local domain)
675         !                                            ! ================ !
676         !
677         IF( ln_zco )   THEN                          ! zco : read level bathymetry
678            CALL iom_open ( 'bathy_level.nc', inum ) 
679            CALL iom_get  ( inum, jpdom_data, 'Bathy_level', bathy )
680            CALL iom_close( inum )
681            mbathy(:,:) = INT( bathy(:,:) )
682            ! initialisation isf variables
683            risfdep(:,:)=0._wp ; misfdep(:,:)=1             
684            !                                                ! =====================
685            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
686               !                                             ! =====================
687               !
688               ii0 = 140   ;   ii1 = 140                  ! Gibraltar Strait open
689               ij0 = 102   ;   ij1 = 102                  ! (Thomson, Ocean Modelling, 1995)
690               DO ji = mi0(ii0), mi1(ii1)
691                  DO jj = mj0(ij0), mj1(ij1)
692                     mbathy(ji,jj) = 15
693                  END DO
694               END DO
695               IF(lwp) WRITE(numout,*)
696               IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
697               !
698               ii0 = 160   ;   ii1 = 160                  ! Bab el mandeb Strait open
699               ij0 = 88    ;   ij1 = 88                   ! (Thomson, Ocean Modelling, 1995)
700               DO ji = mi0(ii0), mi1(ii1)
701                  DO jj = mj0(ij0), mj1(ij1)
702                     mbathy(ji,jj) = 12
703                  END DO
704               END DO
705               IF(lwp) WRITE(numout,*)
706               IF(lwp) WRITE(numout,*) '      orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
707               !
708            ENDIF
709            !
710         ENDIF
711         IF( ln_zps .OR. ln_sco )   THEN              ! zps or sco : read meter bathymetry
712#if defined key_agrif
713            IF (agrif_root()) THEN
714#endif
715            IF( ntopo == 1) THEN
716               CALL iom_open ( 'bathy_meter.nc', inum ) 
717               IF ( ln_isfcav ) THEN
718                  CALL iom_get  ( inum, jpdom_data, 'Bathymetry_isf', bathy, lrowattr=.false. )
719               ELSE
720                  CALL iom_get  ( inum, jpdom_data, 'Bathymetry'    , bathy, lrowattr=ln_use_jattr  )
721               END IF
722               CALL iom_close( inum )
723            ELSE
724               CALL dom_bat
725            ENDIF       
726#if defined key_agrif
727            ELSE
728               IF( ntopo == 1) THEN
729                  CALL agrif_create_bathy_meter()
730               ELSE
731                  CALL dom_bat 
732               ENDIF   
733            ENDIF
734#endif
735            !                                               
736            ! initialisation isf variables
737            risfdep(:,:)=0._wp ; misfdep(:,:)=1             
738            !
739            IF ( ln_isfcav ) THEN
740               CALL iom_open ( 'isf_draft_meter.nc', inum ) 
741               CALL iom_get  ( inum, jpdom_data, 'isf_draft', risfdep )
742               CALL iom_close( inum )
743            END IF
744            !       
745            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
746            !
747              ii0 = 140   ;   ii1 = 140                   ! Gibraltar Strait open
748              ij0 = 102   ;   ij1 = 102                   ! (Thomson, Ocean Modelling, 1995)
749              DO ji = mi0(ii0), mi1(ii1)
750                 DO jj = mj0(ij0), mj1(ij1)
751                    bathy(ji,jj) = 284._wp
752                 END DO
753               END DO
754              IF(lwp) WRITE(numout,*)     
755              IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
756              !
757              ii0 = 160   ;   ii1 = 160                   ! Bab el mandeb Strait open
758              ij0 = 88    ;   ij1 = 88                    ! (Thomson, Ocean Modelling, 1995)
759               DO ji = mi0(ii0), mi1(ii1)
760                 DO jj = mj0(ij0), mj1(ij1)
761                    bathy(ji,jj) = 137._wp
762                 END DO
763              END DO
764              IF(lwp) WRITE(numout,*)
765               IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
766              !
767           ENDIF
768            !
769        ENDIF
770         !                                            ! =============== !
771      ELSE                                            !      error      !
772         !                                            ! =============== !
773         WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
774         CALL ctl_stop( '    zgr_bat : '//trim(ctmp1) )
775      ENDIF
776      !
777      IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==!
778         IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level
779         ELSE                          ;   ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 )  ! from a depth
780         ENDIF
781         zhmin = gdepw_1d(ik+1)                                                        ! minimum depth = ik+1 w-levels
782         WHERE( bathy(:,:) <= 0._wp )   ;   bathy(:,:) = 0._wp                         ! min=0     over the lands
783         ELSE WHERE ( risfdep == 0._wp );   bathy(:,:) = MAX(  zhmin , bathy(:,:)  )   ! min=zhmin over the oceans
784         END WHERE
785         IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
786      ENDIF
787      !
788   END SUBROUTINE zgr_bat
789
790   SUBROUTINE zgr_bat_ctl
791      !!----------------------------------------------------------------------
792      !!                    ***  ROUTINE zgr_bat_ctl  ***
793      !!
794      !! ** Purpose :   check the bathymetry in levels
795      !!
796      !! ** Method  :   The array mbathy is checked to verified its consistency
797      !!      with the model options. in particular:
798      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
799      !!                  along closed boundary.
800      !!            mbathy must be cyclic IF jperio=1.
801      !!            mbathy must be lower or equal to jpk-1.
802      !!            isolated ocean grid points are suppressed from mbathy
803      !!                  since they are only connected to remaining
804      !!                  ocean through vertical diffusion.
805      !!      C A U T I O N : mbathy will be modified during the initializa-
806      !!      tion phase to become the number of non-zero w-levels of a water
807      !!      column, with a minimum value of 1.
808      !!
809      !! ** Action  : - update mbathy: level bathymetry (in level index)
810      !!              - update bathy : meter bathymetry (in meters)
811      !!----------------------------------------------------------------------
812      INTEGER ::   ji, jj, jl                    ! dummy loop indices
813      INTEGER ::   icompt, ibtest, ikmax         ! temporary integers
814      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  zbathy
815      !!----------------------------------------------------------------------
816      !
817      ALLOCATE(zbathy(jpi,jpj))
818      !
819      IF(lwp) WRITE(numout,*)
820      IF(lwp) WRITE(numout,*) '    zgr_bat_ctl : check the bathymetry'
821      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
822      !                                          ! Suppress isolated ocean grid points
823      IF(lwp) WRITE(numout,*)
824      IF(lwp) WRITE(numout,*)'                   suppress isolated ocean grid points'
825      IF(lwp) WRITE(numout,*)'                   -----------------------------------'
826      icompt = 0
827      DO jl = 1, 2
828         IF( l_Iperio ) THEN
829            mbathy( 1 ,:) = mbathy(jpim1,:)           ! local domain is cyclic east-west
830            mbathy(jpi,:) = mbathy(  2  ,:)
831         ENDIF
832         zbathy(:,:) = FLOAT( mbathy(:,:) )
833         CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp )
834         mbathy(:,:) = INT( zbathy(:,:) )
835         
836         DO jj = 2, jpjm1
837            DO ji = 2, jpim1
838               ibtest = MAX(  mbathy(ji-1,jj), mbathy(ji+1,jj),   &
839                  &           mbathy(ji,jj-1), mbathy(ji,jj+1)  )
840               IF( ibtest < mbathy(ji,jj) ) THEN
841                  IF(lwp) WRITE(numout,*) ' the number of ocean level at ',   &
842                     &   'grid-point (i,j) =  ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
843                  mbathy(ji,jj) = ibtest
844                  icompt = icompt + 1
845               ENDIF
846            END DO
847         END DO
848      END DO
849
850      IF( lk_mpp )   CALL mpp_sum( 'domzgr', icompt )
851      IF( icompt == 0 ) THEN
852         IF(lwp) WRITE(numout,*)'     no isolated ocean grid points'
853      ELSE
854         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed'
855      ENDIF
856
857      IF( lk_mpp ) THEN
858         zbathy(:,:) = FLOAT( mbathy(:,:) )
859         CALL lbc_lnk( 'toto',zbathy, 'T', 1._wp )
860         mbathy(:,:) = INT( zbathy(:,:) )
861      ENDIF
862      !                                          ! East-west cyclic boundary conditions
863      IF( jperio == 0 ) THEN
864         IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: jperio = ', jperio
865         IF( lk_mpp ) THEN
866            IF( nbondi == -1 .OR. nbondi == 2 ) THEN
867               IF( jperio /= 1 )   mbathy(1,:) = 0
868            ENDIF
869            IF( nbondi == 1 .OR. nbondi == 2 ) THEN
870               IF( jperio /= 1 )   mbathy(nlci,:) = 0
871            ENDIF
872         ELSE
873            IF( ln_zco .OR. ln_zps ) THEN
874               mbathy( 1 ,:) = 0
875               mbathy(jpi,:) = 0
876            ELSE
877               mbathy( 1 ,:) = jpkm1
878               mbathy(jpi,:) = jpkm1
879            ENDIF
880         ENDIF
881      ELSEIF( l_Iperio ) THEN
882         IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: jperio = ', jperio
883         mbathy( 1 ,:) = mbathy(jpim1,:)
884         mbathy(jpi,:) = mbathy(  2  ,:)
885      ELSEIF( jperio == 2 ) THEN
886         IF(lwp) WRITE(numout,*) '   equatorial boundary conditions on mbathy: jperio = ', jperio
887      ELSE
888         IF(lwp) WRITE(numout,*) '    e r r o r'
889         IF(lwp) WRITE(numout,*) '    parameter , jperio = ', jperio
890         !         STOP 'dom_mba'
891      ENDIF
892
893      !  Boundary condition on mbathy
894      IF( .NOT.lk_mpp ) THEN 
895!!gm     !!bug ???  think about it !
896         !   ... mono- or macro-tasking: T-point, >0, 2D array, no slab
897         zbathy(:,:) = FLOAT( mbathy(:,:) )
898         CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp )
899         mbathy(:,:) = INT( zbathy(:,:) )
900      ENDIF
901
902      ! Number of ocean level inferior or equal to jpkm1
903      zbathy(:,:) = FLOAT( mbathy(:,:) )
904      ikmax = glob_max( 'domzgr', zbathy(:,:) )
905
906      IF( ikmax > jpkm1 ) THEN
907         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' >  jpk-1'
908         IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
909      ELSE IF( ikmax < jpkm1 ) THEN
910         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1' 
911         IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
912      ENDIF
913      !
914      DEALLOCATE( zbathy )
915      !
916   END SUBROUTINE zgr_bat_ctl
917
918
919   SUBROUTINE zgr_bot_level
920      !!----------------------------------------------------------------------
921      !!                    ***  ROUTINE zgr_bot_level  ***
922      !!
923      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
924      !!
925      !! ** Method  :   computes from mbathy with a minimum value of 1 over land
926      !!
927      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest
928      !!                                     ocean level at t-, u- & v-points
929      !!                                     (min value = 1 over land)
930      !!----------------------------------------------------------------------
931      INTEGER ::   ji, jj   ! dummy loop indices
932      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  zmbk
933      !!----------------------------------------------------------------------
934      !
935      ALLOCATE( zmbk(jpi,jpj) )
936      !
937      IF(lwp) WRITE(numout,*)
938      IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
939      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
940      !
941      mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
942 
943      !                                     ! bottom k-index of W-level = mbkt+1
944      DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level
945         DO ji = 1, jpim1
946            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
947            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
948         END DO
949      END DO
950      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
951      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk('domzgr',zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
952      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk('domzgr',zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
953      !
954      DEALLOCATE( zmbk )
955      !
956   END SUBROUTINE zgr_bot_level
957
958
959   SUBROUTINE zgr_top_level
960      !!----------------------------------------------------------------------
961      !!                    ***  ROUTINE zgr_top_level  ***
962      !!
963      !! ** Purpose :   defines the vertical index of ocean top (mik. arrays)
964      !!
965      !! ** Method  :   computes from misfdep with a minimum value of 1
966      !!
967      !! ** Action  :   mikt, miku, mikv :   vertical indices of the shallowest
968      !!                                     ocean level at t-, u- & v-points
969      !!                                     (min value = 1)
970      !!----------------------------------------------------------------------
971      INTEGER ::   ji, jj   ! dummy loop indices
972      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  zmik
973      !!----------------------------------------------------------------------
974      !
975      ALLOCATE( zmik(jpi,jpj) )
976      !
977      IF(lwp) WRITE(numout,*)
978      IF(lwp) WRITE(numout,*) '    zgr_top_level : ocean top k-index of T-, U-, V- and W-levels '
979      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
980      !
981      mikt(:,:) = MAX( misfdep(:,:) , 1 )    ! top k-index of T-level (=1)
982      !                                      ! top k-index of W-level (=mikt)
983      DO jj = 1, jpjm1                       ! top k-index of U- (U-) level
984         DO ji = 1, jpim1
985            miku(ji,jj) = MAX(  mikt(ji+1,jj  ) , mikt(ji,jj)  )
986            mikv(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj)  )
987            mikf(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj), mikt(ji+1,jj  ), mikt(ji+1,jj+1)  )
988         END DO
989      END DO
990
991      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
992      zmik(:,:) = REAL( miku(:,:), wp )   ;   CALL lbc_lnk('domzgr',zmik,'U',1.)   ;   miku  (:,:) = MAX( INT( zmik(:,:) ), 1 )
993      zmik(:,:) = REAL( mikv(:,:), wp )   ;   CALL lbc_lnk('domzgr',zmik,'V',1.)   ;   mikv  (:,:) = MAX( INT( zmik(:,:) ), 1 )
994      zmik(:,:) = REAL( mikf(:,:), wp )   ;   CALL lbc_lnk('domzgr',zmik,'F',1.)   ;   mikf  (:,:) = MAX( INT( zmik(:,:) ), 1 )
995      !
996      DEALLOCATE( zmik )
997      !
998   END SUBROUTINE zgr_top_level
999
1000
1001   SUBROUTINE zgr_zco
1002      !!----------------------------------------------------------------------
1003      !!                  ***  ROUTINE zgr_zco  ***
1004      !!
1005      !! ** Purpose :   define the reference z-coordinate system
1006      !!
1007      !! ** Method  :   set 3D coord. arrays to reference 1D array
1008      !!----------------------------------------------------------------------
1009      INTEGER  ::   jk
1010      !!----------------------------------------------------------------------
1011      !
1012      DO jk = 1, jpk
1013         gdept_0(:,:,jk) = gdept_1d(jk)
1014         gdepw_0(:,:,jk) = gdepw_1d(jk)
1015         e3t_0  (:,:,jk) = e3t_1d  (jk)
1016         e3u_0  (:,:,jk) = e3t_1d  (jk)
1017         e3v_0  (:,:,jk) = e3t_1d  (jk)
1018         e3f_0  (:,:,jk) = e3t_1d  (jk)
1019         e3w_0  (:,:,jk) = e3w_1d  (jk)
1020         e3uw_0 (:,:,jk) = e3w_1d  (jk)
1021         e3vw_0 (:,:,jk) = e3w_1d  (jk)
1022      END DO
1023      !
1024   END SUBROUTINE zgr_zco
1025
1026
1027   SUBROUTINE zgr_zps
1028      !!----------------------------------------------------------------------
1029      !!                  ***  ROUTINE zgr_zps  ***
1030      !!                     
1031      !! ** Purpose :   the depth and vertical scale factor in partial step
1032      !!              reference z-coordinate case
1033      !!
1034      !! ** Method  :   Partial steps : computes the 3D vertical scale factors
1035      !!      of T-, U-, V-, W-, UW-, VW and F-points that are associated with
1036      !!      a partial step representation of bottom topography.
1037      !!
1038      !!        The reference depth of model levels is defined from an analytical
1039      !!      function the derivative of which gives the reference vertical
1040      !!      scale factors.
1041      !!        From  depth and scale factors reference, we compute there new value
1042      !!      with partial steps  on 3d arrays ( i, j, k ).
1043      !!
1044      !!              w-level: gdepw_0(i,j,k)  = gdep(k)
1045      !!                       e3w_0(i,j,k) = dk(gdep)(k)     = e3(i,j,k)
1046      !!              t-level: gdept_0(i,j,k)  = gdep(k+0.5)
1047      !!                       e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5)
1048      !!
1049      !!        With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
1050      !!      we find the mbathy index of the depth at each grid point.
1051      !!      This leads us to three cases:
1052      !!
1053      !!              - bathy = 0 => mbathy = 0
1054      !!              - 1 < mbathy < jpkm1   
1055      !!              - bathy > gdepw_0(jpk) => mbathy = jpkm1 
1056      !!
1057      !!        Then, for each case, we find the new depth at t- and w- levels
1058      !!      and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
1059      !!      and f-points.
1060      !!
1061      !!        This routine is given as an example, it must be modified
1062      !!      following the user s desiderata. nevertheless, the output as
1063      !!      well as the way to compute the model levels and scale factors
1064      !!      must be respected in order to insure second order accuracy
1065      !!      schemes.
1066      !!
1067      !!         c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives
1068      !!         - - - - - - -   gdept_0, gdepw_0 and e3. are positives
1069      !!     
1070      !!  Reference :   Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
1071      !!----------------------------------------------------------------------
1072      INTEGER  ::   ji, jj, jk       ! dummy loop indices
1073      INTEGER  ::   ik, it, ikb, ikt ! temporary integers
1074      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
1075      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
1076      REAL(wp) ::   zdiff            ! temporary scalar
1077      REAL(wp) ::   zmax             ! temporary scalar
1078      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  zprt
1079      !!---------------------------------------------------------------------
1080      !
1081      ALLOCATE( zprt(jpi,jpj,jpk) )
1082      !
1083      IF(lwp) WRITE(numout,*)
1084      IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps'
1085      IF(lwp) WRITE(numout,*) '    ~~~~~~~ '
1086      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used'
1087
1088      ! compute position of the ice shelf grounding line
1089      ! set bathy and isfdraft to 0 where grounded
1090      IF ( ln_isfcav ) CALL zgr_isf_zspace
1091
1092      ! bathymetry in level (from bathy_meter)
1093      ! ===================
1094      zmax = gdepw_1d(jpk) + e3t_1d(jpk)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) )
1095      bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat)
1096      WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0
1097      ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level
1098      END WHERE
1099
1100      ! Compute mbathy for ocean points (i.e. the number of ocean levels)
1101      ! find the number of ocean levels such that the last level thickness
1102      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where
1103      ! e3t_1d is the reference level thickness
1104      DO jk = jpkm1, 1, -1
1105         zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1106         WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1
1107      END DO
1108
1109      ! Check compatibility between bathy and iceshelf draft
1110      ! insure at least 2 wet level on the vertical under an ice shelf
1111      ! compute misfdep and adjust isf draft if needed
1112      IF ( ln_isfcav ) CALL zgr_isf_kspace
1113
1114      ! Scale factors and depth at T- and W-points
1115      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate
1116         gdept_0(:,:,jk) = gdept_1d(jk)
1117         gdepw_0(:,:,jk) = gdepw_1d(jk)
1118         e3t_0  (:,:,jk) = e3t_1d  (jk)
1119         e3w_0  (:,:,jk) = e3w_1d  (jk)
1120      END DO
1121     
1122      ! Scale factors and depth at T- and W-points
1123      DO jj = 1, jpj
1124         DO ji = 1, jpi
1125            ik = mbathy(ji,jj)
1126            IF( ik > 0 ) THEN               ! ocean point only
1127               ! max ocean level case
1128               IF( ik == jpkm1 ) THEN
1129                  zdepwp = bathy(ji,jj)
1130                  ze3tp  = bathy(ji,jj) - gdepw_1d(ik)
1131                  ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
1132                  e3t_0(ji,jj,ik  ) = ze3tp
1133                  e3t_0(ji,jj,ik+1) = ze3tp
1134                  e3w_0(ji,jj,ik  ) = ze3wp
1135                  e3w_0(ji,jj,ik+1) = ze3tp
1136                  gdepw_0(ji,jj,ik+1) = zdepwp
1137                  gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp
1138                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
1139                  !
1140               ELSE                         ! standard case
1141                  IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
1142                  ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
1143                  ENDIF
1144!gm Bug?  check the gdepw_1d
1145                  !       ... on ik
1146                  gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   &
1147                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   &
1148                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) ))
1149                  e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   & 
1150                     &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ) 
1151                  e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   &
1152                     &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
1153                  !       ... on ik+1
1154                  e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
1155                  e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
1156                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
1157               ENDIF
1158            ENDIF
1159         END DO
1160      END DO
1161      !
1162      it = 0
1163      DO jj = 1, jpj
1164         DO ji = 1, jpi
1165            ik = mbathy(ji,jj)
1166            IF( ik > 0 ) THEN               ! ocean point only
1167               e3tp (ji,jj) = e3t_0(ji,jj,ik)
1168               e3wp (ji,jj) = e3w_0(ji,jj,ik)
1169               ! test
1170               zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  )
1171               IF( zdiff <= 0._wp .AND. lwp ) THEN
1172                  it = it + 1
1173                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj
1174                  WRITE(numout,*) ' bathy = ', bathy(ji,jj)
1175                  WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
1176                  WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  )
1177               ENDIF
1178            ENDIF
1179         END DO
1180      END DO
1181      !
1182      ! compute top scale factor if ice shelf
1183      IF (ln_isfcav) CALL zps_isf
1184      !
1185      ! Scale factors and depth at U-, V-, UW and VW-points
1186      DO jk = 1, jpk                        ! initialisation to z-scale factors
1187         e3u_0 (:,:,jk) = e3t_1d(jk)
1188         e3v_0 (:,:,jk) = e3t_1d(jk)
1189         e3uw_0(:,:,jk) = e3w_1d(jk)
1190         e3vw_0(:,:,jk) = e3w_1d(jk)
1191      END DO
1192
1193      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
1194         DO jj = 1, jpjm1
1195            DO ji = 1, jpim1   ! vector opt.
1196               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
1197               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
1198               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
1199               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
1200            END DO
1201         END DO
1202      END DO
1203
1204      ! update e3uw in case only 2 cells in the water column
1205      IF ( ln_isfcav ) CALL zps_isf_e3uv_w
1206      !
1207      CALL lbc_lnk('domzgr', e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk('domzgr', e3uw_0, 'U', 1._wp )   ! lateral boundary conditions
1208      CALL lbc_lnk('domzgr', e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp )
1209      !
1210      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
1211         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk)
1212         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk)
1213         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk)
1214         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk)
1215      END DO
1216     
1217      ! Scale factor at F-point
1218      DO jk = 1, jpk                        ! initialisation to z-scale factors
1219         e3f_0(:,:,jk) = e3t_1d(jk)
1220      END DO
1221      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors
1222         DO jj = 1, jpjm1
1223            DO ji = 1, jpim1   ! vector opt.
1224               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
1225            END DO
1226         END DO
1227      END DO
1228      CALL lbc_lnk('domzgr', e3f_0, 'F', 1._wp )       ! Lateral boundary conditions
1229      !
1230      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
1231         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk)
1232      END DO
1233!!gm  bug ? :  must be a do loop with mj0,mj1
1234      !
1235      e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2
1236      e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 
1237      e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 
1238      e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 
1239      e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 
1240
1241      ! Control of the sign
1242      IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' )
1243      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' )
1244      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' )
1245      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' )
1246      !
1247      ! if in the future gde3w_0 need to be compute, use the function defined in NEMO
1248      ! for now gde3w_0 computation is removed as not an output of domcfg
1249
1250      DEALLOCATE( zprt )
1251      !
1252   END SUBROUTINE zgr_zps
1253
1254
1255   SUBROUTINE zgr_sco
1256      !!----------------------------------------------------------------------
1257      !!                  ***  ROUTINE zgr_sco  ***
1258      !!                     
1259      !! ** Purpose :   define the s-coordinate system
1260      !!
1261      !! ** Method  :   s-coordinate
1262      !!         The depth of model levels is defined as the product of an
1263      !!      analytical function by the local bathymetry, while the vertical
1264      !!      scale factors are defined as the product of the first derivative
1265      !!      of the analytical function by the bathymetry.
1266      !!      (this solution save memory as depth and scale factors are not
1267      !!      3d fields)
1268      !!          - Read bathymetry (in meters) at t-point and compute the
1269      !!         bathymetry at u-, v-, and f-points.
1270      !!            hbatu = mi( hbatt )
1271      !!            hbatv = mj( hbatt )
1272      !!            hbatf = mi( mj( hbatt ) )
1273      !!          - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
1274      !!         function and its derivative given as function.
1275      !!            z_gsigt(k) = fssig (k    )
1276      !!            z_gsigw(k) = fssig (k-0.5)
1277      !!            z_esigt(k) = fsdsig(k    )
1278      !!            z_esigw(k) = fsdsig(k-0.5)
1279      !!      Three options for stretching are give, and they can be modified
1280      !!      following the users requirements. Nevertheless, the output as
1281      !!      well as the way to compute the model levels and scale factors
1282      !!      must be respected in order to insure second order accuracy
1283      !!      schemes.
1284      !!
1285      !!      The three methods for stretching available are:
1286      !!
1287      !!           s_sh94 (Song and Haidvogel 1994)
1288      !!                a sinh/tanh function that allows sigma and stretched sigma
1289      !!
1290      !!           s_sf12 (Siddorn and Furner 2012?)
1291      !!                allows the maintenance of fixed surface and or
1292      !!                bottom cell resolutions (cf. geopotential coordinates)
1293      !!                within an analytically derived stretched S-coordinate framework.
1294      !!
1295      !!          s_tanh  (Madec et al 1996)
1296      !!                a cosh/tanh function that gives stretched coordinates       
1297      !!
1298      !!----------------------------------------------------------------------
1299      INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
1300      INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
1301      INTEGER  ::   ios                      ! Local integer output status for namelist read
1302      REAL(wp) ::   zrmax, ztaper   ! temporary scalars
1303      REAL(wp) ::   zrfact
1304      !
1305      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
1306      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
1307      !!
1308      NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
1309         &                rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
1310     !!----------------------------------------------------------------------
1311      !
1312      ALLOCATE( zenv(jpi,jpj), ztmp(jpi,jpj), zmsk(jpi,jpj), zri(jpi,jpj), zrj(jpi,jpj), zhbat(jpi,jpj) , ztmpi1(jpi,jpj), ztmpi2(jpi,jpj), ztmpj1(jpi,jpj), ztmpj2(jpi,jpj) )
1313      !
1314      REWIND( numnam_ref )              ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
1315      READ  ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901)
1316901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist', lwp )
1317
1318      REWIND( numnam_cfg )              ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters
1319      READ  ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 )
1320902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist', lwp )
1321      IF(lwm) WRITE ( numond, namzgr_sco )
1322
1323      IF(lwp) THEN                           ! control print
1324         WRITE(numout,*)
1325         WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
1326         WRITE(numout,*) '~~~~~~~~~~~'
1327         WRITE(numout,*) '   Namelist namzgr_sco'
1328         WRITE(numout,*) '     stretching coeffs '
1329         WRITE(numout,*) '        maximum depth of s-bottom surface (>0)       rn_sbot_max   = ',rn_sbot_max
1330         WRITE(numout,*) '        minimum depth of s-bottom surface (>0)       rn_sbot_min   = ',rn_sbot_min
1331         WRITE(numout,*) '        Critical depth                               rn_hc         = ',rn_hc
1332         WRITE(numout,*) '        maximum cut-off r-value allowed              rn_rmax       = ',rn_rmax
1333         WRITE(numout,*) '     Song and Haidvogel 1994 stretching              ln_s_sh94     = ',ln_s_sh94
1334         WRITE(numout,*) '        Song and Haidvogel 1994 stretching coefficients'
1335         WRITE(numout,*) '        surface control parameter (0<=rn_theta<=20)  rn_theta      = ',rn_theta
1336         WRITE(numout,*) '        bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ',rn_thetb
1337         WRITE(numout,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb
1338         WRITE(numout,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12
1339         WRITE(numout,*) '        switching to sigma (T) or Z (F) at H<Hc      ln_sigcrit    = ',ln_sigcrit
1340         WRITE(numout,*) '        Siddorn and Furner 2012 stretching coefficients'
1341         WRITE(numout,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha
1342         WRITE(numout,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold
1343         WRITE(numout,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs
1344         WRITE(numout,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a
1345         WRITE(numout,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b
1346         WRITE(numout,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b'
1347      ENDIF
1348
1349      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate
1350      hifu(:,:) = rn_sbot_min
1351      hifv(:,:) = rn_sbot_min
1352      hiff(:,:) = rn_sbot_min
1353
1354      !                                        ! set maximum ocean depth
1355      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1356
1357         DO jj = 1, jpj
1358            DO ji = 1, jpi
1359              IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1360            END DO
1361         END DO
1362      !                                        ! =============================
1363      !                                        ! Define the envelop bathymetry   (hbatt)
1364      !                                        ! =============================
1365      ! use r-value to create hybrid coordinates
1366      zenv(:,:) = bathy(:,:)
1367      !
1368      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
1369         DO jj = 1, jpj
1370            DO ji = 1, jpi
1371               IF( bathy(ji,jj) == 0._wp ) THEN
1372                  iip1 = MIN( ji+1, jpi )
1373                  ijp1 = MIN( jj+1, jpj )
1374                  iim1 = MAX( ji-1, 1 )
1375                  ijm1 = MAX( jj-1, 1 )
1376!!gm BUG fix see ticket #1617
1377                  IF( ( + bathy(iim1,ijm1) + bathy(ji,ijp1) + bathy(iip1,ijp1)              &
1378                     &  + bathy(iim1,jj  )                  + bathy(iip1,jj  )              &
1379                     &  + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1)  ) > 0._wp ) &
1380                     &    zenv(ji,jj) = rn_sbot_min
1381!!gm
1382!!gm               IF( ( bathy(iip1,jj  ) + bathy(iim1,jj  ) + bathy(ji,ijp1  ) + bathy(ji,ijm1) +         &
1383!!gm                  &  bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
1384!!gm               zenv(ji,jj) = rn_sbot_min
1385!!gm             ENDIF
1386!!gm end
1387              ENDIF
1388            END DO
1389         END DO
1390
1391      ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1392      CALL lbc_lnk( 'domzgr',zenv, 'T', 1._wp, 'no0' )
1393      !
1394      ! smooth the bathymetry (if required)
1395      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
1396      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
1397      !
1398      jl = 0
1399      zrmax = 1._wp
1400      !   
1401      !     
1402      ! set scaling factor used in reducing vertical gradients
1403      zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
1404      !
1405      ! initialise temporary evelope depth arrays
1406      ztmpi1(:,:) = zenv(:,:)
1407      ztmpi2(:,:) = zenv(:,:)
1408      ztmpj1(:,:) = zenv(:,:)
1409      ztmpj2(:,:) = zenv(:,:)
1410      !
1411      ! initialise temporary r-value arrays
1412      zri(:,:) = 1._wp
1413      zrj(:,:) = 1._wp
1414      !                                                            ! ================ !
1415      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) !  Iterative loop  !
1416         !                                                         ! ================ !
1417         jl = jl + 1
1418         zrmax = 0._wp
1419         ! we set zrmax from previous r-values (zri and zrj) first
1420         ! if set after current r-value calculation (as previously)
1421         ! we could exit DO WHILE prematurely before checking r-value
1422         ! of current zenv
1423         DO jj = 1, nlcj
1424            DO ji = 1, nlci
1425               zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
1426            END DO
1427         END DO
1428         zri(:,:) = 0._wp
1429         zrj(:,:) = 0._wp
1430         DO jj = 1, nlcj
1431            DO ji = 1, nlci
1432               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
1433               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
1434               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
1435                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
1436               END IF
1437               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
1438                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
1439               END IF
1440               IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact
1441               IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact
1442               IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact
1443               IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact
1444            END DO
1445         END DO
1446  !       IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
1447         !
1448         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax
1449         !
1450         DO jj = 1, nlcj
1451            DO ji = 1, nlci
1452               zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
1453            END DO
1454         END DO
1455         ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1456         CALL lbc_lnk( 'toto',zenv, 'T', 1._wp, 'no0' )
1457         !                                                  ! ================ !
1458      END DO                                                !     End loop     !
1459      !                                                     ! ================ !
1460      DO jj = 1, jpj
1461         DO ji = 1, jpi
1462            zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
1463         END DO
1464      END DO
1465      !
1466      ! Envelope bathymetry saved in hbatt
1467      hbatt(:,:) = zenv(:,:) 
1468      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
1469         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1470         DO jj = 1, jpj
1471            DO ji = 1, jpi
1472               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
1473               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
1474            END DO
1475         END DO
1476      ENDIF
1477      !
1478      !                                        ! ==============================
1479      !                                        !   hbatu, hbatv, hbatf fields
1480      !                                        ! ==============================
1481      IF(lwp) THEN
1482         WRITE(numout,*)
1483           WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
1484      ENDIF
1485      hbatu(:,:) = rn_sbot_min
1486      hbatv(:,:) = rn_sbot_min
1487      hbatf(:,:) = rn_sbot_min
1488      DO jj = 1, jpjm1
1489        DO ji = 1, jpim1   ! NO vector opt.
1490           hbatu(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
1491           hbatv(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
1492           hbatf(ji,jj) = 0.25_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
1493              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
1494        END DO
1495      END DO
1496
1497      !
1498      ! Apply lateral boundary condition
1499!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
1500      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk('domzgr', hbatu, 'U', 1._wp )
1501      DO jj = 1, jpj
1502         DO ji = 1, jpi
1503            IF( hbatu(ji,jj) == 0._wp ) THEN
1504               !No worries about the following line when ln_wd == .true.
1505               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
1506               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
1507            ENDIF
1508         END DO
1509      END DO
1510      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk('domzgr', hbatv, 'V', 1._wp )
1511      DO jj = 1, jpj
1512         DO ji = 1, jpi
1513            IF( hbatv(ji,jj) == 0._wp ) THEN
1514               IF( zhbat(ji,jj) == 0._wp )   hbatv(ji,jj) = rn_sbot_min
1515               IF( zhbat(ji,jj) /= 0._wp )   hbatv(ji,jj) = zhbat(ji,jj)
1516            ENDIF
1517         END DO
1518      END DO
1519      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk('domzgr', hbatf, 'F', 1._wp )
1520      DO jj = 1, jpj
1521         DO ji = 1, jpi
1522            IF( hbatf(ji,jj) == 0._wp ) THEN
1523               IF( zhbat(ji,jj) == 0._wp )   hbatf(ji,jj) = rn_sbot_min
1524               IF( zhbat(ji,jj) /= 0._wp )   hbatf(ji,jj) = zhbat(ji,jj)
1525            ENDIF
1526         END DO
1527      END DO
1528
1529!!bug:  key_helsinki a verifer
1530        hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1531        hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1532        hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1533        hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1534
1535      IF( lwp )   THEN
1536         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
1537            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
1538         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
1539            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
1540         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
1541            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1542         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
1543            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1544      ENDIF
1545!! helsinki
1546
1547      !                                            ! =======================
1548      !                                            !   s-ccordinate fields     (gdep., e3.)
1549      !                                            ! =======================
1550      !
1551      ! non-dimensional "sigma" for model level depth at w- and t-levels
1552
1553
1554!========================================================================
1555! Song and Haidvogel  1994 (ln_s_sh94=T)
1556! Siddorn and Furner 2012 (ln_sf12=T)
1557! or  tanh function       (both false)                   
1558!========================================================================
1559      IF      ( ln_s_sh94 ) THEN
1560                           CALL s_sh94()
1561      ELSE IF ( ln_s_sf12 ) THEN
1562                           CALL s_sf12()
1563      ELSE                 
1564                           CALL s_tanh()
1565      ENDIF
1566
1567      CALL lbc_lnk( 'domzgr',e3t_0 , 'T', 1._wp )
1568      CALL lbc_lnk( 'domzgr',e3u_0 , 'U', 1._wp )
1569      CALL lbc_lnk( 'domzgr',e3v_0 , 'V', 1._wp )
1570      CALL lbc_lnk( 'domzgr',e3f_0 , 'F', 1._wp )
1571      CALL lbc_lnk( 'domzgr',e3w_0 , 'W', 1._wp )
1572      CALL lbc_lnk( 'domzgr',e3uw_0, 'U', 1._wp )
1573      CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp )
1574      !
1575        WHERE( e3t_0 (:,:,:) == 0._wp )   e3t_0 (:,:,:) = 1._wp
1576        WHERE( e3u_0 (:,:,:) == 0._wp )   e3u_0 (:,:,:) = 1._wp
1577        WHERE( e3v_0 (:,:,:) == 0._wp )   e3v_0 (:,:,:) = 1._wp
1578        WHERE( e3f_0 (:,:,:) == 0._wp )   e3f_0 (:,:,:) = 1._wp
1579        WHERE( e3w_0 (:,:,:) == 0._wp )   e3w_0 (:,:,:) = 1._wp
1580        WHERE( e3uw_0(:,:,:) == 0._wp )   e3uw_0(:,:,:) = 1._wp
1581        WHERE( e3vw_0(:,:,:) == 0._wp )   e3vw_0(:,:,:) = 1._wp
1582!!
1583      ! HYBRID :
1584      DO jj = 1, jpj
1585         DO ji = 1, jpi
1586            DO jk = 1, jpkm1
1587               IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
1588            END DO
1589         END DO
1590      END DO
1591      IF(lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
1592         &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
1593
1594      IF( lwp )   THEN         ! min max values over the local domain
1595         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) )
1596         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
1597            &                          ' w ', MINVAL( gdepw_0(:,:,:) )
1598         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0  (:,:,:) ),   &
1599            &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0  (:,:,:) ),   &
1600            &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0 (:,:,:) ),   &
1601            &                          ' w ', MINVAL( e3w_0  (:,:,:) )
1602
1603         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
1604            &                          ' w ', MAXVAL( gdepw_0(:,:,:) )
1605         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0  (:,:,:) ),   &
1606            &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0  (:,:,:) ),   &
1607            &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0 (:,:,:) ),   &
1608            &                          ' w ', MAXVAL( e3w_0  (:,:,:) )
1609      ENDIF
1610      !  END DO
1611      IF(lwp) THEN                                  ! selected vertical profiles
1612         WRITE(numout,*)
1613         WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
1614         WRITE(numout,*) ' ~~~~~~  --------------------'
1615         WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
1616         WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     &
1617            &                                 e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
1618         DO jj = mj0(20), mj1(20)
1619            DO ji = mi0(20), mi1(20)
1620               WRITE(numout,*)
1621               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
1622               WRITE(numout,*) ' ~~~~~~  --------------------'
1623               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
1624               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
1625                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
1626            END DO
1627         END DO
1628         DO jj = mj0(74), mj1(74)
1629            DO ji = mi0(100), mi1(100)
1630               WRITE(numout,*)
1631               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
1632               WRITE(numout,*) ' ~~~~~~  --------------------'
1633               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
1634               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
1635                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
1636            END DO
1637         END DO
1638      ENDIF
1639      !
1640      !================================================================================
1641      ! check the coordinate makes sense
1642      !================================================================================
1643      DO ji = 1, jpi
1644         DO jj = 1, jpj
1645            !
1646            IF( hbatt(ji,jj) > 0._wp) THEN
1647               DO jk = 1, mbathy(ji,jj)
1648                 ! check coordinate is monotonically increasing
1649                 IF (e3w_0(ji,jj,jk) <= 0._wp .OR. e3t_0(ji,jj,jk) <= 0._wp ) THEN
1650                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
1651                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
1652                    WRITE(numout,*) 'e3w',e3w_0(ji,jj,:)
1653                    WRITE(numout,*) 'e3t',e3t_0(ji,jj,:)
1654                    CALL ctl_stop( ctmp1 )
1655                 ENDIF
1656                 ! and check it has never gone negative
1657                 IF( gdepw_0(ji,jj,jk) < 0._wp .OR. gdept_0(ji,jj,jk) < 0._wp ) THEN
1658                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
1659                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk
1660                    WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:)
1661                    WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
1662                    CALL ctl_stop( ctmp1 )
1663                 ENDIF
1664                 ! and check it never exceeds the total depth
1665                 IF( gdepw_0(ji,jj,jk) > hbatt(ji,jj) ) THEN
1666                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
1667                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
1668                    WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:)
1669                    CALL ctl_stop( ctmp1 )
1670                 ENDIF
1671               END DO
1672               !
1673               DO jk = 1, mbathy(ji,jj)-1
1674                 ! and check it never exceeds the total depth
1675                IF( gdept_0(ji,jj,jk) > hbatt(ji,jj) ) THEN
1676                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
1677                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
1678                    WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
1679                    CALL ctl_stop( ctmp1 )
1680                 ENDIF
1681               END DO
1682            ENDIF
1683         END DO
1684      END DO
1685      !
1686      DEALLOCATE( zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
1687      !
1688   END SUBROUTINE zgr_sco
1689
1690
1691   SUBROUTINE s_sh94()
1692      !!----------------------------------------------------------------------
1693      !!                  ***  ROUTINE s_sh94  ***
1694      !!                     
1695      !! ** Purpose :   stretch the s-coordinate system
1696      !!
1697      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994
1698      !!                mixed S/sigma coordinate
1699      !!
1700      !! Reference : Song and Haidvogel 1994.
1701      !!----------------------------------------------------------------------
1702      INTEGER  ::   ji, jj, jk           ! dummy loop argument
1703      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
1704      REAL(wp) ::   ztmpu,  ztmpv,  ztmpf
1705      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1
1706      !
1707      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3
1708      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
1709      !!----------------------------------------------------------------------
1710
1711      ALLOCATE( z_gsigw3 (jpi,jpj,jpk), z_gsigt3 (jpi,jpj,jpk) )
1712      ALLOCATE( z_esigt3 (jpi,jpj,jpk), z_esigw3 (jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk), z_esigtv3(jpi,jpj,jpk) )
1713      ALLOCATE( z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk), z_esigwv3(jpi,jpj,jpk) )
1714
1715      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp
1716      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
1717      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
1718      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
1719      !
1720      DO ji = 1, jpi
1721         DO jj = 1, jpj
1722            !
1723            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
1724               DO jk = 1, jpk
1725                  z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
1726                  z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
1727               END DO
1728            ELSE ! shallow water, uniform sigma
1729               DO jk = 1, jpk
1730                  z_gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
1731                  z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
1732                  END DO
1733            ENDIF
1734            !
1735            DO jk = 1, jpkm1
1736               z_esigt3(ji,jj,jk  ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
1737               z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
1738            END DO
1739            z_esigw3(ji,jj,1  ) = 2._wp * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) )
1740            z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
1741            !
1742            DO jk = 1, jpk
1743               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1744               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1745               gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
1746               gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
1747            END DO
1748           !
1749         END DO   ! for all jj's
1750      END DO    ! for all ji's
1751
1752      DO ji = 1, jpim1
1753         DO jj = 1, jpjm1
1754            ! extended for Wetting/Drying case
1755            ztmpu  = hbatt(ji,jj)+hbatt(ji+1,jj)
1756            ztmpv  = hbatt(ji,jj)+hbatt(ji,jj+1)
1757            ztmpf  = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
1758            ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
1759            ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
1760            ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
1761                   & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
1762            DO jk = 1, jpk
1763                 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
1764                        &              / ztmpu
1765                 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
1766                        &              / ztmpu
1767
1768                 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
1769                        &              / ztmpv
1770                 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
1771                        &              / ztmpv
1772
1773                 z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                               &
1774                        &            +   hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                               &
1775                        &            +   hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                               &
1776                        &            +   hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
1777
1778               !
1779               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1780               e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1781               e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1782               e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1783               !
1784               e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1785               e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1786               e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1787            END DO
1788        END DO
1789      END DO
1790      !
1791      DEALLOCATE( z_gsigw3, z_gsigt3                                                        )
1792      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1793      !
1794   END SUBROUTINE s_sh94
1795
1796
1797   SUBROUTINE s_sf12
1798      !!----------------------------------------------------------------------
1799      !!                  ***  ROUTINE s_sf12 ***
1800      !!                     
1801      !! ** Purpose :   stretch the s-coordinate system
1802      !!
1803      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
1804      !!                mixed S/sigma/Z coordinate
1805      !!
1806      !!                This method allows the maintenance of fixed surface and or
1807      !!                bottom cell resolutions (cf. geopotential coordinates)
1808      !!                within an analytically derived stretched S-coordinate framework.
1809      !!
1810      !!
1811      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
1812      !!----------------------------------------------------------------------
1813      INTEGER  ::   ji, jj, jk           ! dummy loop argument
1814      REAL(wp) ::   zsmth               ! smoothing around critical depth
1815      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
1816      REAL(wp) ::   ztmpu, ztmpv, ztmpf
1817      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1
1818      !
1819      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3
1820      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
1821      !!----------------------------------------------------------------------
1822      !
1823      ALLOCATE( z_gsigw3 (jpi,jpj,jpk), z_gsigt3 (jpi,jpj,jpk) )
1824      ALLOCATE( z_esigt3 (jpi,jpj,jpk), z_esigw3 (jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk), z_esigtv3(jpi,jpj,jpk))
1825      ALLOCATE( z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk), z_esigwv3(jpi,jpj,jpk) )
1826
1827      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp
1828      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
1829      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
1830      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
1831
1832      DO ji = 1, jpi
1833         DO jj = 1, jpj
1834
1835          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
1836             
1837              zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
1838                                                     ! could be changed by users but care must be taken to do so carefully
1839              zzb = 1.0_wp-(zzb/hbatt(ji,jj))
1840           
1841              zzs = rn_zs / hbatt(ji,jj) 
1842             
1843              IF (rn_efold /= 0.0_wp) THEN
1844                zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
1845              ELSE
1846                zsmth = 1.0_wp 
1847              ENDIF
1848               
1849              DO jk = 1, jpk
1850                z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
1851                z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
1852              ENDDO
1853              z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth  )
1854              z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth  )
1855 
1856          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
1857
1858            DO jk = 1, jpk
1859              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
1860              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
1861            END DO
1862
1863          ELSE  ! shallow water, z coordinates
1864
1865            DO jk = 1, jpk
1866              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
1867              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
1868            END DO
1869
1870          ENDIF
1871
1872          DO jk = 1, jpkm1
1873             z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
1874             z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
1875          END DO
1876          z_esigw3(ji,jj,1  ) = 2.0_wp * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ))
1877          z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
1878
1879          DO jk = 1, jpk
1880             gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
1881             gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
1882          END DO
1883
1884        ENDDO   ! for all jj's
1885      ENDDO    ! for all ji's
1886
1887      DO ji=1,jpi-1
1888        DO jj=1,jpj-1
1889
1890           ! extend to suit for Wetting/Drying case
1891           ztmpu  = hbatt(ji,jj)+hbatt(ji+1,jj)
1892           ztmpv  = hbatt(ji,jj)+hbatt(ji,jj+1)
1893           ztmpf  = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
1894           ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
1895           ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
1896           ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
1897                  & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
1898           DO jk = 1, jpk
1899                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
1900                       &              / ztmpu
1901                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
1902                       &              / ztmpu
1903                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
1904                       &              / ztmpv
1905                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
1906                       &              / ztmpv
1907
1908                z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                               &
1909                       &              + hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                               &
1910                       &              + hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                               &
1911                       &              + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
1912
1913             ! Code prior to wetting and drying option (for reference)
1914             !z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
1915             !                     /( hbatt(ji,jj)+hbatt(ji+1,jj) )
1916             !
1917             !z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
1918             !                     /( hbatt(ji,jj)+hbatt(ji+1,jj) )
1919             !
1920             !z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
1921             !                     /( hbatt(ji,jj)+hbatt(ji,jj+1) )
1922             !
1923             !z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
1924             !                     /( hbatt(ji,jj)+hbatt(ji,jj+1) )
1925             !
1926             !z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                                 &
1927             !                    &  +hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                                 &
1928             !                       +hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                                 &
1929             !                    &  +hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )                               &
1930             !                     /( hbatt(ji  ,jj  )+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
1931
1932             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
1933             e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
1934             e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
1935             e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
1936             !
1937             e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
1938             e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
1939             e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
1940          END DO
1941
1942        ENDDO
1943      ENDDO
1944      !
1945      CALL lbc_lnk('domzgr',e3t_0 ,'T',1.) ; CALL lbc_lnk('domzgr',e3u_0 ,'T',1.)
1946      CALL lbc_lnk('domzgr',e3v_0 ,'T',1.) ; CALL lbc_lnk('domzgr',e3f_0 ,'T',1.)
1947      CALL lbc_lnk('domzgr',e3w_0 ,'T',1.)
1948      CALL lbc_lnk('domzgr',e3uw_0,'T',1.) ; CALL lbc_lnk('domzgr',e3vw_0,'T',1.)
1949      !
1950      DEALLOCATE( z_gsigw3, z_gsigt3                                                        )
1951      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1952      !
1953   END SUBROUTINE s_sf12
1954
1955
1956   SUBROUTINE s_tanh()
1957      !!----------------------------------------------------------------------
1958      !!                  ***  ROUTINE s_tanh***
1959      !!                     
1960      !! ** Purpose :   stretch the s-coordinate system
1961      !!
1962      !! ** Method  :   s-coordinate stretch
1963      !!
1964      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
1965      !!----------------------------------------------------------------------
1966      INTEGER  ::   ji, jj, jk       ! dummy loop argument
1967      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
1968      REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_gsigw, z_gsigt
1969      REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_esigt, z_esigw
1970      !!----------------------------------------------------------------------
1971
1972      ALLOCATE( z_gsigw(jpk), z_gsigt(jpk) )
1973      ALLOCATE( z_esigt(jpk), z_esigw(jpk) )
1974
1975      z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp
1976      z_esigt  = 0._wp   ;   z_esigw  = 0._wp 
1977
1978      DO jk = 1, jpk
1979        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
1980        z_gsigt(jk) = -fssig( REAL(jk,wp)        )
1981      END DO
1982      IF( lwp )   WRITE(numout,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk)
1983      !
1984      ! Coefficients for vertical scale factors at w-, t- levels
1985!!gm bug :  define it from analytical function, not like juste bellow....
1986!!gm        or betteroffer the 2 possibilities....
1987      DO jk = 1, jpkm1
1988         z_esigt(jk  ) = z_gsigw(jk+1) - z_gsigw(jk)
1989         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
1990      END DO
1991      z_esigw( 1 ) = 2._wp * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
1992      z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
1993      !
1994      DO jk = 1, jpk
1995         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1996         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1997         gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
1998         gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
1999      END DO
2000
2001      DO jj = 1, jpj
2002         DO ji = 1, jpi
2003            DO jk = 1, jpk
2004              e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2005              e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2006              e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2007              e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
2008              !
2009              e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2010              e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2011              e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2012            END DO
2013         END DO
2014      END DO
2015      !
2016      DEALLOCATE( z_gsigw, z_gsigt )
2017      DEALLOCATE( z_esigt, z_esigw )
2018      !
2019   END SUBROUTINE s_tanh
2020
2021
2022   FUNCTION fssig( pk ) RESULT( pf )
2023      !!----------------------------------------------------------------------
2024      !!                 ***  ROUTINE fssig ***
2025      !!       
2026      !! ** Purpose :   provide the analytical function in s-coordinate
2027      !!         
2028      !! ** Method  :   the function provide the non-dimensional position of
2029      !!                T and W (i.e. between 0 and 1)
2030      !!                T-points at integer values (between 1 and jpk)
2031      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2032      !!----------------------------------------------------------------------
2033      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
2034      REAL(wp)             ::   pf   ! sigma value
2035      !!----------------------------------------------------------------------
2036      !
2037      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   &
2038         &     - TANH( rn_thetb * rn_theta                                )  )   &
2039         & * (   COSH( rn_theta                           )                      &
2040         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              &
2041         & / ( 2._wp * SINH( rn_theta ) )
2042      !
2043   END FUNCTION fssig
2044
2045
2046   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
2047      !!----------------------------------------------------------------------
2048      !!                 ***  ROUTINE fssig1 ***
2049      !!
2050      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
2051      !!
2052      !! ** Method  :   the function provides the non-dimensional position of
2053      !!                T and W (i.e. between 0 and 1)
2054      !!                T-points at integer values (between 1 and jpk)
2055      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2056      !!----------------------------------------------------------------------
2057      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
2058      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
2059      REAL(wp)             ::   pf1   ! sigma value
2060      !!----------------------------------------------------------------------
2061      !
2062      IF ( rn_theta == 0 ) then      ! uniform sigma
2063         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
2064      ELSE                        ! stretched sigma
2065         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
2066            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
2067            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
2068      ENDIF
2069      !
2070   END FUNCTION fssig1
2071
2072
2073   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
2074      !!----------------------------------------------------------------------
2075      !!                 ***  ROUTINE fgamma  ***
2076      !!
2077      !! ** Purpose :   provide analytical function for the s-coordinate
2078      !!
2079      !! ** Method  :   the function provides the non-dimensional position of
2080      !!                T and W (i.e. between 0 and 1)
2081      !!                T-points at integer values (between 1 and jpk)
2082      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2083      !!
2084      !!                This method allows the maintenance of fixed surface and or
2085      !!                bottom cell resolutions (cf. geopotential coordinates)
2086      !!                within an analytically derived stretched S-coordinate framework.
2087      !!
2088      !! Reference  :   Siddorn and Furner, in prep
2089      !!----------------------------------------------------------------------
2090      REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate
2091      REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate
2092      REAL(wp), INTENT(in   ) ::   pzb            ! Bottom box depth
2093      REAL(wp), INTENT(in   ) ::   pzs            ! surface box depth
2094      REAL(wp), INTENT(in   ) ::   psmth          ! Smoothing parameter
2095      !
2096      INTEGER  ::   jk             ! dummy loop index
2097      REAL(wp) ::   za1,za2,za3    ! local scalar
2098      REAL(wp) ::   zn1,zn2        !   -      -
2099      REAL(wp) ::   za,zb,zx       !   -      -
2100      !!----------------------------------------------------------------------
2101      !
2102      zn1  =  1._wp / REAL( jpkm1, wp )
2103      zn2  =  1._wp -  zn1
2104      !
2105      za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 
2106      za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp)
2107      za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
2108      !
2109      za = pzb - za3*(pzs-za1)-za2
2110      za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) )
2111      zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
2112      zx = 1.0_wp-za/2.0_wp-zb
2113      !
2114      DO jk = 1, jpk
2115         p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp +  &
2116            &          zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &
2117            &               (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
2118        p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth)
2119      END DO
2120      !
2121   END FUNCTION fgamma
2122
2123   !!======================================================================
2124END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.