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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domzgr.F90 in utils/tools/DOMAINcfg/src – NEMO

source: utils/tools/DOMAINcfg/src/domzgr.F90 @ 15348

Last change on this file since 15348 was 15348, checked in by jchanut, 8 months ago

#2638, set bathymetry minimum on child grids prior boundary update from parent.

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