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 trunk/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 3967

Last change on this file since 3967 was 3967, checked in by cbricaud, 11 years ago

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