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

source: branches/2013/dev_MERCATOR_UKMO_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 4240

Last change on this file since 4240 was 4240, checked in by cbricaud, 10 years ago

correct an error and remove key vectop_loop in AMM cpp file

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