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 @ 15279

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

#2222 and #2638: Enable creating agrif meshes with different vertical grids (geopotential only as a start)

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 defined key_agrif
748      IF ( .NOT.Agrif_Root() ) CALL agrif_bathymetry_connect
749#endif
750      !
751      IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==!
752         IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level
753         ELSE                          ;   ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 )  ! from a depth
754         ENDIF
755         zhmin = gdepw_1d(ik+1)                                                        ! minimum depth = ik+1 w-levels
756         WHERE( bathy(:,:) <= 0._wp )   ;   bathy(:,:) = 0._wp                         ! min=0     over the lands
757         ELSE WHERE ( risfdep == 0._wp );   bathy(:,:) = MAX(  zhmin , bathy(:,:)  )   ! min=zhmin over the oceans
758         END WHERE
759         IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
760      ENDIF
761      !
762   END SUBROUTINE zgr_bat
763
764   SUBROUTINE zgr_bat_ctl
765      !!----------------------------------------------------------------------
766      !!                    ***  ROUTINE zgr_bat_ctl  ***
767      !!
768      !! ** Purpose :   check the bathymetry in levels
769      !!
770      !! ** Method  :   The array mbathy is checked to verified its consistency
771      !!      with the model options. in particular:
772      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
773      !!                  along closed boundary.
774      !!            mbathy must be cyclic IF jperio=1.
775      !!            mbathy must be lower or equal to jpk-1.
776      !!            isolated ocean grid points are suppressed from mbathy
777      !!                  since they are only connected to remaining
778      !!                  ocean through vertical diffusion.
779      !!      C A U T I O N : mbathy will be modified during the initializa-
780      !!      tion phase to become the number of non-zero w-levels of a water
781      !!      column, with a minimum value of 1.
782      !!
783      !! ** Action  : - update mbathy: level bathymetry (in level index)
784      !!              - update bathy : meter bathymetry (in meters)
785      !!----------------------------------------------------------------------
786      INTEGER ::   ji, jj, jl                    ! dummy loop indices
787      INTEGER ::   icompt, ibtest, ikmax         ! temporary integers
788      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  zbathy
789      !!----------------------------------------------------------------------
790      !
791      ALLOCATE(zbathy(jpi,jpj))
792      !
793      IF(lwp) WRITE(numout,*)
794      IF(lwp) WRITE(numout,*) '    zgr_bat_ctl : check the bathymetry'
795      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
796      !                                          ! Suppress isolated ocean grid points
797      IF(lwp) WRITE(numout,*)
798      IF(lwp) WRITE(numout,*)'                   suppress isolated ocean grid points'
799      IF(lwp) WRITE(numout,*)'                   -----------------------------------'
800      icompt = 0
801      DO jl = 1, 2
802         IF( l_Iperio ) THEN
803            mbathy( 1 ,:) = mbathy(jpim1,:)           ! local domain is cyclic east-west
804            mbathy(jpi,:) = mbathy(  2  ,:)
805         ENDIF
806         zbathy(:,:) = FLOAT( mbathy(:,:) )
807         CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp )
808         mbathy(:,:) = INT( zbathy(:,:) )
809         
810         DO jj = 2, jpjm1
811            DO ji = 2, jpim1
812               ibtest = MAX(  mbathy(ji-1,jj), mbathy(ji+1,jj),   &
813                  &           mbathy(ji,jj-1), mbathy(ji,jj+1)  )
814               IF( ibtest < mbathy(ji,jj) ) THEN
815                  IF(lwp) WRITE(numout,*) ' the number of ocean level at ',   &
816                     &   'grid-point (i,j) =  ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
817                  mbathy(ji,jj) = ibtest
818                  icompt = icompt + 1
819               ENDIF
820            END DO
821         END DO
822      END DO
823
824      IF( lk_mpp )   CALL mpp_sum( 'domzgr', icompt )
825      IF( icompt == 0 ) THEN
826         IF(lwp) WRITE(numout,*)'     no isolated ocean grid points'
827      ELSE
828         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed'
829      ENDIF
830
831      IF( lk_mpp ) THEN
832         zbathy(:,:) = FLOAT( mbathy(:,:) )
833         CALL lbc_lnk( 'toto',zbathy, 'T', 1._wp )
834         mbathy(:,:) = INT( zbathy(:,:) )
835      ENDIF
836      !                                          ! East-west cyclic boundary conditions
837
838      IF( jperio == 0 ) THEN
839         IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: jperio = ', jperio
840         IF( ln_zco .OR. ln_zps ) THEN
841           mbathy(  mi0(     1+nn_hls):mi1(     1+nn_hls),:) = 0
842           mbathy(  mi0(jpiglo-nn_hls):mi1(jpiglo-nn_hls),:) = 0
843         ELSE
844           mbathy(  mi0(     1+nn_hls):mi1(     1+nn_hls),:) = jpkm1
845           mbathy(  mi0(jpiglo-nn_hls):mi1(jpiglo-nn_hls),:) = jpkm1
846         ENDIF
847      ELSEIF( jperio == 1 .OR. jperio == 4 .OR. jperio ==  6 ) THEN
848         IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: jperio = ', jperio
849!         mbathy( 1 ,:) = mbathy(jpim1,:)
850!         mbathy(jpi,:) = mbathy(  2  ,:)
851      ELSEIF( jperio == 2 ) THEN
852         IF(lwp) WRITE(numout,*) '   equatorial boundary conditions on mbathy: jperio = ', jperio
853      ELSE
854         IF(lwp) WRITE(numout,*) '    e r r o r'
855         IF(lwp) WRITE(numout,*) '    parameter , jperio = ', jperio
856         !         STOP 'dom_mba'
857      ENDIF
858
859      !  Boundary condition on mbathy
860      IF( .NOT.lk_mpp ) THEN 
861!!gm     !!bug ???  think about it !
862         !   ... mono- or macro-tasking: T-point, >0, 2D array, no slab
863         zbathy(:,:) = FLOAT( mbathy(:,:) )
864         CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp )
865         mbathy(:,:) = INT( zbathy(:,:) )
866      ENDIF
867
868      ! Number of ocean level inferior or equal to jpkm1
869      zbathy(:,:) = FLOAT( mbathy(:,:) )
870      ikmax = MAXVAL(zbathy(:,:))
871      CALL mpp_max( 'domzgr',ikmax)
872
873      IF( ikmax > jpkm1 ) THEN
874         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' >  jpk-1'
875         IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
876      ELSE IF( ikmax < jpkm1 ) THEN
877         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1' 
878         IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
879      ENDIF
880      !
881      DEALLOCATE( zbathy )
882      !
883   END SUBROUTINE zgr_bat_ctl
884
885
886   SUBROUTINE zgr_bot_level
887      !!----------------------------------------------------------------------
888      !!                    ***  ROUTINE zgr_bot_level  ***
889      !!
890      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
891      !!
892      !! ** Method  :   computes from mbathy with a minimum value of 1 over land
893      !!
894      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest
895      !!                                     ocean level at t-, u- & v-points
896      !!                                     (min value = 1 over land)
897      !!----------------------------------------------------------------------
898      INTEGER ::   ji, jj   ! dummy loop indices
899      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  zmbk
900      !!----------------------------------------------------------------------
901      !
902      ALLOCATE( zmbk(jpi,jpj) )
903      !
904      IF(lwp) WRITE(numout,*)
905      IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
906      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
907      !
908      mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
909 
910      !                                     ! bottom k-index of W-level = mbkt+1
911      DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level
912         DO ji = 1, jpim1
913            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
914            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
915         END DO
916      END DO
917      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
918      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk('domzgr',zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
919      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk('domzgr',zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
920      !
921      DEALLOCATE( zmbk )
922      !
923   END SUBROUTINE zgr_bot_level
924
925
926   SUBROUTINE zgr_top_level
927      !!----------------------------------------------------------------------
928      !!                    ***  ROUTINE zgr_top_level  ***
929      !!
930      !! ** Purpose :   defines the vertical index of ocean top (mik. arrays)
931      !!
932      !! ** Method  :   computes from misfdep with a minimum value of 1
933      !!
934      !! ** Action  :   mikt, miku, mikv :   vertical indices of the shallowest
935      !!                                     ocean level at t-, u- & v-points
936      !!                                     (min value = 1)
937      !!----------------------------------------------------------------------
938      INTEGER ::   ji, jj   ! dummy loop indices
939      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  zmik
940      !!----------------------------------------------------------------------
941      !
942      ALLOCATE( zmik(jpi,jpj) )
943      !
944      IF(lwp) WRITE(numout,*)
945      IF(lwp) WRITE(numout,*) '    zgr_top_level : ocean top k-index of T-, U-, V- and W-levels '
946      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
947      !
948      mikt(:,:) = MAX( misfdep(:,:) , 1 )    ! top k-index of T-level (=1)
949      !                                      ! top k-index of W-level (=mikt)
950      DO jj = 1, jpjm1                       ! top k-index of U- (U-) level
951         DO ji = 1, jpim1
952            miku(ji,jj) = MAX(  mikt(ji+1,jj  ) , mikt(ji,jj)  )
953            mikv(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj)  )
954            mikf(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj), mikt(ji+1,jj  ), mikt(ji+1,jj+1)  )
955         END DO
956      END DO
957
958      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
959      zmik(:,:) = REAL( miku(:,:), wp )   ;   CALL lbc_lnk('domzgr',zmik,'U',1.)   ;   miku  (:,:) = MAX( INT( zmik(:,:) ), 1 )
960      zmik(:,:) = REAL( mikv(:,:), wp )   ;   CALL lbc_lnk('domzgr',zmik,'V',1.)   ;   mikv  (:,:) = MAX( INT( zmik(:,:) ), 1 )
961      zmik(:,:) = REAL( mikf(:,:), wp )   ;   CALL lbc_lnk('domzgr',zmik,'F',1.)   ;   mikf  (:,:) = MAX( INT( zmik(:,:) ), 1 )
962      !
963      DEALLOCATE( zmik )
964      !
965   END SUBROUTINE zgr_top_level
966
967
968   SUBROUTINE zgr_zco
969      !!----------------------------------------------------------------------
970      !!                  ***  ROUTINE zgr_zco  ***
971      !!
972      !! ** Purpose :   define the reference z-coordinate system
973      !!
974      !! ** Method  :   set 3D coord. arrays to reference 1D array
975      !!----------------------------------------------------------------------
976      INTEGER  ::   jk
977      !!----------------------------------------------------------------------
978      !
979      DO jk = 1, jpk
980         gdept_0(:,:,jk) = gdept_1d(jk)
981         gdepw_0(:,:,jk) = gdepw_1d(jk)
982         e3t_0  (:,:,jk) = e3t_1d  (jk)
983         e3u_0  (:,:,jk) = e3t_1d  (jk)
984         e3v_0  (:,:,jk) = e3t_1d  (jk)
985         e3f_0  (:,:,jk) = e3t_1d  (jk)
986         e3w_0  (:,:,jk) = e3w_1d  (jk)
987         e3uw_0 (:,:,jk) = e3w_1d  (jk)
988         e3vw_0 (:,:,jk) = e3w_1d  (jk)
989      END DO
990      !
991   END SUBROUTINE zgr_zco
992
993
994   SUBROUTINE zgr_zps
995      !!----------------------------------------------------------------------
996      !!                  ***  ROUTINE zgr_zps  ***
997      !!                     
998      !! ** Purpose :   the depth and vertical scale factor in partial step
999      !!              reference z-coordinate case
1000      !!
1001      !! ** Method  :   Partial steps : computes the 3D vertical scale factors
1002      !!      of T-, U-, V-, W-, UW-, VW and F-points that are associated with
1003      !!      a partial step representation of bottom topography.
1004      !!
1005      !!        The reference depth of model levels is defined from an analytical
1006      !!      function the derivative of which gives the reference vertical
1007      !!      scale factors.
1008      !!        From  depth and scale factors reference, we compute there new value
1009      !!      with partial steps  on 3d arrays ( i, j, k ).
1010      !!
1011      !!              w-level: gdepw_0(i,j,k)  = gdep(k)
1012      !!                       e3w_0(i,j,k) = dk(gdep)(k)     = e3(i,j,k)
1013      !!              t-level: gdept_0(i,j,k)  = gdep(k+0.5)
1014      !!                       e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5)
1015      !!
1016      !!        With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
1017      !!      we find the mbathy index of the depth at each grid point.
1018      !!      This leads us to three cases:
1019      !!
1020      !!              - bathy = 0 => mbathy = 0
1021      !!              - 1 < mbathy < jpkm1   
1022      !!              - bathy > gdepw_0(jpk) => mbathy = jpkm1 
1023      !!
1024      !!        Then, for each case, we find the new depth at t- and w- levels
1025      !!      and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
1026      !!      and f-points.
1027      !!
1028      !!        This routine is given as an example, it must be modified
1029      !!      following the user s desiderata. nevertheless, the output as
1030      !!      well as the way to compute the model levels and scale factors
1031      !!      must be respected in order to insure second order accuracy
1032      !!      schemes.
1033      !!
1034      !!         c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives
1035      !!         - - - - - - -   gdept_0, gdepw_0 and e3. are positives
1036      !!     
1037      !!  Reference :   Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
1038      !!----------------------------------------------------------------------
1039      INTEGER  ::   ji, jj, jk       ! dummy loop indices
1040      INTEGER  ::   ik, it, ikb, ikt ! temporary integers
1041      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
1042      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
1043      REAL(wp) ::   zdiff            ! temporary scalar
1044      REAL(wp) ::   zmax             ! temporary scalar
1045      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  zprt
1046      !!---------------------------------------------------------------------
1047      !
1048      ALLOCATE( zprt(jpi,jpj,jpk) )
1049      !
1050      IF(lwp) WRITE(numout,*)
1051      IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps'
1052      IF(lwp) WRITE(numout,*) '    ~~~~~~~ '
1053      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used'
1054
1055      ! compute position of the ice shelf grounding line
1056      ! set bathy and isfdraft to 0 where grounded
1057      IF ( ln_isfcav ) CALL zgr_isf_zspace
1058
1059      ! bathymetry in level (from bathy_meter)
1060      ! ===================
1061      zmax = gdepw_1d(jpk) + e3t_1d(jpk)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) )
1062      bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat)
1063      WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0
1064      ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level
1065      END WHERE
1066
1067      ! Compute mbathy for ocean points (i.e. the number of ocean levels)
1068      ! find the number of ocean levels such that the last level thickness
1069      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where
1070      ! e3t_1d is the reference level thickness
1071      DO jk = jpkm1, 1, -1
1072         zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1073         WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1
1074      END DO
1075
1076      ! Check compatibility between bathy and iceshelf draft
1077      ! insure at least 2 wet level on the vertical under an ice shelf
1078      ! compute misfdep and adjust isf draft if needed
1079      IF ( ln_isfcav ) CALL zgr_isf_kspace
1080
1081      ! Scale factors and depth at T- and W-points
1082      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate
1083         gdept_0(:,:,jk) = gdept_1d(jk)
1084         gdepw_0(:,:,jk) = gdepw_1d(jk)
1085         e3t_0  (:,:,jk) = e3t_1d  (jk)
1086         e3w_0  (:,:,jk) = e3w_1d  (jk)
1087      END DO
1088     
1089      ! Scale factors and depth at T- and W-points
1090      DO jj = 1, jpj
1091         DO ji = 1, jpi
1092            ik = mbathy(ji,jj)
1093            IF( ik > 0 ) THEN               ! ocean point only
1094               ! max ocean level case
1095               IF( ik == jpkm1 ) THEN
1096                  zdepwp = bathy(ji,jj)
1097                  ze3tp  = bathy(ji,jj) - gdepw_1d(ik)
1098                  ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
1099                  e3t_0(ji,jj,ik  ) = ze3tp
1100                  e3t_0(ji,jj,ik+1) = ze3tp
1101                  e3w_0(ji,jj,ik  ) = ze3wp
1102                  e3w_0(ji,jj,ik+1) = ze3tp
1103                  gdepw_0(ji,jj,ik+1) = zdepwp
1104                  gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp
1105                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
1106                  !
1107               ELSE                         ! standard case
1108                  IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
1109                  ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
1110                  ENDIF
1111!gm Bug?  check the gdepw_1d
1112                  !       ... on ik
1113                  gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   &
1114                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   &
1115                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) ))
1116                  e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   & 
1117                     &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ) 
1118                  e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   &
1119                     &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
1120                  !       ... on ik+1
1121                  e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
1122                  e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
1123                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
1124               ENDIF
1125            ENDIF
1126         END DO
1127      END DO
1128      !
1129      it = 0
1130      DO jj = 1, jpj
1131         DO ji = 1, jpi
1132            ik = mbathy(ji,jj)
1133            IF( ik > 0 ) THEN               ! ocean point only
1134               e3tp (ji,jj) = e3t_0(ji,jj,ik)
1135               e3wp (ji,jj) = e3w_0(ji,jj,ik)
1136               ! test
1137               zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  )
1138               IF( zdiff <= 0._wp .AND. lwp ) THEN
1139                  it = it + 1
1140                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj
1141                  WRITE(numout,*) ' bathy = ', bathy(ji,jj)
1142                  WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
1143                  WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  )
1144               ENDIF
1145            ENDIF
1146         END DO
1147      END DO
1148      !
1149      ! compute top scale factor if ice shelf
1150      IF (ln_isfcav) CALL zps_isf
1151      !
1152      ! Scale factors and depth at U-, V-, UW and VW-points
1153      DO jk = 1, jpk                        ! initialisation to z-scale factors
1154         e3u_0 (:,:,jk) = e3t_1d(jk)
1155         e3v_0 (:,:,jk) = e3t_1d(jk)
1156         e3uw_0(:,:,jk) = e3w_1d(jk)
1157         e3vw_0(:,:,jk) = e3w_1d(jk)
1158      END DO
1159
1160      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
1161         DO jj = 1, jpjm1
1162            DO ji = 1, jpim1   ! vector opt.
1163               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
1164               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
1165               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
1166               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
1167            END DO
1168         END DO
1169      END DO
1170
1171      ! update e3uw in case only 2 cells in the water column
1172      IF ( ln_isfcav ) CALL zps_isf_e3uv_w
1173      !
1174      CALL lbc_lnk('domzgr', e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk('domzgr', e3uw_0, 'U', 1._wp )   ! lateral boundary conditions
1175      CALL lbc_lnk('domzgr', e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp )
1176      !
1177      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
1178         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk)
1179         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk)
1180         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk)
1181         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk)
1182      END DO
1183     
1184      ! Scale factor at F-point
1185      DO jk = 1, jpk                        ! initialisation to z-scale factors
1186         e3f_0(:,:,jk) = e3t_1d(jk)
1187      END DO
1188      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors
1189         DO jj = 1, jpjm1
1190            DO ji = 1, jpim1   ! vector opt.
1191               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
1192            END DO
1193         END DO
1194      END DO
1195      CALL lbc_lnk('domzgr', e3f_0, 'F', 1._wp )       ! Lateral boundary conditions
1196      !
1197      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
1198         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk)
1199      END DO
1200!!gm  bug ? :  must be a do loop with mj0,mj1
1201      !
1202      e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2
1203      e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 
1204      e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 
1205      e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 
1206      e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 
1207
1208      ! Control of the sign
1209      IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' )
1210      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' )
1211      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' )
1212      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' )
1213      !
1214      ! if in the future gde3w_0 need to be compute, use the function defined in NEMO
1215      ! for now gde3w_0 computation is removed as not an output of domcfg
1216
1217      DEALLOCATE( zprt )
1218      !
1219   END SUBROUTINE zgr_zps
1220
1221
1222   SUBROUTINE zgr_sco
1223      !!----------------------------------------------------------------------
1224      !!                  ***  ROUTINE zgr_sco  ***
1225      !!                     
1226      !! ** Purpose :   define the s-coordinate system
1227      !!
1228      !! ** Method  :   s-coordinate
1229      !!         The depth of model levels is defined as the product of an
1230      !!      analytical function by the local bathymetry, while the vertical
1231      !!      scale factors are defined as the product of the first derivative
1232      !!      of the analytical function by the bathymetry.
1233      !!      (this solution save memory as depth and scale factors are not
1234      !!      3d fields)
1235      !!          - Read bathymetry (in meters) at t-point and compute the
1236      !!         bathymetry at u-, v-, and f-points.
1237      !!            hbatu = mi( hbatt )
1238      !!            hbatv = mj( hbatt )
1239      !!            hbatf = mi( mj( hbatt ) )
1240      !!          - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
1241      !!         function and its derivative given as function.
1242      !!            z_gsigt(k) = fssig (k    )
1243      !!            z_gsigw(k) = fssig (k-0.5)
1244      !!            z_esigt(k) = fsdsig(k    )
1245      !!            z_esigw(k) = fsdsig(k-0.5)
1246      !!      Three options for stretching are give, and they can be modified
1247      !!      following the users requirements. Nevertheless, the output as
1248      !!      well as the way to compute the model levels and scale factors
1249      !!      must be respected in order to insure second order accuracy
1250      !!      schemes.
1251      !!
1252      !!      The three methods for stretching available are:
1253      !!
1254      !!           s_sh94 (Song and Haidvogel 1994)
1255      !!                a sinh/tanh function that allows sigma and stretched sigma
1256      !!
1257      !!           s_sf12 (Siddorn and Furner 2012?)
1258      !!                allows the maintenance of fixed surface and or
1259      !!                bottom cell resolutions (cf. geopotential coordinates)
1260      !!                within an analytically derived stretched S-coordinate framework.
1261      !!
1262      !!          s_tanh  (Madec et al 1996)
1263      !!                a cosh/tanh function that gives stretched coordinates       
1264      !!
1265      !!----------------------------------------------------------------------
1266      INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
1267      INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
1268      INTEGER  ::   ios                      ! Local integer output status for namelist read
1269      REAL(wp) ::   zrmax, ztaper   ! temporary scalars
1270      REAL(wp) ::   zrfact
1271      !
1272      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
1273      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
1274      !!
1275      NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
1276         &                rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
1277     !!----------------------------------------------------------------------
1278      !
1279      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) )
1280      !
1281      !REWIND( numnam_ref )              ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
1282      READ  ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901)
1283901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist')
1284
1285      !REWIND( numnam_cfg )              ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters
1286      READ  ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 )
1287902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist')
1288      IF(lwm) WRITE ( numond, namzgr_sco )
1289
1290      IF(lwp) THEN                           ! control print
1291         WRITE(numout,*)
1292         WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
1293         WRITE(numout,*) '~~~~~~~~~~~'
1294         WRITE(numout,*) '   Namelist namzgr_sco'
1295         WRITE(numout,*) '     stretching coeffs '
1296         WRITE(numout,*) '        maximum depth of s-bottom surface (>0)       rn_sbot_max   = ',rn_sbot_max
1297         WRITE(numout,*) '        minimum depth of s-bottom surface (>0)       rn_sbot_min   = ',rn_sbot_min
1298         WRITE(numout,*) '        Critical depth                               rn_hc         = ',rn_hc
1299         WRITE(numout,*) '        maximum cut-off r-value allowed              rn_rmax       = ',rn_rmax
1300         WRITE(numout,*) '     Song and Haidvogel 1994 stretching              ln_s_sh94     = ',ln_s_sh94
1301         WRITE(numout,*) '        Song and Haidvogel 1994 stretching coefficients'
1302         WRITE(numout,*) '        surface control parameter (0<=rn_theta<=20)  rn_theta      = ',rn_theta
1303         WRITE(numout,*) '        bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ',rn_thetb
1304         WRITE(numout,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb
1305         WRITE(numout,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12
1306         WRITE(numout,*) '        switching to sigma (T) or Z (F) at H<Hc      ln_sigcrit    = ',ln_sigcrit
1307         WRITE(numout,*) '        Siddorn and Furner 2012 stretching coefficients'
1308         WRITE(numout,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha
1309         WRITE(numout,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold
1310         WRITE(numout,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs
1311         WRITE(numout,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a
1312         WRITE(numout,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b
1313         WRITE(numout,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b'
1314      ENDIF
1315
1316      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate
1317      hifu(:,:) = rn_sbot_min
1318      hifv(:,:) = rn_sbot_min
1319      hiff(:,:) = rn_sbot_min
1320
1321      !                                        ! set maximum ocean depth
1322      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1323
1324         DO jj = 1, jpj
1325            DO ji = 1, jpi
1326              IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1327            END DO
1328         END DO
1329      !                                        ! =============================
1330      !                                        ! Define the envelop bathymetry   (hbatt)
1331      !                                        ! =============================
1332      ! use r-value to create hybrid coordinates
1333      zenv(:,:) = bathy(:,:)
1334      !
1335      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
1336         DO jj = 1, jpj
1337            DO ji = 1, jpi
1338               IF( bathy(ji,jj) == 0._wp ) THEN
1339                  iip1 = MIN( ji+1, jpi )
1340                  ijp1 = MIN( jj+1, jpj )
1341                  iim1 = MAX( ji-1, 1 )
1342                  ijm1 = MAX( jj-1, 1 )
1343!!gm BUG fix see ticket #1617
1344                  IF( ( + bathy(iim1,ijm1) + bathy(ji,ijp1) + bathy(iip1,ijp1)              &
1345                     &  + bathy(iim1,jj  )                  + bathy(iip1,jj  )              &
1346                     &  + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1)  ) > 0._wp ) &
1347                     &    zenv(ji,jj) = rn_sbot_min
1348!!gm
1349!!gm               IF( ( bathy(iip1,jj  ) + bathy(iim1,jj  ) + bathy(ji,ijp1  ) + bathy(ji,ijm1) +         &
1350!!gm                  &  bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
1351!!gm               zenv(ji,jj) = rn_sbot_min
1352!!gm             ENDIF
1353!!gm end
1354              ENDIF
1355            END DO
1356         END DO
1357
1358      ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1359      CALL lbc_lnk( 'domzgr',zenv, 'T', 1._wp, kfillmode=jpfillnothing )
1360      !
1361      ! smooth the bathymetry (if required)
1362      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
1363      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
1364      !
1365      jl = 0
1366      zrmax = 1._wp
1367      !   
1368      !     
1369      ! set scaling factor used in reducing vertical gradients
1370      zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
1371      !
1372      ! initialise temporary evelope depth arrays
1373      ztmpi1(:,:) = zenv(:,:)
1374      ztmpi2(:,:) = zenv(:,:)
1375      ztmpj1(:,:) = zenv(:,:)
1376      ztmpj2(:,:) = zenv(:,:)
1377      !
1378      ! initialise temporary r-value arrays
1379      zri(:,:) = 1._wp
1380      zrj(:,:) = 1._wp
1381      !                                                            ! ================ !
1382      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) !  Iterative loop  !
1383         !                                                         ! ================ !
1384         jl = jl + 1
1385         zrmax = 0._wp
1386         ! we set zrmax from previous r-values (zri and zrj) first
1387         ! if set after current r-value calculation (as previously)
1388         ! we could exit DO WHILE prematurely before checking r-value
1389         ! of current zenv
1390          DO_2D( 0, 0, 0, 0 )
1391               zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
1392         END_2D
1393         zri(:,:) = 0._wp
1394         zrj(:,:) = 0._wp
1395          DO_2D( 0, 0, 0, 0 )
1396               iip1 = MIN( ji+1, jpi )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
1397               ijp1 = MIN( jj+1, jpj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
1398               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
1399                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
1400               END IF
1401               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
1402                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
1403               END IF
1404               IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact
1405               IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact
1406               IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact
1407               IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact
1408         END_2D
1409  !       IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
1410         !
1411         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax
1412         !
1413         DO_2D( 0, 0, 0, 0 )
1414               zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
1415         END_2D
1416         ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1417         CALL lbc_lnk( 'toto',zenv, 'T', 1._wp, kfillmode=jpfillnothing)
1418         !                                                  ! ================ !
1419      END DO                                                !     End loop     !
1420      !                                                     ! ================ !
1421      DO jj = 1, jpj
1422         DO ji = 1, jpi
1423            zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
1424         END DO
1425      END DO
1426      !
1427      ! Envelope bathymetry saved in hbatt
1428      hbatt(:,:) = zenv(:,:) 
1429      IF ((ntopo>0).AND.MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
1430         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1431         DO jj = 1, jpj
1432            DO ji = 1, jpi
1433               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
1434               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
1435            END DO
1436         END DO
1437      ENDIF
1438      !
1439      !                                        ! ==============================
1440      !                                        !   hbatu, hbatv, hbatf fields
1441      !                                        ! ==============================
1442      IF(lwp) THEN
1443         WRITE(numout,*)
1444           WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
1445      ENDIF
1446      hbatu(:,:) = rn_sbot_min
1447      hbatv(:,:) = rn_sbot_min
1448      hbatf(:,:) = rn_sbot_min
1449      DO jj = 1, jpjm1
1450        DO ji = 1, jpim1   ! NO vector opt.
1451           hbatu(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
1452           hbatv(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
1453           hbatf(ji,jj) = 0.25_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
1454              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
1455        END DO
1456      END DO
1457
1458      !
1459      ! Apply lateral boundary condition
1460!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
1461      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk('domzgr', hbatu, 'U', 1._wp )
1462      DO jj = 1, jpj
1463         DO ji = 1, jpi
1464            IF( hbatu(ji,jj) == 0._wp ) THEN
1465               !No worries about the following line when ln_wd == .true.
1466               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
1467               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
1468            ENDIF
1469         END DO
1470      END DO
1471      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk('domzgr', hbatv, 'V', 1._wp )
1472      DO jj = 1, jpj
1473         DO ji = 1, jpi
1474            IF( hbatv(ji,jj) == 0._wp ) THEN
1475               IF( zhbat(ji,jj) == 0._wp )   hbatv(ji,jj) = rn_sbot_min
1476               IF( zhbat(ji,jj) /= 0._wp )   hbatv(ji,jj) = zhbat(ji,jj)
1477            ENDIF
1478         END DO
1479      END DO
1480      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk('domzgr', hbatf, 'F', 1._wp )
1481      DO jj = 1, jpj
1482         DO ji = 1, jpi
1483            IF( hbatf(ji,jj) == 0._wp ) THEN
1484               IF( zhbat(ji,jj) == 0._wp )   hbatf(ji,jj) = rn_sbot_min
1485               IF( zhbat(ji,jj) /= 0._wp )   hbatf(ji,jj) = zhbat(ji,jj)
1486            ENDIF
1487         END DO
1488      END DO
1489
1490!!bug:  key_helsinki a verifer
1491        hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1492        hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1493        hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1494        hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1495
1496      IF( lwp )   THEN
1497         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
1498            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
1499         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
1500            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
1501         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
1502            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1503         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
1504            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1505      ENDIF
1506!! helsinki
1507
1508      !                                            ! =======================
1509      !                                            !   s-ccordinate fields     (gdep., e3.)
1510      !                                            ! =======================
1511      !
1512      ! non-dimensional "sigma" for model level depth at w- and t-levels
1513
1514
1515!========================================================================
1516! Song and Haidvogel  1994 (ln_s_sh94=T)
1517! Siddorn and Furner 2012 (ln_sf12=T)
1518! or  tanh function       (both false)                   
1519!========================================================================
1520      IF      ( ln_s_sh94 ) THEN
1521                           CALL s_sh94()
1522      ELSE IF ( ln_s_sf12 ) THEN
1523                           CALL s_sf12()
1524      ELSE                 
1525                           CALL s_tanh()
1526      ENDIF
1527
1528      CALL lbc_lnk( 'domzgr',e3t_0 , 'T', 1._wp )
1529      CALL lbc_lnk( 'domzgr',e3u_0 , 'U', 1._wp )
1530      CALL lbc_lnk( 'domzgr',e3v_0 , 'V', 1._wp )
1531      CALL lbc_lnk( 'domzgr',e3f_0 , 'F', 1._wp )
1532      CALL lbc_lnk( 'domzgr',e3w_0 , 'W', 1._wp )
1533      CALL lbc_lnk( 'domzgr',e3uw_0, 'U', 1._wp )
1534      CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp )
1535      !
1536        WHERE( e3t_0 (:,:,:) == 0._wp )   e3t_0 (:,:,:) = 1._wp
1537        WHERE( e3u_0 (:,:,:) == 0._wp )   e3u_0 (:,:,:) = 1._wp
1538        WHERE( e3v_0 (:,:,:) == 0._wp )   e3v_0 (:,:,:) = 1._wp
1539        WHERE( e3f_0 (:,:,:) == 0._wp )   e3f_0 (:,:,:) = 1._wp
1540        WHERE( e3w_0 (:,:,:) == 0._wp )   e3w_0 (:,:,:) = 1._wp
1541        WHERE( e3uw_0(:,:,:) == 0._wp )   e3uw_0(:,:,:) = 1._wp
1542        WHERE( e3vw_0(:,:,:) == 0._wp )   e3vw_0(:,:,:) = 1._wp
1543!!
1544      ! HYBRID :
1545      DO jj = 1, jpj
1546         DO ji = 1, jpi
1547            DO jk = 1, jpkm1
1548               IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
1549            END DO
1550         END DO
1551      END DO
1552
1553      WHERE (bathy(:,:)<=0) mbathy(:,:) = 0
1554
1555      IF(lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
1556         &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
1557
1558      IF( lwp )   THEN         ! min max values over the local domain
1559         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) )
1560         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
1561            &                          ' w ', MINVAL( gdepw_0(:,:,:) )
1562         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0  (:,:,:) ),   &
1563            &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0  (:,:,:) ),   &
1564            &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0 (:,:,:) ),   &
1565            &                          ' w ', MINVAL( e3w_0  (:,:,:) )
1566
1567         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
1568            &                          ' w ', MAXVAL( gdepw_0(:,:,:) )
1569         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0  (:,:,:) ),   &
1570            &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0  (:,:,:) ),   &
1571            &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0 (:,:,:) ),   &
1572            &                          ' w ', MAXVAL( e3w_0  (:,:,:) )
1573      ENDIF
1574      !  END DO
1575      IF(lwp) THEN                                  ! selected vertical profiles
1576         WRITE(numout,*)
1577         WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
1578         WRITE(numout,*) ' ~~~~~~  --------------------'
1579         WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
1580         WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     &
1581            &                                 e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
1582         DO jj = mj0(20), mj1(20)
1583            DO ji = mi0(20), mi1(20)
1584               WRITE(numout,*)
1585               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
1586               WRITE(numout,*) ' ~~~~~~  --------------------'
1587               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
1588               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
1589                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
1590            END DO
1591         END DO
1592         DO jj = mj0(74), mj1(74)
1593            DO ji = mi0(100), mi1(100)
1594               WRITE(numout,*)
1595               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
1596               WRITE(numout,*) ' ~~~~~~  --------------------'
1597               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
1598               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
1599                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
1600            END DO
1601         END DO
1602      ENDIF
1603      !
1604      !================================================================================
1605      ! check the coordinate makes sense
1606      !================================================================================
1607      DO ji = 1, jpi
1608         DO jj = 1, jpj
1609            !
1610            IF( hbatt(ji,jj) > 0._wp) THEN
1611               DO jk = 1, mbathy(ji,jj)
1612                 ! check coordinate is monotonically increasing
1613                 IF (e3w_0(ji,jj,jk) <= 0._wp .OR. e3t_0(ji,jj,jk) <= 0._wp ) THEN
1614                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
1615                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
1616                    WRITE(numout,*) 'e3w',e3w_0(ji,jj,:)
1617                    WRITE(numout,*) 'e3t',e3t_0(ji,jj,:)
1618                    CALL ctl_stop( ctmp1 )
1619                 ENDIF
1620                 ! and check it has never gone negative
1621                 IF( gdepw_0(ji,jj,jk) < 0._wp .OR. gdept_0(ji,jj,jk) < 0._wp ) THEN
1622                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
1623                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk
1624                    WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:)
1625                    WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
1626                    CALL ctl_stop( ctmp1 )
1627                 ENDIF
1628                 ! and check it never exceeds the total depth
1629                 IF( gdepw_0(ji,jj,jk) > hbatt(ji,jj) ) THEN
1630                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
1631                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
1632                    WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:)
1633                    CALL ctl_stop( ctmp1 )
1634                 ENDIF
1635               END DO
1636               !
1637               DO jk = 1, mbathy(ji,jj)-1
1638                 ! and check it never exceeds the total depth
1639                IF( gdept_0(ji,jj,jk) > hbatt(ji,jj) ) THEN
1640                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
1641                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
1642                    WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
1643                    CALL ctl_stop( ctmp1 )
1644                 ENDIF
1645               END DO
1646            ENDIF
1647         END DO
1648      END DO
1649      !
1650      DEALLOCATE( zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
1651      !
1652   END SUBROUTINE zgr_sco
1653
1654
1655   SUBROUTINE s_sh94()
1656      !!----------------------------------------------------------------------
1657      !!                  ***  ROUTINE s_sh94  ***
1658      !!                     
1659      !! ** Purpose :   stretch the s-coordinate system
1660      !!
1661      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994
1662      !!                mixed S/sigma coordinate
1663      !!
1664      !! Reference : Song and Haidvogel 1994.
1665      !!----------------------------------------------------------------------
1666      INTEGER  ::   ji, jj, jk           ! dummy loop argument
1667      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
1668      REAL(wp) ::   ztmpu,  ztmpv,  ztmpf
1669      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1
1670      !
1671      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3
1672      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
1673      !!----------------------------------------------------------------------
1674
1675      ALLOCATE( z_gsigw3 (jpi,jpj,jpk), z_gsigt3 (jpi,jpj,jpk) )
1676      ALLOCATE( z_esigt3 (jpi,jpj,jpk), z_esigw3 (jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk), z_esigtv3(jpi,jpj,jpk) )
1677      ALLOCATE( z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk), z_esigwv3(jpi,jpj,jpk) )
1678
1679      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp
1680      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
1681      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
1682      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
1683      !
1684      DO ji = 1, jpi
1685         DO jj = 1, jpj
1686            !
1687            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
1688               DO jk = 1, jpk
1689                  z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
1690                  z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
1691               END DO
1692            ELSE ! shallow water, uniform sigma
1693               DO jk = 1, jpk
1694                  z_gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
1695                  z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
1696                  END DO
1697            ENDIF
1698            !
1699            DO jk = 1, jpkm1
1700               z_esigt3(ji,jj,jk  ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
1701               z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
1702            END DO
1703            z_esigw3(ji,jj,1  ) = 2._wp * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) )
1704            z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
1705            !
1706            DO jk = 1, jpk
1707               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1708               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1709               gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
1710               gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
1711            END DO
1712           !
1713         END DO   ! for all jj's
1714      END DO    ! for all ji's
1715
1716      DO ji = 1, jpim1
1717         DO jj = 1, jpjm1
1718            ! extended for Wetting/Drying case
1719            ztmpu  = hbatt(ji,jj)+hbatt(ji+1,jj)
1720            ztmpv  = hbatt(ji,jj)+hbatt(ji,jj+1)
1721            ztmpf  = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
1722            ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
1723            ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
1724            ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
1725                   & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
1726            DO jk = 1, jpk
1727                 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
1728                        &              / ztmpu
1729                 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
1730                        &              / ztmpu
1731
1732                 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
1733                        &              / ztmpv
1734                 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
1735                        &              / ztmpv
1736
1737                 z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                               &
1738                        &            +   hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                               &
1739                        &            +   hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                               &
1740                        &            +   hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
1741
1742               !
1743               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1744               e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1745               e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1746               e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1747               !
1748               e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1749               e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1750               e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1751            END DO
1752        END DO
1753      END DO
1754      !
1755      DEALLOCATE( z_gsigw3, z_gsigt3                                                        )
1756      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1757      !
1758   END SUBROUTINE s_sh94
1759
1760
1761   SUBROUTINE s_sf12
1762      !!----------------------------------------------------------------------
1763      !!                  ***  ROUTINE s_sf12 ***
1764      !!                     
1765      !! ** Purpose :   stretch the s-coordinate system
1766      !!
1767      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
1768      !!                mixed S/sigma/Z coordinate
1769      !!
1770      !!                This method allows the maintenance of fixed surface and or
1771      !!                bottom cell resolutions (cf. geopotential coordinates)
1772      !!                within an analytically derived stretched S-coordinate framework.
1773      !!
1774      !!
1775      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
1776      !!----------------------------------------------------------------------
1777      INTEGER  ::   ji, jj, jk           ! dummy loop argument
1778      REAL(wp) ::   zsmth               ! smoothing around critical depth
1779      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
1780      REAL(wp) ::   ztmpu, ztmpv, ztmpf
1781      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1
1782      !
1783      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3
1784      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
1785      !!----------------------------------------------------------------------
1786      !
1787      ALLOCATE( z_gsigw3 (jpi,jpj,jpk), z_gsigt3 (jpi,jpj,jpk) )
1788      ALLOCATE( z_esigt3 (jpi,jpj,jpk), z_esigw3 (jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk), z_esigtv3(jpi,jpj,jpk))
1789      ALLOCATE( z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk), z_esigwv3(jpi,jpj,jpk) )
1790
1791      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp
1792      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
1793      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
1794      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
1795
1796      DO ji = 1, jpi
1797         DO jj = 1, jpj
1798
1799          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
1800             
1801              zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
1802                                                     ! could be changed by users but care must be taken to do so carefully
1803              zzb = 1.0_wp-(zzb/hbatt(ji,jj))
1804           
1805              zzs = rn_zs / hbatt(ji,jj) 
1806             
1807              IF (rn_efold /= 0.0_wp) THEN
1808                zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
1809              ELSE
1810                zsmth = 1.0_wp 
1811              ENDIF
1812               
1813              DO jk = 1, jpk
1814                z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
1815                z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
1816              ENDDO
1817              z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth  )
1818              z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth  )
1819 
1820          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
1821
1822            DO jk = 1, jpk
1823              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
1824              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
1825            END DO
1826
1827          ELSE  ! shallow water, z coordinates
1828
1829            DO jk = 1, jpk
1830              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
1831              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
1832            END DO
1833
1834          ENDIF
1835
1836          DO jk = 1, jpkm1
1837             z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
1838             z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
1839          END DO
1840          z_esigw3(ji,jj,1  ) = 2.0_wp * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ))
1841          z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
1842
1843          DO jk = 1, jpk
1844             gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
1845             gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
1846          END DO
1847
1848        ENDDO   ! for all jj's
1849      ENDDO    ! for all ji's
1850
1851      DO ji=1,jpi-1
1852        DO jj=1,jpj-1
1853
1854           ! extend to suit for Wetting/Drying case
1855           ztmpu  = hbatt(ji,jj)+hbatt(ji+1,jj)
1856           ztmpv  = hbatt(ji,jj)+hbatt(ji,jj+1)
1857           ztmpf  = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
1858           ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
1859           ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
1860           ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
1861                  & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
1862           DO jk = 1, jpk
1863                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
1864                       &              / ztmpu
1865                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
1866                       &              / ztmpu
1867                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
1868                       &              / ztmpv
1869                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
1870                       &              / ztmpv
1871
1872                z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                               &
1873                       &              + hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                               &
1874                       &              + hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                               &
1875                       &              + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
1876
1877             ! Code prior to wetting and drying option (for reference)
1878             !z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
1879             !                     /( hbatt(ji,jj)+hbatt(ji+1,jj) )
1880             !
1881             !z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
1882             !                     /( hbatt(ji,jj)+hbatt(ji+1,jj) )
1883             !
1884             !z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
1885             !                     /( hbatt(ji,jj)+hbatt(ji,jj+1) )
1886             !
1887             !z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
1888             !                     /( hbatt(ji,jj)+hbatt(ji,jj+1) )
1889             !
1890             !z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                                 &
1891             !                    &  +hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                                 &
1892             !                       +hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                                 &
1893             !                    &  +hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )                               &
1894             !                     /( hbatt(ji  ,jj  )+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
1895
1896             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
1897             e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
1898             e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
1899             e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
1900             !
1901             e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
1902             e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
1903             e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
1904          END DO
1905
1906        ENDDO
1907      ENDDO
1908      !
1909      CALL lbc_lnk('domzgr',e3t_0 ,'T',1.) ; CALL lbc_lnk('domzgr',e3u_0 ,'T',1.)
1910      CALL lbc_lnk('domzgr',e3v_0 ,'T',1.) ; CALL lbc_lnk('domzgr',e3f_0 ,'T',1.)
1911      CALL lbc_lnk('domzgr',e3w_0 ,'T',1.)
1912      CALL lbc_lnk('domzgr',e3uw_0,'T',1.) ; CALL lbc_lnk('domzgr',e3vw_0,'T',1.)
1913      !
1914      DEALLOCATE( z_gsigw3, z_gsigt3                                                        )
1915      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1916      !
1917   END SUBROUTINE s_sf12
1918
1919
1920   SUBROUTINE s_tanh()
1921      !!----------------------------------------------------------------------
1922      !!                  ***  ROUTINE s_tanh***
1923      !!                     
1924      !! ** Purpose :   stretch the s-coordinate system
1925      !!
1926      !! ** Method  :   s-coordinate stretch
1927      !!
1928      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
1929      !!----------------------------------------------------------------------
1930      INTEGER  ::   ji, jj, jk       ! dummy loop argument
1931      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
1932      REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_gsigw, z_gsigt
1933      REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_esigt, z_esigw
1934      !!----------------------------------------------------------------------
1935
1936      ALLOCATE( z_gsigw(jpk), z_gsigt(jpk) )
1937      ALLOCATE( z_esigt(jpk), z_esigw(jpk) )
1938
1939      z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp
1940      z_esigt  = 0._wp   ;   z_esigw  = 0._wp 
1941
1942      DO jk = 1, jpk
1943        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
1944        z_gsigt(jk) = -fssig( REAL(jk,wp)        )
1945      END DO
1946      IF( lwp )   WRITE(numout,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk)
1947      !
1948      ! Coefficients for vertical scale factors at w-, t- levels
1949!!gm bug :  define it from analytical function, not like juste bellow....
1950!!gm        or betteroffer the 2 possibilities....
1951      DO jk = 1, jpkm1
1952         z_esigt(jk  ) = z_gsigw(jk+1) - z_gsigw(jk)
1953         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
1954      END DO
1955      z_esigw( 1 ) = 2._wp * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
1956      z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
1957      !
1958      DO jk = 1, jpk
1959         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1960         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1961         gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
1962         gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
1963      END DO
1964
1965      DO jj = 1, jpj
1966         DO ji = 1, jpi
1967            DO jk = 1, jpk
1968              e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1969              e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1970              e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1971              e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
1972              !
1973              e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1974              e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1975              e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1976            END DO
1977         END DO
1978      END DO
1979      !
1980      DEALLOCATE( z_gsigw, z_gsigt )
1981      DEALLOCATE( z_esigt, z_esigw )
1982      !
1983   END SUBROUTINE s_tanh
1984
1985
1986   FUNCTION fssig( pk ) RESULT( pf )
1987      !!----------------------------------------------------------------------
1988      !!                 ***  ROUTINE fssig ***
1989      !!       
1990      !! ** Purpose :   provide the analytical function in s-coordinate
1991      !!         
1992      !! ** Method  :   the function provide the non-dimensional position of
1993      !!                T and W (i.e. between 0 and 1)
1994      !!                T-points at integer values (between 1 and jpk)
1995      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1996      !!----------------------------------------------------------------------
1997      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
1998      REAL(wp)             ::   pf   ! sigma value
1999      !!----------------------------------------------------------------------
2000      !
2001      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   &
2002         &     - TANH( rn_thetb * rn_theta                                )  )   &
2003         & * (   COSH( rn_theta                           )                      &
2004         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              &
2005         & / ( 2._wp * SINH( rn_theta ) )
2006      !
2007   END FUNCTION fssig
2008
2009
2010   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
2011      !!----------------------------------------------------------------------
2012      !!                 ***  ROUTINE fssig1 ***
2013      !!
2014      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
2015      !!
2016      !! ** Method  :   the function provides the non-dimensional position of
2017      !!                T and W (i.e. between 0 and 1)
2018      !!                T-points at integer values (between 1 and jpk)
2019      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2020      !!----------------------------------------------------------------------
2021      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
2022      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
2023      REAL(wp)             ::   pf1   ! sigma value
2024      !!----------------------------------------------------------------------
2025      !
2026      IF ( rn_theta == 0 ) then      ! uniform sigma
2027         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
2028      ELSE                        ! stretched sigma
2029         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
2030            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
2031            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
2032      ENDIF
2033      !
2034   END FUNCTION fssig1
2035
2036
2037   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
2038      !!----------------------------------------------------------------------
2039      !!                 ***  ROUTINE fgamma  ***
2040      !!
2041      !! ** Purpose :   provide analytical function for the s-coordinate
2042      !!
2043      !! ** Method  :   the function provides the non-dimensional position of
2044      !!                T and W (i.e. between 0 and 1)
2045      !!                T-points at integer values (between 1 and jpk)
2046      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2047      !!
2048      !!                This method allows the maintenance of fixed surface and or
2049      !!                bottom cell resolutions (cf. geopotential coordinates)
2050      !!                within an analytically derived stretched S-coordinate framework.
2051      !!
2052      !! Reference  :   Siddorn and Furner, in prep
2053      !!----------------------------------------------------------------------
2054      REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate
2055      REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate
2056      REAL(wp), INTENT(in   ) ::   pzb            ! Bottom box depth
2057      REAL(wp), INTENT(in   ) ::   pzs            ! surface box depth
2058      REAL(wp), INTENT(in   ) ::   psmth          ! Smoothing parameter
2059      !
2060      INTEGER  ::   jk             ! dummy loop index
2061      REAL(wp) ::   za1,za2,za3    ! local scalar
2062      REAL(wp) ::   zn1,zn2        !   -      -
2063      REAL(wp) ::   za,zb,zx       !   -      -
2064      !!----------------------------------------------------------------------
2065      !
2066      zn1  =  1._wp / REAL( jpkm1, wp )
2067      zn2  =  1._wp -  zn1
2068      !
2069      za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 
2070      za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp)
2071      za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
2072      !
2073      za = pzb - za3*(pzs-za1)-za2
2074      za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) )
2075      zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
2076      zx = 1.0_wp-za/2.0_wp-zb
2077      !
2078      DO jk = 1, jpk
2079         p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp +  &
2080            &          zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &
2081            &               (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
2082        p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth)
2083      END DO
2084      !
2085   END FUNCTION fgamma
2086
2087   !!======================================================================
2088END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.