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

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

source: branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5727

Last change on this file since 5727 was 5727, checked in by rfurner, 9 years ago

some bug fixes for wetting and drying elements...still not working though

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