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

source: branches/2013/dev_LOCEAN_CMCC_INGV_MERC_UKMO_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 4245

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

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