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

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

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

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

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