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

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

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 4148

Last change on this file since 4148 was 4148, checked in by cetlod, 10 years ago

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