source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 18 months ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

File size: 130.9 KB
Line 
1MODULE domzgr
2   !!==============================================================================
3   !!                       ***  MODULE domzgr   ***
4   !! Ocean initialization : domain initialization
5   !!==============================================================================
6   !! History :  OPA  ! 1995-12  (G. Madec)  Original code : s vertical coordinate
7   !!                 ! 1997-07  (G. Madec)  lbc_lnk call
8   !!                 ! 1997-04  (J.-O. Beismann)
9   !!            8.5  ! 2002-09  (A. Bozec, G. Madec)  F90: Free form and module
10   !!             -   ! 2002-09  (A. de Miranda)  rigid-lid + islands
11   !!  NEMO      1.0  ! 2003-08  (G. Madec)  F90: Free form and module
12   !!             -   ! 2005-10  (A. Beckmann)  modifications for hybrid s-ccordinates & new stretching function
13   !!            2.0  ! 2006-04  (R. Benshila, G. Madec)  add zgr_zco
14   !!            3.0  ! 2008-06  (G. Madec)  insertion of domzgr_zps.h90 & conding style
15   !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option
16   !!            3.3  ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level
17   !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn and Furner stretching function
18   !!            3.4  ! 2012-12  (R. Bourdalle-Badie and G. Reffray)  modify C1D case 
19   !!            3.6  ! 2014-11  (P. Mathiot and C. Harris) add ice shelf capabilitye 
20   !!----------------------------------------------------------------------
21
22   !!----------------------------------------------------------------------
23   !!   dom_zgr          : defined the ocean vertical coordinate system
24   !!       zgr_bat      : bathymetry fields (levels and meters)
25   !!       zgr_bat_zoom : modify the bathymetry field if zoom domain
26   !!       zgr_bat_ctl  : check the bathymetry files
27   !!       zgr_bot_level: deepest ocean level for t-, u, and v-points
28   !!       zgr_z        : reference z-coordinate
29   !!       zgr_zco      : z-coordinate
30   !!       zgr_zps      : z-coordinate with partial steps
31   !!       zgr_sco      : s-coordinate
32   !!       fssig        : tanh stretch function
33   !!       fssig1       : Song and Haidvogel 1994 stretch function
34   !!       fgamma       : Siddorn and Furner 2012 stretching function
35   !!---------------------------------------------------------------------
36   USE oce               ! ocean variables
37   USE dom_oce           ! ocean domain
38   USE closea            ! closed seas
39   USE c1d               ! 1D vertical configuration
40   USE in_out_manager    ! I/O manager
41   USE iom               ! I/O library
42   USE lbclnk            ! ocean lateral boundary conditions (or mpp link)
43   USE lib_mpp           ! distributed memory computing library
44   USE wrk_nemo          ! Memory allocation
45   USE timing            ! Timing
46
47   IMPLICIT NONE
48   PRIVATE
49
50   PUBLIC   dom_zgr        ! called by dom_init.F90
51
52   !                              !!* Namelist namzgr_sco *
53   LOGICAL  ::   ln_s_sh94         ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T)
54   LOGICAL  ::   ln_s_sf12         ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T)
55   !
56   REAL(wp) ::   rn_sbot_min       ! minimum depth of s-bottom surface (>0) (m)
57   REAL(wp) ::   rn_sbot_max       ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
58   REAL(wp) ::   rn_rmax           ! maximum cut-off r-value allowed (0<rn_rmax<1)
59   REAL(wp) ::   rn_hc             ! Critical depth for transition from sigma to stretched coordinates
60   ! Song and Haidvogel 1994 stretching parameters
61   REAL(wp) ::   rn_theta          ! surface control parameter (0<=rn_theta<=20)
62   REAL(wp) ::   rn_thetb          ! bottom control parameter  (0<=rn_thetb<= 1)
63   REAL(wp) ::   rn_bb             ! stretching parameter
64   !                                        ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
65   ! Siddorn and Furner stretching parameters
66   LOGICAL  ::   ln_sigcrit        ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
67   REAL(wp) ::   rn_alpha          ! control parameter ( > 1 stretch towards surface, < 1 towards seabed)
68   REAL(wp) ::   rn_efold          !  efold length scale for transition to stretched coord
69   REAL(wp) ::   rn_zs             !  depth of surface grid box
70                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
71   REAL(wp) ::   rn_zb_a           !  bathymetry scaling factor for calculating Zb
72   REAL(wp) ::   rn_zb_b           !  offset for calculating Zb
73
74  !! * Substitutions
75#  include "domzgr_substitute.h90"
76#  include "vectopt_loop_substitute.h90"
77   !!----------------------------------------------------------------------
78   !! NEMO/OPA 3.3.1 , NEMO Consortium (2011)
79   !! $Id$
80   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
81   !!----------------------------------------------------------------------
82CONTAINS       
83
84   SUBROUTINE dom_zgr
85      !!----------------------------------------------------------------------
86      !!                ***  ROUTINE dom_zgr  ***
87      !!                   
88      !! ** Purpose :   set the depth of model levels and the resulting
89      !!              vertical scale factors.
90      !!
91      !! ** Method  : - reference 1D vertical coordinate (gdep._1d, e3._1d)
92      !!              - read/set ocean depth and ocean levels (bathy, mbathy)
93      !!              - vertical coordinate (gdep., e3.) depending on the
94      !!                coordinate chosen :
95      !!                   ln_zco=T   z-coordinate   
96      !!                   ln_zps=T   z-coordinate with partial steps
97      !!                   ln_zco=T   s-coordinate
98      !!
99      !! ** Action  :   define gdep., e3., mbathy and bathy
100      !!----------------------------------------------------------------------
101      INTEGER ::   ioptio, ibat   ! local integer
102      INTEGER ::   ios
103      !
104      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav
105      !!----------------------------------------------------------------------
106      !
107      IF( nn_timing == 1 )   CALL timing_start('dom_zgr')
108      !
109      REWIND( numnam_ref )              ! Namelist namzgr in reference namelist : Vertical coordinate
110      READ  ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901 )
111901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist', lwp )
112
113      REWIND( numnam_cfg )              ! Namelist namzgr in configuration namelist : Vertical coordinate
114      READ  ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 )
115902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp )
116      IF(lwm .AND. nprint > 2) WRITE ( numond, namzgr )
117
118      IF(lwp) THEN                     ! Control print
119         WRITE(numout,*)
120         WRITE(numout,*) 'dom_zgr : vertical coordinate'
121         WRITE(numout,*) '~~~~~~~'
122         WRITE(numout,*) '          Namelist namzgr : set vertical coordinate'
123         WRITE(numout,*) '             z-coordinate - full steps      ln_zco    = ', ln_zco
124         WRITE(numout,*) '             z-coordinate - partial steps   ln_zps    = ', ln_zps
125         WRITE(numout,*) '             s- or hybrid z-s-coordinate    ln_sco    = ', ln_sco
126         WRITE(numout,*) '             ice shelf cavities             ln_isfcav = ', ln_isfcav
127      ENDIF
128
129      ioptio = 0                       ! Check Vertical coordinate options
130      IF( ln_zco      )   ioptio = ioptio + 1
131      IF( ln_zps      )   ioptio = ioptio + 1
132      IF( ln_sco      )   ioptio = ioptio + 1
133      IF( ioptio /= 1 )   CALL ctl_stop( ' none or several vertical coordinate options used' )
134      !
135      ! Build the vertical coordinate system
136      ! ------------------------------------
137                          CALL zgr_z            ! Reference z-coordinate system (always called)
138                          CALL zgr_bat          ! Bathymetry fields (levels and meters)
139      IF( lk_c1d      )   CALL lbc_lnk( bathy , 'T', 1._wp )   ! 1D config.: same bathy value over the 3x3 domain
140      IF( ln_zco      )   CALL zgr_zco          ! z-coordinate
141      IF( ln_zps      )   CALL zgr_zps          ! Partial step z-coordinate
142      IF( ln_sco      )   CALL zgr_sco          ! s-coordinate or hybrid z-s coordinate
143      !
144      ! final adjustment of mbathy & check
145      ! -----------------------------------
146      IF( lzoom       )   CALL zgr_bat_zoom     ! correct mbathy in case of zoom subdomain
147      IF( .NOT.lk_c1d )   CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isolated ocean points
148                          CALL zgr_bot_level    ! deepest ocean level for t-, u- and v-points
149                          CALL zgr_top_level    ! shallowest ocean level for T-, U-, V- points
150      !
151      IF( lk_c1d ) THEN                         ! 1D config.: same mbathy value over the 3x3 domain
152         ibat = mbathy(2,2)
153         mbathy(:,:) = ibat
154      END IF
155      !
156      IF( nprint > 3 .AND. lwp )   THEN
157         WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
158         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
159            &                   ' w ',   MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gdep3w_0(:,:,:) )
160         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL( e3f_0(:,:,:) ),  &
161            &                   ' u ',   MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL( e3v_0(:,:,:) ),  &
162            &                   ' uw',   MINVAL( e3uw_0(:,:,:)), ' vw', MINVAL( e3vw_0(:,:,:)),   &
163            &                   ' w ',   MINVAL( e3w_0(:,:,:) )
164
165         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
166            &                   ' w ',   MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gdep3w_0(:,:,:) )
167         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL( e3f_0(:,:,:) ),  &
168            &                   ' u ',   MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL( e3v_0(:,:,:) ),  &
169            &                   ' uw',   MAXVAL( e3uw_0(:,:,:)), ' vw', MAXVAL( e3vw_0(:,:,:)),   &
170            &                   ' w ',   MAXVAL( e3w_0(:,:,:) )
171      ENDIF
172      !
173      IF(lwp .AND. lflush) CALL flush(numout)
174      !
175      IF( nn_timing == 1 )  CALL timing_stop('dom_zgr')
176      !
177   END SUBROUTINE dom_zgr
178
179
180   SUBROUTINE zgr_z
181      !!----------------------------------------------------------------------
182      !!                   ***  ROUTINE zgr_z  ***
183      !!                   
184      !! ** Purpose :   set the depth of model levels and the resulting
185      !!      vertical scale factors.
186      !!
187      !! ** Method  :   z-coordinate system (use in all type of coordinate)
188      !!        The depth of model levels is defined from an analytical
189      !!      function the derivative of which gives the scale factors.
190      !!        both depth and scale factors only depend on k (1d arrays).
191      !!              w-level: gdepw_1d  = gdep(k)
192      !!                       e3w_1d(k) = dk(gdep)(k)     = e3(k)
193      !!              t-level: gdept_1d  = gdep(k+0.5)
194      !!                       e3t_1d(k) = dk(gdep)(k+0.5) = e3(k+0.5)
195      !!
196      !! ** Action  : - gdept_1d, gdepw_1d : depth of T- and W-point (m)
197      !!              - e3t_1d  , e3w_1d   : scale factors at T- and W-levels (m)
198      !!
199      !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
200      !!----------------------------------------------------------------------
201      INTEGER  ::   jk                     ! dummy loop indices
202      REAL(wp) ::   zt, zw                 ! temporary scalars
203      REAL(wp) ::   zsur, za0, za1, zkth   ! Values set from parameters in
204      REAL(wp) ::   zacr, zdzmin, zhmax    ! par_CONFIG_Rxx.h90
205      REAL(wp) ::   zrefdep                ! depth of the reference level (~10m)
206      REAL(wp) ::   za2, zkth2, zacr2      ! Values for optional double tanh function set from parameters
207      !!----------------------------------------------------------------------
208      !
209      IF( nn_timing == 1 )  CALL timing_start('zgr_z')
210      !
211      ! Set variables from parameters
212      ! ------------------------------
213       zkth = ppkth       ;   zacr = ppacr
214       zdzmin = ppdzmin   ;   zhmax = pphmax
215       zkth2 = ppkth2     ;   zacr2 = ppacr2   ! optional (ldbletanh=T) double tanh parameters
216
217      ! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
218      !  za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
219      IF(   ppa1  == pp_to_be_computed  .AND.  &
220         &  ppa0  == pp_to_be_computed  .AND.  &
221         &  ppsur == pp_to_be_computed           ) THEN
222         !
223#if defined key_agrif
224         za1  = (  ppdzmin - pphmax / FLOAT(jpkdta-1)  )                                                   &
225            & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpkdta-1) * (  LOG( COSH( (jpkdta - ppkth) / ppacr) )&
226            &                                                      - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
227#else
228         za1  = (  ppdzmin - pphmax / FLOAT(jpkm1)  )                                                      &
229            & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * (  LOG( COSH( (jpk - ppkth) / ppacr) )      &
230            &                                                   - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
231#endif
232         za0  = ppdzmin - za1 *              TANH( (1-ppkth) / ppacr )
233         zsur =   - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
234      ELSE
235         za1 = ppa1 ;       za0 = ppa0 ;          zsur = ppsur
236         za2 = ppa2                            ! optional (ldbletanh=T) double tanh parameter
237      ENDIF
238
239      IF(lwp) THEN                         ! Parameter print
240         WRITE(numout,*)
241         WRITE(numout,*) '    zgr_z   : Reference vertical z-coordinates'
242         WRITE(numout,*) '    ~~~~~~~'
243         IF(  ppkth == 0._wp ) THEN             
244              WRITE(numout,*) '            Uniform grid with ',jpk-1,' layers'
245              WRITE(numout,*) '            Total depth    :', zhmax
246#if defined key_agrif
247              WRITE(numout,*) '            Layer thickness:', zhmax/(jpkdta-1)
248#else
249              WRITE(numout,*) '            Layer thickness:', zhmax/(jpk-1)
250#endif
251         ELSE
252            IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN
253               WRITE(numout,*) '         zsur, za0, za1 computed from '
254               WRITE(numout,*) '                 zdzmin = ', zdzmin
255               WRITE(numout,*) '                 zhmax  = ', zhmax
256            ENDIF
257            WRITE(numout,*) '           Value of coefficients for vertical mesh:'
258            WRITE(numout,*) '                 zsur = ', zsur
259            WRITE(numout,*) '                 za0  = ', za0
260            WRITE(numout,*) '                 za1  = ', za1
261            WRITE(numout,*) '                 zkth = ', zkth
262            WRITE(numout,*) '                 zacr = ', zacr
263            IF( ldbletanh ) THEN
264               WRITE(numout,*) ' (Double tanh    za2  = ', za2
265               WRITE(numout,*) '  parameters)    zkth2= ', zkth2
266               WRITE(numout,*) '                 zacr2= ', zacr2
267            ENDIF
268         ENDIF
269         IF(lflush) CALL flush(numout)
270      ENDIF
271
272
273      ! Reference z-coordinate (depth - scale factor at T- and W-points)
274      ! ======================
275      IF( ppkth == 0._wp ) THEN            !  uniform vertical grid
276#if defined key_agrif
277         za1 = zhmax / FLOAT(jpkdta-1) 
278#else
279         za1 = zhmax / FLOAT(jpk-1) 
280#endif
281         DO jk = 1, jpk
282            zw = FLOAT( jk )
283            zt = FLOAT( jk ) + 0.5_wp
284            gdepw_1d(jk) = ( zw - 1 ) * za1
285            gdept_1d(jk) = ( zt - 1 ) * za1
286            e3w_1d  (jk) =  za1
287            e3t_1d  (jk) =  za1
288         END DO
289      ELSE                                ! Madec & Imbard 1996 function
290         IF( .NOT. ldbletanh ) THEN
291            DO jk = 1, jpk
292               zw = REAL( jk , wp )
293               zt = REAL( jk , wp ) + 0.5_wp
294               gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) )  )
295               gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) )  )
296               e3w_1d  (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   )
297               e3t_1d  (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   )
298            END DO
299         ELSE
300            DO jk = 1, jpk
301               zw = FLOAT( jk )
302               zt = FLOAT( jk ) + 0.5_wp
303               ! Double tanh function
304               gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr  ) )    &
305                  &                             + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) )  )
306               gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr  ) )    &
307                  &                             + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) )  )
308               e3w_1d  (jk) =          za0      + za1        * TANH(       (zw-zkth ) / zacr  )      &
309                  &                             + za2        * TANH(       (zw-zkth2) / zacr2 )
310               e3t_1d  (jk) =          za0      + za1        * TANH(       (zt-zkth ) / zacr  )      &
311                  &                             + za2        * TANH(       (zt-zkth2) / zacr2 )
312            END DO
313         ENDIF
314         gdepw_1d(1) = 0._wp                    ! force first w-level to be exactly at zero
315      ENDIF
316
317      IF ( ln_isfcav ) THEN
318! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth)
319! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively
320         DO jk = 1, jpkm1
321            e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) 
322         END DO
323         e3t_1d(jpk) = e3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO
324
325         DO jk = 2, jpk
326            e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1) 
327         END DO
328         e3w_1d(1  ) = 2._wp * (gdept_1d(1) - gdepw_1d(1)) 
329      END IF
330
331!!gm BUG in s-coordinate this does not work!
332      ! deepest/shallowest W level Above/Below ~10m
333      zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d )                   ! ref. depth with tolerance (10% of minimum layer thickness)
334      nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
335      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m
336!!gm end bug
337
338      IF(lwp) THEN                        ! control print
339         WRITE(numout,*)
340         WRITE(numout,*) '              Reference z-coordinate depth and scale factors:'
341         WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
342         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk )
343         IF(lflush) CALL flush(numout)
344      ENDIF
345      DO jk = 1, jpk                      ! control positivity
346         IF( e3w_1d  (jk) <= 0._wp .OR. e3t_1d  (jk) <= 0._wp )   CALL ctl_stop( 'dom:zgr_z: e3w_1d or e3t_1d =< 0 '    )
347         IF( gdepw_1d(jk) <  0._wp .OR. gdept_1d(jk) <  0._wp )   CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' )
348      END DO
349      !
350      IF( nn_timing == 1 )  CALL timing_stop('zgr_z')
351      !
352   END SUBROUTINE zgr_z
353
354
355   SUBROUTINE zgr_bat
356      !!----------------------------------------------------------------------
357      !!                    ***  ROUTINE zgr_bat  ***
358      !!
359      !! ** Purpose :   set bathymetry both in levels and meters
360      !!
361      !! ** Method  :   read or define mbathy and bathy arrays
362      !!       * level bathymetry:
363      !!      The ocean basin geometry is given by a two-dimensional array,
364      !!      mbathy, which is defined as follow :
365      !!            mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
366      !!                              at t-point (ji,jj).
367      !!                            = 0  over the continental t-point.
368      !!      The array mbathy is checked to verified its consistency with
369      !!      model option. in particular:
370      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
371      !!                  along closed boundary.
372      !!            mbathy must be cyclic IF jperio=1.
373      !!            mbathy must be lower or equal to jpk-1.
374      !!            isolated ocean grid points are suppressed from mbathy
375      !!                  since they are only connected to remaining
376      !!                  ocean through vertical diffusion.
377      !!      ntopo=-1 :   rectangular channel or bassin with a bump
378      !!      ntopo= 0 :   flat rectangular channel or basin
379      !!      ntopo= 1 :   mbathy is read in 'bathy_level.nc' NetCDF file
380      !!                   bathy  is read in 'bathy_meter.nc' NetCDF file
381      !!
382      !! ** Action  : - mbathy: level bathymetry (in level index)
383      !!              - bathy : meter bathymetry (in meters)
384      !!----------------------------------------------------------------------
385      INTEGER  ::   ji, jj, jl, jk            ! dummy loop indices
386      INTEGER  ::   inum                      ! temporary logical unit
387      INTEGER  ::   ierror                    ! error flag
388      INTEGER  ::   ii_bump, ij_bump, ih      ! bump center position
389      INTEGER  ::   ii0, ii1, ij0, ij1, ik    ! local indices
390      REAL(wp) ::   r_bump , h_bump , h_oce   ! bump characteristics
391      REAL(wp) ::   zi, zj, zh, zhmin         ! local scalars
392      INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   idta   ! global domain integer data
393      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdta   ! global domain scalar data
394      !!----------------------------------------------------------------------
395      !
396      IF( nn_timing == 1 )  CALL timing_start('zgr_bat')
397      !
398      IF(lwp) WRITE(numout,*)
399      IF(lwp) WRITE(numout,*) '    zgr_bat : defines level and meter bathymetry'
400      IF(lwp) WRITE(numout,*) '    ~~~~~~~'
401      !                                               ! ================== !
402      IF( ntopo == 0 .OR. ntopo == -1 ) THEN          !   defined by hand  !
403         !                                            ! ================== !
404         !                                            ! global domain level and meter bathymetry (idta,zdta)
405         !
406         ALLOCATE( idta(jpidta,jpjdta), STAT=ierror )
407         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate idta array' )
408         ALLOCATE( zdta(jpidta,jpjdta), STAT=ierror )
409         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zdta array' )
410         !
411         IF( ntopo == 0 ) THEN                        ! flat basin
412            IF(lwp) WRITE(numout,*)
413            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin'
414            IF( rn_bathy > 0.01 ) THEN
415               IF(lwp) WRITE(numout,*) '         Depth = rn_bathy read in namelist'
416               zdta(:,:) = rn_bathy
417               IF( ln_sco ) THEN                                   ! s-coordinate (zsc       ): idta()=jpk
418                  idta(:,:) = jpkm1
419               ELSE                                                ! z-coordinate (zco or zps): step-like topography
420                  idta(:,:) = jpkm1
421                  DO jk = 1, jpkm1
422                     WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) )   idta(:,:) = jk
423                  END DO
424               ENDIF
425            ELSE
426               IF(lwp) WRITE(numout,*) '         Depth = depthw(jpkm1)'
427               idta(:,:) = jpkm1                            ! before last level
428               zdta(:,:) = gdepw_1d(jpk)                     ! last w-point depth
429               h_oce     = gdepw_1d(jpk)
430            ENDIF
431         ELSE                                         ! bump centered in the basin
432            IF(lwp) WRITE(numout,*)
433            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin with a bump'
434            ii_bump = jpidta / 2                           ! i-index of the bump center
435            ij_bump = jpjdta / 2                           ! j-index of the bump center
436            r_bump  = 50000._wp                            ! bump radius (meters)       
437            h_bump  =  2700._wp                            ! bump height (meters)
438            h_oce   = gdepw_1d(jpk)                        ! background ocean depth (meters)
439            IF(lwp) WRITE(numout,*) '            bump characteristics: '
440            IF(lwp) WRITE(numout,*) '               bump center (i,j)   = ', ii_bump, ii_bump
441            IF(lwp) WRITE(numout,*) '               bump height         = ', h_bump , ' meters'
442            IF(lwp) WRITE(numout,*) '               bump radius         = ', r_bump , ' index'
443            IF(lwp) WRITE(numout,*) '            background ocean depth = ', h_oce  , ' meters'
444            !                                       
445            DO jj = 1, jpjdta                              ! zdta :
446               DO ji = 1, jpidta
447                  zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
448                  zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
449                  zdta(ji,jj) = h_oce - h_bump * EXP( -( zi*zi + zj*zj ) )
450               END DO
451            END DO
452            !                                              ! idta :
453            IF( ln_sco ) THEN                                   ! s-coordinate (zsc       ): idta()=jpk
454               idta(:,:) = jpkm1
455            ELSE                                                ! z-coordinate (zco or zps): step-like topography
456               idta(:,:) = jpkm1
457               DO jk = 1, jpkm1
458                  WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) )   idta(:,:) = jk
459               END DO
460            ENDIF
461         ENDIF
462         !                                            ! set GLOBAL boundary conditions
463         !                                            ! Caution : idta on the global domain: use of jperio, not nperio
464         IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
465            idta( :    , 1    ) = -1                ;      zdta( :    , 1    ) = -1._wp
466            idta( :    ,jpjdta) =  0                ;      zdta( :    ,jpjdta) =  0._wp
467         ELSEIF( jperio == 2 ) THEN
468            idta( :    , 1    ) = idta( : ,  3  )   ;      zdta( :    , 1    ) = zdta( : ,  3  )
469            idta( :    ,jpjdta) = 0                 ;      zdta( :    ,jpjdta) =  0._wp
470            idta( 1    , :    ) = 0                 ;      zdta( 1    , :    ) =  0._wp
471            idta(jpidta, :    ) = 0                 ;      zdta(jpidta, :    ) =  0._wp
472         ELSE
473            ih = 0                                  ;      zh = 0._wp
474            IF( ln_sco )   ih = jpkm1               ;      IF( ln_sco )   zh = h_oce
475            idta( :    , 1    ) = ih                ;      zdta( :    , 1    ) =  zh
476            idta( :    ,jpjdta) = ih                ;      zdta( :    ,jpjdta) =  zh
477            idta( 1    , :    ) = ih                ;      zdta( 1    , :    ) =  zh
478            idta(jpidta, :    ) = ih                ;      zdta(jpidta, :    ) =  zh
479         ENDIF
480
481         !                                            ! local domain level and meter bathymetries (mbathy,bathy)
482         mbathy(:,:) = 0                                   ! set to zero extra halo points
483         bathy (:,:) = 0._wp                               ! (require for mpp case)
484         DO jj = 1, nlcj                                   ! interior values
485            DO ji = 1, nlci
486               mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
487               bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
488            END DO
489         END DO
490         risfdep(:,:)=0.e0
491         misfdep(:,:)=1
492         !
493         DEALLOCATE( idta, zdta )
494         !
495         !                                            ! ================ !
496      ELSEIF( ntopo == 1 ) THEN                       !   read in file   ! (over the local domain)
497         !                                            ! ================ !
498         !
499         IF( ln_zco )   THEN                          ! zco : read level bathymetry
500            CALL iom_open ( 'bathy_level.nc', inum ) 
501            CALL iom_get  ( inum, jpdom_data, 'Bathy_level', bathy )
502            CALL iom_close( inum )
503            mbathy(:,:) = INT( bathy(:,:) )
504            !                                                ! =====================
505            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
506               !                                             ! =====================
507               IF( nn_cla == 0 ) THEN
508                  ii0 = 140   ;   ii1 = 140                  ! Gibraltar Strait open
509                  ij0 = 102   ;   ij1 = 102                  ! (Thomson, Ocean Modelling, 1995)
510                  DO ji = mi0(ii0), mi1(ii1)
511                     DO jj = mj0(ij0), mj1(ij1)
512                        mbathy(ji,jj) = 15
513                     END DO
514                  END DO
515                  IF(lwp) WRITE(numout,*)
516                  IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
517                  !
518                  ii0 = 160   ;   ii1 = 160                  ! Bab el mandeb Strait open
519                  ij0 = 88    ;   ij1 = 88                   ! (Thomson, Ocean Modelling, 1995)
520                  DO ji = mi0(ii0), mi1(ii1)
521                     DO jj = mj0(ij0), mj1(ij1)
522                        mbathy(ji,jj) = 12
523                     END DO
524                  END DO
525                  IF(lwp) WRITE(numout,*)
526                  IF(lwp) WRITE(numout,*) '      orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
527               ENDIF
528               !
529            ENDIF
530            !
531         ENDIF
532         IF( ln_zps .OR. ln_sco )   THEN              ! zps or sco : read meter bathymetry
533            CALL iom_open ( 'bathy_meter.nc', inum ) 
534            IF ( ln_isfcav ) THEN
535               CALL iom_get  ( inum, jpdom_data, 'Bathymetry_isf', bathy, lrowattr=.false. )
536            ELSE
537               CALL iom_get  ( inum, jpdom_data, 'Bathymetry'    , bathy, lrowattr=ln_use_jattr  )
538            END IF
539            CALL iom_close( inum )
540            !                                               
541            risfdep(:,:)=0._wp         
542            misfdep(:,:)=1             
543            IF ( ln_isfcav ) THEN
544               CALL iom_open ( 'isf_draft_meter.nc', inum ) 
545               CALL iom_get  ( inum, jpdom_data, 'isf_draft', risfdep )
546               CALL iom_close( inum )
547               WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp
548            END IF
549            !       
550            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
551               !
552              IF( nn_cla == 0 ) THEN
553                 ii0 = 140   ;   ii1 = 140                   ! Gibraltar Strait open
554                 ij0 = 102   ;   ij1 = 102                   ! (Thomson, Ocean Modelling, 1995)
555                 DO ji = mi0(ii0), mi1(ii1)
556                    DO jj = mj0(ij0), mj1(ij1)
557                       bathy(ji,jj) = 284._wp
558                    END DO
559                 END DO
560                 IF(lwp) WRITE(numout,*)     
561                 IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
562                 !
563                 ii0 = 160   ;   ii1 = 160                   ! Bab el mandeb Strait open
564                 ij0 = 88    ;   ij1 = 88                    ! (Thomson, Ocean Modelling, 1995)
565                 DO ji = mi0(ii0), mi1(ii1)
566                    DO jj = mj0(ij0), mj1(ij1)
567                       bathy(ji,jj) = 137._wp
568                    END DO
569                 END DO
570                 IF(lwp) WRITE(numout,*)
571                 IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
572              ENDIF
573              !
574           ENDIF
575            !
576        ENDIF
577         !                                            ! =============== !
578      ELSE                                            !      error      !
579         !                                            ! =============== !
580         WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
581         CALL ctl_stop( '    zgr_bat : '//trim(ctmp1) )
582      ENDIF
583      !
584      IF( nn_closea == 0 )   CALL clo_bat( bathy, mbathy )    !==  NO closed seas or lakes  ==!
585      !                       
586      IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==!
587         ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin
588         IF ( ln_isfcav ) THEN
589            WHERE (bathy == risfdep)
590               bathy   = 0.0_wp ; risfdep = 0.0_wp
591            END WHERE
592         END IF
593         ! end patch
594         IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level
595         ELSE                          ;   ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 )  ! from a depth
596         ENDIF
597         zhmin = gdepw_1d(ik+1)                                                         ! minimum depth = ik+1 w-levels
598         WHERE( bathy(:,:) <= 0._wp )   ;   bathy(:,:) = 0._wp                         ! min=0     over the lands
599         ELSE WHERE                     ;   bathy(:,:) = MAX(  zhmin , bathy(:,:)  )   ! min=zhmin over the oceans
600         END WHERE
601         IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
602      ENDIF
603      !
604      IF(lwp .AND. lflush) CALL flush(numout)
605      !
606      IF( nn_timing == 1 )  CALL timing_stop('zgr_bat')
607      !
608   END SUBROUTINE zgr_bat
609
610
611   SUBROUTINE zgr_bat_zoom
612      !!----------------------------------------------------------------------
613      !!                    ***  ROUTINE zgr_bat_zoom  ***
614      !!
615      !! ** Purpose : - Close zoom domain boundary if necessary
616      !!              - Suppress Med Sea from ORCA R2 and R05 arctic zoom
617      !!
618      !! ** Method  :
619      !!
620      !! ** Action  : - update mbathy: level bathymetry (in level index)
621      !!----------------------------------------------------------------------
622      INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integers
623      !!----------------------------------------------------------------------
624      !
625      IF(lwp) WRITE(numout,*)
626      IF(lwp) WRITE(numout,*) '    zgr_bat_zoom : modify the level bathymetry for zoom domain'
627      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~'
628      !
629      ! Zoom domain
630      ! ===========
631      !
632      ! Forced closed boundary if required
633      IF( lzoom_s )   mbathy(  : , mj0(jpjzoom):mj1(jpjzoom) )      = 0
634      IF( lzoom_w )   mbathy(      mi0(jpizoom):mi1(jpizoom) , :  ) = 0
635      IF( lzoom_e )   mbathy(      mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , :  ) = 0
636      IF( lzoom_n )   mbathy(  : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) )      = 0
637      !
638      ! Configuration specific domain modifications
639      ! (here, ORCA arctic configuration: suppress Med Sea)
640      IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN
641         SELECT CASE ( jp_cfg )
642         !                                        ! =======================
643         CASE ( 2 )                               !  ORCA_R2 configuration
644            !                                     ! =======================
645            IF(lwp) WRITE(numout,*) '                   ORCA R2 arctic zoom: suppress the Med Sea'
646            ii0 = 141   ;   ii1 = 162      ! Sea box i,j indices
647            ij0 =  98   ;   ij1 = 110
648            !                                     ! =======================
649         CASE ( 05 )                              !  ORCA_R05 configuration
650            !                                     ! =======================
651            IF(lwp) WRITE(numout,*) '                   ORCA R05 arctic zoom: suppress the Med Sea'
652            ii0 = 563   ;   ii1 = 642      ! zero over the Med Sea boxe
653            ij0 = 314   ;   ij1 = 370 
654         END SELECT
655         !
656         mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0   ! zero over the Med Sea boxe
657         !
658      ENDIF
659      !
660      IF(lwp .AND. lflush) CALL flush(numout)
661      !
662   END SUBROUTINE zgr_bat_zoom
663
664
665   SUBROUTINE zgr_bat_ctl
666      !!----------------------------------------------------------------------
667      !!                    ***  ROUTINE zgr_bat_ctl  ***
668      !!
669      !! ** Purpose :   check the bathymetry in levels
670      !!
671      !! ** Method  :   The array mbathy is checked to verified its consistency
672      !!      with the model options. in particular:
673      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
674      !!                  along closed boundary.
675      !!            mbathy must be cyclic IF jperio=1.
676      !!            mbathy must be lower or equal to jpk-1.
677      !!            isolated ocean grid points are suppressed from mbathy
678      !!                  since they are only connected to remaining
679      !!                  ocean through vertical diffusion.
680      !!      C A U T I O N : mbathy will be modified during the initializa-
681      !!      tion phase to become the number of non-zero w-levels of a water
682      !!      column, with a minimum value of 1.
683      !!
684      !! ** Action  : - update mbathy: level bathymetry (in level index)
685      !!              - update bathy : meter bathymetry (in meters)
686      !!----------------------------------------------------------------------
687      !!
688      INTEGER ::   ji, jj, jl                    ! dummy loop indices
689      INTEGER ::   icompt, ibtest, ikmax         ! temporary integers
690      REAL(wp), POINTER, DIMENSION(:,:) ::  zbathy
691
692      !!----------------------------------------------------------------------
693      !
694      IF( nn_timing == 1 )  CALL timing_start('zgr_bat_ctl')
695      !
696      CALL wrk_alloc( jpi, jpj, zbathy )
697      !
698      IF(lwp) WRITE(numout,*)
699      IF(lwp) WRITE(numout,*) '    zgr_bat_ctl : check the bathymetry'
700      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
701      !                                          ! Suppress isolated ocean grid points
702      IF(lwp) WRITE(numout,*)
703      IF(lwp) WRITE(numout,*)'                   suppress isolated ocean grid points'
704      IF(lwp) WRITE(numout,*)'                   -----------------------------------'
705      icompt = 0
706      DO jl = 1, 2
707         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
708            mbathy( 1 ,:) = mbathy(jpim1,:)           ! local domain is cyclic east-west
709            mbathy(jpi,:) = mbathy(  2  ,:)
710         ENDIF
711         DO jj = 2, jpjm1
712            DO ji = 2, jpim1
713               ibtest = MAX(  mbathy(ji-1,jj), mbathy(ji+1,jj),   &
714                  &           mbathy(ji,jj-1), mbathy(ji,jj+1)  )
715               IF( ibtest < mbathy(ji,jj) ) THEN
716                  IF(lwp) WRITE(numout,*) ' the number of ocean level at ',   &
717                     &   'grid-point (i,j) =  ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
718                  mbathy(ji,jj) = ibtest
719                  icompt = icompt + 1
720               ENDIF
721            END DO
722         END DO
723      END DO
724      IF( lk_mpp )   CALL mpp_sum( icompt )
725      IF( icompt == 0 ) THEN
726         IF(lwp) WRITE(numout,*)'     no isolated ocean grid points'
727      ELSE
728         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed'
729      ENDIF
730      IF( lk_mpp ) THEN
731         zbathy(:,:) = FLOAT( mbathy(:,:) )
732         CALL lbc_lnk( zbathy, 'T', 1._wp )
733         mbathy(:,:) = INT( zbathy(:,:) )
734      ENDIF
735      !                                          ! East-west cyclic boundary conditions
736      IF( nperio == 0 ) THEN
737         IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio
738         IF( lk_mpp ) THEN
739            IF( nbondi == -1 .OR. nbondi == 2 ) THEN
740               IF( jperio /= 1 )   mbathy(1,:) = 0
741            ENDIF
742            IF( nbondi == 1 .OR. nbondi == 2 ) THEN
743               IF( jperio /= 1 )   mbathy(nlci,:) = 0
744            ENDIF
745         ELSE
746            IF( ln_zco .OR. ln_zps ) THEN
747               mbathy( 1 ,:) = 0
748               mbathy(jpi,:) = 0
749            ELSE
750               mbathy( 1 ,:) = jpkm1
751               mbathy(jpi,:) = jpkm1
752            ENDIF
753         ENDIF
754      ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio ==  6 ) THEN
755         IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio
756         mbathy( 1 ,:) = mbathy(jpim1,:)
757         mbathy(jpi,:) = mbathy(  2  ,:)
758      ELSEIF( nperio == 2 ) THEN
759         IF(lwp) WRITE(numout,*) '   equatorial boundary conditions on mbathy: nperio = ', nperio
760      ELSE
761         IF(lwp) WRITE(numout,*) '    e r r o r'
762         IF(lwp) WRITE(numout,*) '    parameter , nperio = ', nperio
763         !         STOP 'dom_mba'
764      ENDIF
765      !  Boundary condition on mbathy
766      IF( .NOT.lk_mpp ) THEN 
767!!gm     !!bug ???  think about it !
768         !   ... mono- or macro-tasking: T-point, >0, 2D array, no slab
769         zbathy(:,:) = FLOAT( mbathy(:,:) )
770         CALL lbc_lnk( zbathy, 'T', 1._wp )
771         mbathy(:,:) = INT( zbathy(:,:) )
772      ENDIF
773      ! Number of ocean level inferior or equal to jpkm1
774      ikmax = 0
775      DO jj = 1, jpj
776         DO ji = 1, jpi
777            ikmax = MAX( ikmax, mbathy(ji,jj) )
778         END DO
779      END DO
780!!gm  !!! test to do:   ikmax = MAX( mbathy(:,:) )   ???
781      IF( ikmax > jpkm1 ) THEN
782         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' >  jpk-1'
783         IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
784      ELSE IF( ikmax < jpkm1 ) THEN
785         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1' 
786         IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
787      ENDIF
788
789      IF( lwp .AND. nprint > 3 ) THEN      ! control print
790         WRITE(numout,*)
791         WRITE(numout,*) ' bathymetric field :   number of non-zero T-levels '
792         WRITE(numout,*) ' ------------------'
793         CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )
794         WRITE(numout,*)
795      ENDIF
796      !
797      CALL wrk_dealloc( jpi, jpj, zbathy )
798      !
799      IF(lwp .AND. lflush) CALL flush(numout)
800      !
801      IF( nn_timing == 1 )  CALL timing_stop('zgr_bat_ctl')
802      !
803   END SUBROUTINE zgr_bat_ctl
804
805
806   SUBROUTINE zgr_bot_level
807      !!----------------------------------------------------------------------
808      !!                    ***  ROUTINE zgr_bot_level  ***
809      !!
810      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
811      !!
812      !! ** Method  :   computes from mbathy with a minimum value of 1 over land
813      !!
814      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest
815      !!                                     ocean level at t-, u- & v-points
816      !!                                     (min value = 1 over land)
817      !!----------------------------------------------------------------------
818      !!
819      INTEGER ::   ji, jj   ! dummy loop indices
820      REAL(wp), POINTER, DIMENSION(:,:) ::  zmbk
821      !!----------------------------------------------------------------------
822      !
823      IF( nn_timing == 1 )  CALL timing_start('zgr_bot_level')
824      !
825      CALL wrk_alloc( jpi, jpj, zmbk )
826      !
827      IF(lwp) WRITE(numout,*)
828      IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
829      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
830      IF(lwp .AND. lflush) CALL flush(numout)
831      !
832      mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
833 
834      !                                     ! bottom k-index of W-level = mbkt+1
835      DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level
836         DO ji = 1, jpim1
837            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
838            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
839         END DO
840      END DO
841      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
842      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
843      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
844      !
845      CALL wrk_dealloc( jpi, jpj, zmbk )
846      !
847      IF( nn_timing == 1 )  CALL timing_stop('zgr_bot_level')
848      !
849   END SUBROUTINE zgr_bot_level
850
851      SUBROUTINE zgr_top_level
852      !!----------------------------------------------------------------------
853      !!                    ***  ROUTINE zgr_bot_level  ***
854      !!
855      !! ** Purpose :   defines the vertical index of ocean top (mik. arrays)
856      !!
857      !! ** Method  :   computes from misfdep with a minimum value of 1
858      !!
859      !! ** Action  :   mikt, miku, mikv :   vertical indices of the shallowest
860      !!                                     ocean level at t-, u- & v-points
861      !!                                     (min value = 1)
862      !!----------------------------------------------------------------------
863      !!
864      INTEGER ::   ji, jj   ! dummy loop indices
865      REAL(wp), POINTER, DIMENSION(:,:) ::  zmik
866      !!----------------------------------------------------------------------
867      !
868      IF( nn_timing == 1 )  CALL timing_start('zgr_top_level')
869      !
870      CALL wrk_alloc( jpi, jpj, zmik )
871      !
872      IF(lwp) WRITE(numout,*)
873      IF(lwp) WRITE(numout,*) '    zgr_top_level : ocean top k-index of T-, U-, V- and W-levels '
874      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
875      IF(lwp .AND. lflush) CALL flush(numout)
876      !
877      mikt(:,:) = MAX( misfdep(:,:) , 1 )    ! top k-index of T-level (=1)
878      !                                      ! top k-index of W-level (=mikt)
879      DO jj = 1, jpjm1                       ! top k-index of U- (U-) level
880         DO ji = 1, jpim1
881            miku(ji,jj) = MAX(  mikt(ji+1,jj  ) , mikt(ji,jj)  )
882            mikv(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj)  )
883            mikf(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj), mikt(ji+1,jj  ), mikt(ji+1,jj+1)  )
884         END DO
885      END DO
886
887      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
888      zmik(:,:) = REAL( miku(:,:), wp )   ;   CALL lbc_lnk(zmik,'U',1.)   ;   miku  (:,:) = MAX( INT( zmik(:,:) ), 1 )
889      zmik(:,:) = REAL( mikv(:,:), wp )   ;   CALL lbc_lnk(zmik,'V',1.)   ;   mikv  (:,:) = MAX( INT( zmik(:,:) ), 1 )
890      zmik(:,:) = REAL( mikf(:,:), wp )   ;   CALL lbc_lnk(zmik,'F',1.)   ;   mikf  (:,:) = MAX( INT( zmik(:,:) ), 1 )
891      !
892      CALL wrk_dealloc( jpi, jpj, zmik )
893      !
894      IF( nn_timing == 1 )  CALL timing_stop('zgr_top_level')
895      !
896   END SUBROUTINE zgr_top_level
897
898   SUBROUTINE zgr_zco
899      !!----------------------------------------------------------------------
900      !!                  ***  ROUTINE zgr_zco  ***
901      !!
902      !! ** Purpose :   define the z-coordinate system
903      !!
904      !! ** Method  :   set 3D coord. arrays to reference 1D array
905      !!----------------------------------------------------------------------
906      INTEGER  ::   jk
907      !!----------------------------------------------------------------------
908      !
909      IF( nn_timing == 1 )  CALL timing_start('zgr_zco')
910      !
911      DO jk = 1, jpk
912         gdept_0 (:,:,jk) = gdept_1d(jk)
913         gdepw_0 (:,:,jk) = gdepw_1d(jk)
914         gdep3w_0(:,:,jk) = gdepw_1d(jk)
915         e3t_0   (:,:,jk) = e3t_1d  (jk)
916         e3u_0   (:,:,jk) = e3t_1d  (jk)
917         e3v_0   (:,:,jk) = e3t_1d  (jk)
918         e3f_0   (:,:,jk) = e3t_1d  (jk)
919         e3w_0   (:,:,jk) = e3w_1d  (jk)
920         e3uw_0  (:,:,jk) = e3w_1d  (jk)
921         e3vw_0  (:,:,jk) = e3w_1d  (jk)
922      END DO
923      !
924      IF( nn_timing == 1 )  CALL timing_stop('zgr_zco')
925      !
926   END SUBROUTINE zgr_zco
927
928
929   SUBROUTINE zgr_zps
930      !!----------------------------------------------------------------------
931      !!                  ***  ROUTINE zgr_zps  ***
932      !!                     
933      !! ** Purpose :   the depth and vertical scale factor in partial step
934      !!      z-coordinate case
935      !!
936      !! ** Method  :   Partial steps : computes the 3D vertical scale factors
937      !!      of T-, U-, V-, W-, UW-, VW and F-points that are associated with
938      !!      a partial step representation of bottom topography.
939      !!
940      !!        The reference depth of model levels is defined from an analytical
941      !!      function the derivative of which gives the reference vertical
942      !!      scale factors.
943      !!        From  depth and scale factors reference, we compute there new value
944      !!      with partial steps  on 3d arrays ( i, j, k ).
945      !!
946      !!              w-level: gdepw_0(i,j,k)  = gdep(k)
947      !!                       e3w_0(i,j,k) = dk(gdep)(k)     = e3(i,j,k)
948      !!              t-level: gdept_0(i,j,k)  = gdep(k+0.5)
949      !!                       e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5)
950      !!
951      !!        With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
952      !!      we find the mbathy index of the depth at each grid point.
953      !!      This leads us to three cases:
954      !!
955      !!              - bathy = 0 => mbathy = 0
956      !!              - 1 < mbathy < jpkm1   
957      !!              - bathy > gdepw_0(jpk) => mbathy = jpkm1 
958      !!
959      !!        Then, for each case, we find the new depth at t- and w- levels
960      !!      and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
961      !!      and f-points.
962      !!
963      !!        This routine is given as an example, it must be modified
964      !!      following the user s desiderata. nevertheless, the output as
965      !!      well as the way to compute the model levels and scale factors
966      !!      must be respected in order to insure second order accuracy
967      !!      schemes.
968      !!
969      !!         c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives
970      !!         - - - - - - -   gdept_0, gdepw_0 and e3. are positives
971      !!     
972      !!  Reference :   Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
973      !!----------------------------------------------------------------------
974      !!
975      INTEGER  ::   ji, jj, jk       ! dummy loop indices
976      INTEGER  ::   ik, it, ikb, ikt ! temporary integers
977      LOGICAL  ::   ll_print         ! Allow  control print for debugging
978      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
979      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
980      REAL(wp) ::   zmax             ! Maximum depth
981      REAL(wp) ::   zdiff            ! temporary scalar
982      REAL(wp) ::   zrefdep          ! temporary scalar
983      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt
984      !!---------------------------------------------------------------------
985      !
986      IF( nn_timing == 1 )  CALL timing_start('zgr_zps')
987      !
988      CALL wrk_alloc( jpi, jpj, jpk, zprt )
989      !
990      IF(lwp) WRITE(numout,*)
991      IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps'
992      IF(lwp) WRITE(numout,*) '    ~~~~~~~ '
993      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used'
994
995      ll_print = .FALSE.                   ! Local variable for debugging
996     
997      IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth
998         WRITE(numout,*)
999         WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)'
1000         CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
1001      ENDIF
1002
1003
1004      ! bathymetry in level (from bathy_meter)
1005      ! ===================
1006      zmax = gdepw_1d(jpk) + e3t_1d(jpk)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) )
1007      bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat)
1008      WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0
1009      ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level
1010      END WHERE
1011
1012      ! Compute mbathy for ocean points (i.e. the number of ocean levels)
1013      ! find the number of ocean levels such that the last level thickness
1014      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where
1015      ! e3t_1d is the reference level thickness
1016      DO jk = jpkm1, 1, -1
1017         zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1018         WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1
1019      END DO
1020
1021      IF ( ln_isfcav ) CALL zgr_isf
1022
1023      ! Scale factors and depth at T- and W-points
1024      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate
1025         gdept_0(:,:,jk) = gdept_1d(jk)
1026         gdepw_0(:,:,jk) = gdepw_1d(jk)
1027         e3t_0  (:,:,jk) = e3t_1d  (jk)
1028         e3w_0  (:,:,jk) = e3w_1d  (jk)
1029      END DO
1030      !
1031      DO jj = 1, jpj
1032         DO ji = 1, jpi
1033            ik = mbathy(ji,jj)
1034            IF( ik > 0 ) THEN               ! ocean point only
1035               ! max ocean level case
1036               IF( ik == jpkm1 ) THEN
1037                  zdepwp = bathy(ji,jj)
1038                  ze3tp  = bathy(ji,jj) - gdepw_1d(ik)
1039                  ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
1040                  e3t_0(ji,jj,ik  ) = ze3tp
1041                  e3t_0(ji,jj,ik+1) = ze3tp
1042                  e3w_0(ji,jj,ik  ) = ze3wp
1043                  e3w_0(ji,jj,ik+1) = ze3tp
1044                  gdepw_0(ji,jj,ik+1) = zdepwp
1045                  gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp
1046                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
1047                  !
1048               ELSE                         ! standard case
1049                  IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
1050                  ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
1051                  ENDIF
1052!gm Bug?  check the gdepw_1d
1053                  !       ... on ik
1054                  gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   &
1055                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   &
1056                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) ))
1057                  e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   & 
1058                     &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ) 
1059                  e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   &
1060                     &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
1061                  !       ... on ik+1
1062                  e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
1063                  e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
1064                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
1065               ENDIF
1066            ENDIF
1067         END DO
1068      END DO
1069      !
1070      it = 0
1071      DO jj = 1, jpj
1072         DO ji = 1, jpi
1073            ik = mbathy(ji,jj)
1074            IF( ik > 0 ) THEN               ! ocean point only
1075               e3tp (ji,jj) = e3t_0(ji,jj,ik)
1076               e3wp (ji,jj) = e3w_0(ji,jj,ik)
1077               ! test
1078               zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  )
1079               IF( zdiff <= 0._wp .AND. lwp ) THEN
1080                  it = it + 1
1081                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj
1082                  WRITE(numout,*) ' bathy = ', bathy(ji,jj)
1083                  WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
1084                  WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  )
1085               ENDIF
1086            ENDIF
1087         END DO
1088      END DO
1089      !
1090      IF ( ln_isfcav ) THEN
1091      ! (ISF) Definition of e3t, u, v, w for ISF case
1092         DO jj = 1, jpj 
1093            DO ji = 1, jpi 
1094               ik = misfdep(ji,jj) 
1095               IF( ik > 1 ) THEN               ! ice shelf point only
1096                  IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik) 
1097                  gdepw_0(ji,jj,ik) = risfdep(ji,jj) 
1098!gm Bug?  check the gdepw_0
1099               !       ... on ik
1100                  gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   & 
1101                     &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   & 
1102                     &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      ) 
1103                  e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 
1104                  e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik)
1105
1106                  IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)
1107                     e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 
1108                  ENDIF 
1109               !       ... on ik / ik-1
1110                  e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 
1111                  e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
1112! The next line isn't required and doesn't affect results - included for consistency with bathymetry code
1113                  gdept_0(ji,jj,ik-1) = gdept_1d(ik-1)
1114               ENDIF
1115            END DO
1116         END DO 
1117      !
1118         it = 0 
1119         DO jj = 1, jpj 
1120            DO ji = 1, jpi 
1121               ik = misfdep(ji,jj) 
1122               IF( ik > 1 ) THEN               ! ice shelf point only
1123                  e3tp (ji,jj) = e3t_0(ji,jj,ik  ) 
1124                  e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 
1125               ! test
1126                  zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  ) 
1127                  IF( zdiff <= 0. .AND. lwp ) THEN 
1128                     it = it + 1 
1129                     WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
1130                     WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 
1131                     WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
1132                     WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj) 
1133                  ENDIF
1134               ENDIF
1135            END DO
1136         END DO
1137      END IF
1138      ! END (ISF)
1139
1140      ! Scale factors and depth at U-, V-, UW and VW-points
1141      DO jk = 1, jpk                        ! initialisation to z-scale factors
1142         e3u_0 (:,:,jk) = e3t_1d(jk)
1143         e3v_0 (:,:,jk) = e3t_1d(jk)
1144         e3uw_0(:,:,jk) = e3w_1d(jk)
1145         e3vw_0(:,:,jk) = e3w_1d(jk)
1146      END DO
1147      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
1148         DO jj = 1, jpjm1
1149            DO ji = 1, fs_jpim1   ! vector opt.
1150               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
1151               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
1152               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
1153               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
1154            END DO
1155         END DO
1156      END DO
1157      IF ( ln_isfcav ) THEN
1158      ! (ISF) define e3uw (adapted for 2 cells in the water column)
1159         DO jj = 2, jpjm1 
1160            DO ji = 2, fs_jpim1   ! vector opt.
1161               ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))
1162               ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))
1163               IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) &
1164                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) )
1165               ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))
1166               ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))
1167               IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) &
1168                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) )
1169            END DO
1170         END DO
1171      END IF
1172
1173      CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions
1174      CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp )
1175      !
1176      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
1177         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk)
1178         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk)
1179         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk)
1180         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk)
1181      END DO
1182     
1183      ! Scale factor at F-point
1184      DO jk = 1, jpk                        ! initialisation to z-scale factors
1185         e3f_0(:,:,jk) = e3t_1d(jk)
1186      END DO
1187      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors
1188         DO jj = 1, jpjm1
1189            DO ji = 1, fs_jpim1   ! vector opt.
1190               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
1191            END DO
1192         END DO
1193      END DO
1194      CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions
1195      !
1196      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
1197         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk)
1198      END DO
1199!!gm  bug ? :  must be a do loop with mj0,mj1
1200      !
1201      e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2
1202      e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 
1203      e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 
1204      e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 
1205      e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 
1206
1207      ! Control of the sign
1208      IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' )
1209      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' )
1210      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' )
1211      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' )
1212     
1213      ! Compute gdep3w_0 (vertical sum of e3w)
1214      IF ( ln_isfcav ) THEN ! if cavity
1215         WHERE (misfdep == 0) misfdep = 1
1216         DO jj = 1,jpj
1217            DO ji = 1,jpi
1218               gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)
1219               DO jk = 2, misfdep(ji,jj)
1220                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
1221               END DO
1222               IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))
1223               DO jk = misfdep(ji,jj) + 1, jpk
1224                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
1225               END DO
1226            END DO
1227         END DO
1228      ELSE ! no cavity
1229         gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
1230         DO jk = 2, jpk
1231            gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk)
1232         END DO
1233      END IF
1234      !                                               ! ================= !
1235      IF(lwp .AND. ll_print) THEN                     !   Control print   !
1236         !                                            ! ================= !
1237         DO jj = 1,jpj
1238            DO ji = 1, jpi
1239               ik = MAX( mbathy(ji,jj), 1 )
1240               zprt(ji,jj,1) = e3t_0   (ji,jj,ik)
1241               zprt(ji,jj,2) = e3w_0   (ji,jj,ik)
1242               zprt(ji,jj,3) = e3u_0   (ji,jj,ik)
1243               zprt(ji,jj,4) = e3v_0   (ji,jj,ik)
1244               zprt(ji,jj,5) = e3f_0   (ji,jj,ik)
1245               zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)
1246            END DO
1247         END DO
1248         WRITE(numout,*)
1249         WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1250         WRITE(numout,*)
1251         WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1252         WRITE(numout,*)
1253         WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1254         WRITE(numout,*)
1255         WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1256         WRITE(numout,*)
1257         WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1258         WRITE(numout,*)
1259         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1260         IF(lflush) CALL flush(numout)
1261      ENDIF 
1262      !
1263      CALL wrk_dealloc( jpi, jpj, jpk, zprt )
1264      !
1265      IF( nn_timing == 1 )  CALL timing_stop('zgr_zps')
1266      !
1267   END SUBROUTINE zgr_zps
1268
1269   SUBROUTINE zgr_isf
1270      !!----------------------------------------------------------------------
1271      !!                    ***  ROUTINE zgr_isf  ***
1272      !!   
1273      !! ** Purpose :   check the bathymetry in levels
1274      !!   
1275      !! ** Method  :   THe water column have to contained at least 2 cells
1276      !!                Bathymetry and isfdraft are modified (dig/close) to respect
1277      !!                this criterion.
1278      !!                 
1279      !!   
1280      !! ** Action  : - test compatibility between isfdraft and bathy
1281      !!              - bathy and isfdraft are modified
1282      !!----------------------------------------------------------------------
1283      !!   
1284      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
1285      INTEGER  ::   ik, it           ! temporary integers
1286      INTEGER  ::   id, jd, nprocd
1287      INTEGER  ::   icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1   ! (ISF)
1288      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
1289      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
1290      REAL(wp) ::   zmax, zmin       ! Maximum and minimum depth
1291      REAL(wp) ::   zdiff            ! temporary scalar
1292      REAL(wp) ::   zrefdep          ! temporary scalar
1293      REAL(wp) ::   zbathydiff, zrisfdepdiff  ! isf temporary scalar
1294      REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 2D workspace (ISH)
1295      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep         ! 2D workspace (ISH)
1296      !!---------------------------------------------------------------------
1297      !
1298      IF( nn_timing == 1 )  CALL timing_start('zgr_isf')
1299      !
1300      CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep)
1301      CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy )
1302
1303
1304      ! (ISF) compute misfdep
1305      WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1 
1306      ELSEWHERE                      ;                          misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level
1307      END WHERE 
1308
1309      ! Compute misfdep for ocean points (i.e. first wet level)
1310      ! find the first ocean level such that the first level thickness
1311      ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where
1312      ! e3t_0 is the reference level thickness
1313      DO jk = 2, jpkm1 
1314         zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
1315         WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1 
1316      END DO
1317      WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp)
1318         risfdep(:,:) = 0. ; misfdep(:,:) = 1
1319      END WHERE
1320 
1321! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation
1322      icompt = 0 
1323! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together
1324      DO jl = 1, 10     
1325         WHERE (bathy(:,:) .EQ. risfdep(:,:) )
1326            misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp
1327            mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp
1328         END WHERE
1329         WHERE (mbathy(:,:) <= 0) 
1330            misfdep(:,:) = 0; risfdep(:,:) = 0._wp 
1331            mbathy (:,:) = 0; bathy  (:,:) = 0._wp
1332         ENDWHERE
1333         IF( lk_mpp ) THEN
1334            zbathy(:,:) = FLOAT( misfdep(:,:) )
1335            CALL lbc_lnk( zbathy, 'T', 1. )
1336            misfdep(:,:) = INT( zbathy(:,:) )
1337            CALL lbc_lnk( risfdep, 'T', 1. )
1338            CALL lbc_lnk( bathy, 'T', 1. )
1339            zbathy(:,:) = FLOAT( mbathy(:,:) )
1340            CALL lbc_lnk( zbathy, 'T', 1. )
1341            mbathy(:,:) = INT( zbathy(:,:) )
1342         ENDIF
1343         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
1344            misfdep( 1 ,:) = misfdep(jpim1,:)           ! local domain is cyclic east-west
1345            misfdep(jpi,:) = misfdep(  2  ,:) 
1346         ENDIF
1347
1348         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
1349            mbathy( 1 ,:) = mbathy(jpim1,:)             ! local domain is cyclic east-west
1350            mbathy(jpi,:) = mbathy(  2  ,:)
1351         ENDIF
1352
1353         ! split last cell if possible (only where water column is 2 cell or less)
1354         DO jk = jpkm1, 1, -1
1355            zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1356            WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)
1357               mbathy(:,:) = jk
1358               bathy(:,:)  = zmax
1359            END WHERE
1360         END DO
1361 
1362         ! split top cell if possible (only where water column is 2 cell or less)
1363         DO jk = 2, jpkm1
1364            zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1365            WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy)
1366               misfdep(:,:) = jk
1367               risfdep(:,:) = zmax
1368            END WHERE
1369         END DO
1370
1371 
1372 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition
1373         DO jj = 1, jpj
1374            DO ji = 1, jpi
1375               ! find the minimum change option:
1376               ! test bathy
1377               IF (risfdep(ji,jj) .GT. 1) THEN
1378               zbathydiff =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) &
1379                 &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
1380               zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) &
1381                 &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
1382 
1383                  IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN
1384                     IF (zbathydiff .LE. zrisfdepdiff) THEN
1385                        bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )
1386                        mbathy(ji,jj)= mbathy(ji,jj) + 1
1387                     ELSE
1388                        risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )
1389                        misfdep(ji,jj) = misfdep(ji,jj) - 1
1390                     END IF
1391                  END IF
1392               END IF
1393            END DO
1394         END DO
1395 
1396          ! At least 2 levels for water thickness at T, U, and V point.
1397         DO jj = 1, jpj
1398            DO ji = 1, jpi
1399               ! find the minimum change option:
1400               ! test bathy
1401               IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
1402                  zbathydiff  =ABS(bathy(ji,jj)  - (gdepw_1d(mbathy (ji,jj)+1)&
1403                    & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
1404                  zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  )  &
1405                    & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
1406                  IF (zbathydiff .LE. zrisfdepdiff) THEN
1407                     mbathy(ji,jj) = mbathy(ji,jj) + 1
1408                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
1409                  ELSE
1410                     misfdep(ji,jj)= misfdep(ji,jj) - 1
1411                     risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )
1412                  END IF
1413               ENDIF
1414            END DO
1415         END DO
1416 
1417 ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1)
1418         DO jj = 1, jpjm1
1419            DO ji = 1, jpim1
1420               IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
1421                  zbathydiff =ABS(bathy(ji,jj    ) - (gdepw_1d(mbathy (ji,jj)+1)   &
1422                    &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat )))
1423                  zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1))  &
1424                    &  - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))
1425                  IF (zbathydiff .LE. zrisfdepdiff) THEN
1426                     mbathy(ji,jj) = mbathy(ji,jj) + 1
1427                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) &
1428                   &    + MIN( e3zps_min, e3t_1d(mbathy(ji,jj   )+1)*e3zps_rat )
1429                  ELSE
1430                     misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1
1431                     risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) &
1432                   &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )
1433                  END IF
1434               ENDIF
1435            END DO
1436         END DO
1437 
1438         IF( lk_mpp ) THEN
1439            zbathy(:,:) = FLOAT( misfdep(:,:) )
1440            CALL lbc_lnk( zbathy, 'T', 1. )
1441            misfdep(:,:) = INT( zbathy(:,:) )
1442            CALL lbc_lnk( risfdep, 'T', 1. )
1443            CALL lbc_lnk( bathy, 'T', 1. )
1444            zbathy(:,:) = FLOAT( mbathy(:,:) )
1445            CALL lbc_lnk( zbathy, 'T', 1. )
1446            mbathy(:,:) = INT( zbathy(:,:) )
1447         ENDIF
1448 ! point V misdep(ji,jj) EQ mbathy(ji,jj+1)
1449         DO jj = 1, jpjm1
1450            DO ji = 1, jpim1
1451               IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN
1452                  zbathydiff =ABS(  bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) &
1453                   &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))
1454                  zrisfdepdiff=ABS(risfdep(ji,jj  ) - (gdepw_1d(misfdep(ji,jj  )  )  &
1455                   &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat )))
1456                  IF (zbathydiff .LE. zrisfdepdiff) THEN
1457                     mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1
1458                     bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )
1459                  ELSE
1460                     misfdep(ji,jj)   = misfdep(ji,jj) - 1
1461                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat )
1462                  END IF
1463               ENDIF
1464            END DO
1465         END DO
1466 
1467 
1468         IF( lk_mpp ) THEN
1469            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
1470            CALL lbc_lnk( zbathy, 'T', 1. ) 
1471            misfdep(:,:) = INT( zbathy(:,:) ) 
1472            CALL lbc_lnk( risfdep, 'T', 1. ) 
1473            CALL lbc_lnk( bathy, 'T', 1. )
1474            zbathy(:,:) = FLOAT( mbathy(:,:) )
1475            CALL lbc_lnk( zbathy, 'T', 1. )
1476            mbathy(:,:) = INT( zbathy(:,:) )
1477         ENDIF 
1478 
1479 ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1)
1480         DO jj = 1, jpjm1
1481            DO ji = 1, jpim1
1482               IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
1483                  zbathydiff =ABS(  bathy(ji  ,jj) - (gdepw_1d(mbathy (ji,jj)+1) &
1484                    &   + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat )))
1485                  zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) &
1486                    &  - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))
1487                  IF (zbathydiff .LE. zrisfdepdiff) THEN
1488                     mbathy(ji,jj) = mbathy(ji,jj) + 1
1489                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
1490                  ELSE
1491                     misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1
1492                     risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )
1493                  END IF
1494               ENDIF
1495            ENDDO
1496         ENDDO
1497 
1498         IF( lk_mpp ) THEN
1499            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
1500            CALL lbc_lnk( zbathy, 'T', 1. ) 
1501            misfdep(:,:) = INT( zbathy(:,:) ) 
1502            CALL lbc_lnk( risfdep, 'T', 1. ) 
1503            CALL lbc_lnk( bathy, 'T', 1. )
1504            zbathy(:,:) = FLOAT( mbathy(:,:) )
1505            CALL lbc_lnk( zbathy, 'T', 1. )
1506            mbathy(:,:) = INT( zbathy(:,:) )
1507         ENDIF 
1508 
1509 ! point U misfdep(ji,jj) EQ bathy(ji,jj+1)
1510         DO jj = 1, jpjm1
1511            DO ji = 1, jpim1
1512               IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
1513                  zbathydiff =ABS(  bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) &
1514                      &   + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))
1515                  zrisfdepdiff=ABS(risfdep(ji  ,jj) - (gdepw_1d(misfdep(ji  ,jj)  )  &
1516                      &  - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat )))
1517                  IF (zbathydiff .LE. zrisfdepdiff) THEN
1518                     mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1
1519                     bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  )  &
1520                      &   + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat )
1521                  ELSE
1522                     misfdep(ji,jj)   = misfdep(ji  ,jj) - 1
1523                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) &
1524                      &   - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat )
1525                  END IF
1526               ENDIF
1527            ENDDO
1528         ENDDO
1529 
1530         IF( lk_mpp ) THEN
1531            zbathy(:,:) = FLOAT( misfdep(:,:) )
1532            CALL lbc_lnk( zbathy, 'T', 1. )
1533            misfdep(:,:) = INT( zbathy(:,:) )
1534            CALL lbc_lnk( risfdep, 'T', 1. )
1535            CALL lbc_lnk( bathy, 'T', 1. )
1536            zbathy(:,:) = FLOAT( mbathy(:,:) )
1537            CALL lbc_lnk( zbathy, 'T', 1. )
1538            mbathy(:,:) = INT( zbathy(:,:) )
1539         ENDIF
1540      END DO
1541      ! end dig bathy/ice shelf to be compatible
1542      ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness
1543      DO jl = 1,20
1544 
1545 ! remove single point "bay" on isf coast line in the ice shelf draft'
1546         DO jk = 2, jpk
1547            WHERE (misfdep==0) misfdep=jpk
1548            zmask=0
1549            WHERE (misfdep .LE. jk) zmask=1
1550            DO jj = 2, jpjm1
1551               DO ji = 2, jpim1
1552                  IF (misfdep(ji,jj) .EQ. jk) THEN
1553                     ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
1554                     IF (ibtest .LE. 1) THEN
1555                        risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1
1556                        IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) misfdep(ji,jj) = jpk
1557                     END IF
1558                  END IF
1559               END DO
1560            END DO
1561         END DO
1562         WHERE (misfdep==jpk)
1563             misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.
1564         END WHERE
1565         IF( lk_mpp ) THEN
1566            zbathy(:,:) = FLOAT( misfdep(:,:) )
1567            CALL lbc_lnk( zbathy, 'T', 1. )
1568            misfdep(:,:) = INT( zbathy(:,:) )
1569            CALL lbc_lnk( risfdep, 'T', 1. )
1570            CALL lbc_lnk( bathy, 'T', 1. )
1571            zbathy(:,:) = FLOAT( mbathy(:,:) )
1572            CALL lbc_lnk( zbathy, 'T', 1. )
1573            mbathy(:,:) = INT( zbathy(:,:) )
1574         ENDIF
1575 
1576 ! remove single point "bay" on bathy coast line beneath an ice shelf'
1577         DO jk = jpk,1,-1
1578            zmask=0
1579            WHERE (mbathy .GE. jk ) zmask=1
1580            DO jj = 2, jpjm1
1581               DO ji = 2, jpim1
1582                  IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) THEN
1583                     ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
1584                     IF (ibtest .LE. 1) THEN
1585                        bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1
1586                        IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 0
1587                     END IF
1588                  END IF
1589               END DO
1590            END DO
1591         END DO
1592         WHERE (mbathy==0)
1593             misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.
1594         END WHERE
1595         IF( lk_mpp ) THEN
1596            zbathy(:,:) = FLOAT( misfdep(:,:) )
1597            CALL lbc_lnk( zbathy, 'T', 1. )
1598            misfdep(:,:) = INT( zbathy(:,:) )
1599            CALL lbc_lnk( risfdep, 'T', 1. )
1600            CALL lbc_lnk( bathy, 'T', 1. )
1601            zbathy(:,:) = FLOAT( mbathy(:,:) )
1602            CALL lbc_lnk( zbathy, 'T', 1. )
1603            mbathy(:,:) = INT( zbathy(:,:) )
1604         ENDIF
1605 
1606 ! fill hole in ice shelf
1607         zmisfdep = misfdep
1608         zrisfdep = risfdep
1609         WHERE (zmisfdep .LE. 1) zmisfdep=jpk
1610         DO jj = 2, jpjm1
1611            DO ji = 2, jpim1
1612               ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  )
1613               ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1)
1614               IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj  ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj  ) - 1)
1615               IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj  ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj  ) - 1)
1616               IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji  ,jj-1) - 1)
1617               IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji  ,jj+1) - 1)
1618               ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
1619               IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN
1620                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp
1621               END IF
1622               IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN
1623                  misfdep(ji,jj) = ibtest
1624                  risfdep(ji,jj) = gdepw_1d(ibtest)
1625               ENDIF
1626            ENDDO
1627         ENDDO
1628 
1629         IF( lk_mpp ) THEN
1630            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
1631            CALL lbc_lnk( zbathy, 'T', 1. ) 
1632            misfdep(:,:) = INT( zbathy(:,:) ) 
1633            CALL lbc_lnk( risfdep, 'T', 1. ) 
1634            CALL lbc_lnk( bathy, 'T', 1. )
1635            zbathy(:,:) = FLOAT( mbathy(:,:) )
1636            CALL lbc_lnk( zbathy, 'T', 1. )
1637            mbathy(:,:) = INT( zbathy(:,:) )
1638         ENDIF 
1639 !
1640 !! fill hole in bathymetry
1641         zmbathy (:,:)=mbathy (:,:)
1642         DO jj = 2, jpjm1
1643            DO ji = 2, jpim1
1644               ibtestim1 = zmbathy(ji-1,jj  ) ; ibtestip1 = zmbathy(ji+1,jj  )
1645               ibtestjm1 = zmbathy(ji  ,jj-1) ; ibtestjp1 = zmbathy(ji  ,jj+1)
1646               IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj  ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj  ) + 1)
1647               IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj  ) ) ibtestip1 = 0
1648               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj-1) ) ibtestjm1 = 0
1649               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj+1) ) ibtestjp1 = 0
1650               ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
1651               IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN
1652                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ;
1653               END IF
1654               IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN
1655                  mbathy(ji,jj) = ibtest
1656                  bathy(ji,jj)  = gdepw_1d(ibtest+1) 
1657               ENDIF
1658            END DO
1659         END DO
1660         IF( lk_mpp ) THEN
1661            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
1662            CALL lbc_lnk( zbathy, 'T', 1. ) 
1663            misfdep(:,:) = INT( zbathy(:,:) ) 
1664            CALL lbc_lnk( risfdep, 'T', 1. ) 
1665            CALL lbc_lnk( bathy, 'T', 1. )
1666            zbathy(:,:) = FLOAT( mbathy(:,:) )
1667            CALL lbc_lnk( zbathy, 'T', 1. )
1668            mbathy(:,:) = INT( zbathy(:,:) )
1669         ENDIF 
1670 ! if not compatible after all check (ie U point water column less than 2 cells), mask U
1671         DO jj = 1, jpjm1
1672            DO ji = 1, jpim1
1673               IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN
1674                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ;
1675               END IF
1676            END DO
1677         END DO
1678         IF( lk_mpp ) THEN
1679            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
1680            CALL lbc_lnk( zbathy, 'T', 1. ) 
1681            misfdep(:,:) = INT( zbathy(:,:) ) 
1682            CALL lbc_lnk( risfdep, 'T', 1. ) 
1683            CALL lbc_lnk( bathy, 'T', 1. )
1684            zbathy(:,:) = FLOAT( mbathy(:,:) )
1685            CALL lbc_lnk( zbathy, 'T', 1. )
1686            mbathy(:,:) = INT( zbathy(:,:) )
1687         ENDIF 
1688 ! if not compatible after all check (ie U point water column less than 2 cells), mask U
1689         DO jj = 1, jpjm1
1690            DO ji = 1, jpim1
1691               IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN
1692                  mbathy(ji+1,jj)  = mbathy(ji+1,jj) - 1;   bathy(ji+1,jj)   = gdepw_1d(mbathy(ji+1,jj)+1) ;
1693               END IF
1694            END DO
1695         END DO
1696         IF( lk_mpp ) THEN
1697            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
1698            CALL lbc_lnk( zbathy, 'T', 1. ) 
1699            misfdep(:,:) = INT( zbathy(:,:) ) 
1700            CALL lbc_lnk( risfdep, 'T', 1. ) 
1701            CALL lbc_lnk( bathy, 'T', 1. )
1702            zbathy(:,:) = FLOAT( mbathy(:,:) )
1703            CALL lbc_lnk( zbathy, 'T', 1. )
1704            mbathy(:,:) = INT( zbathy(:,:) )
1705         ENDIF 
1706 ! if not compatible after all check (ie V point water column less than 2 cells), mask V
1707         DO jj = 1, jpjm1
1708            DO ji = 1, jpi
1709               IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN
1710                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ;
1711               END IF
1712            END DO
1713         END DO
1714         IF( lk_mpp ) THEN
1715            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
1716            CALL lbc_lnk( zbathy, 'T', 1. ) 
1717            misfdep(:,:) = INT( zbathy(:,:) ) 
1718            CALL lbc_lnk( risfdep, 'T', 1. ) 
1719            CALL lbc_lnk( bathy, 'T', 1. )
1720            zbathy(:,:) = FLOAT( mbathy(:,:) )
1721            CALL lbc_lnk( zbathy, 'T', 1. )
1722            mbathy(:,:) = INT( zbathy(:,:) )
1723         ENDIF 
1724 ! if not compatible after all check (ie V point water column less than 2 cells), mask V
1725         DO jj = 1, jpjm1
1726            DO ji = 1, jpi
1727               IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN
1728                  mbathy(ji,jj+1)  = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1)   = gdepw_1d(mbathy(ji,jj+1)+1) ;
1729               END IF
1730            END DO
1731         END DO
1732         IF( lk_mpp ) THEN
1733            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
1734            CALL lbc_lnk( zbathy, 'T', 1. ) 
1735            misfdep(:,:) = INT( zbathy(:,:) ) 
1736            CALL lbc_lnk( risfdep, 'T', 1. ) 
1737            CALL lbc_lnk( bathy, 'T', 1. )
1738            zbathy(:,:) = FLOAT( mbathy(:,:) )
1739            CALL lbc_lnk( zbathy, 'T', 1. )
1740            mbathy(:,:) = INT( zbathy(:,:) )
1741         ENDIF 
1742 ! if not compatible after all check, mask T
1743         DO jj = 1, jpj
1744            DO ji = 1, jpi
1745               IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN
1746                  misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj)  = 0 ; bathy(ji,jj)   = 0._wp ;
1747               END IF
1748            END DO
1749         END DO
1750 
1751         WHERE (mbathy(:,:) == 1)
1752            mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp
1753         END WHERE
1754      END DO 
1755! end check compatibility ice shelf/bathy
1756      ! remove very shallow ice shelf (less than ~ 10m if 75L)
1757      WHERE (misfdep(:,:) <= 5)
1758         misfdep = 1; risfdep = 0.0_wp;
1759      END WHERE
1760
1761      IF( icompt == 0 ) THEN
1762         IF(lwp) WRITE(numout,*)'     no points with ice shelf too close to bathymetry' 
1763      ELSE
1764         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry' 
1765      ENDIF
1766
1767      IF(lwp .AND. lflush) CALL flush(numout)
1768
1769      CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep )
1770      CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy )
1771
1772      IF( nn_timing == 1 )  CALL timing_stop('zgr_isf')
1773     
1774   END SUBROUTINE
1775
1776   SUBROUTINE zgr_sco
1777      !!----------------------------------------------------------------------
1778      !!                  ***  ROUTINE zgr_sco  ***
1779      !!                     
1780      !! ** Purpose :   define the s-coordinate system
1781      !!
1782      !! ** Method  :   s-coordinate
1783      !!         The depth of model levels is defined as the product of an
1784      !!      analytical function by the local bathymetry, while the vertical
1785      !!      scale factors are defined as the product of the first derivative
1786      !!      of the analytical function by the bathymetry.
1787      !!      (this solution save memory as depth and scale factors are not
1788      !!      3d fields)
1789      !!          - Read bathymetry (in meters) at t-point and compute the
1790      !!         bathymetry at u-, v-, and f-points.
1791      !!            hbatu = mi( hbatt )
1792      !!            hbatv = mj( hbatt )
1793      !!            hbatf = mi( mj( hbatt ) )
1794      !!          - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
1795      !!         function and its derivative given as function.
1796      !!            z_gsigt(k) = fssig (k    )
1797      !!            z_gsigw(k) = fssig (k-0.5)
1798      !!            z_esigt(k) = fsdsig(k    )
1799      !!            z_esigw(k) = fsdsig(k-0.5)
1800      !!      Three options for stretching are give, and they can be modified
1801      !!      following the users requirements. Nevertheless, the output as
1802      !!      well as the way to compute the model levels and scale factors
1803      !!      must be respected in order to insure second order accuracy
1804      !!      schemes.
1805      !!
1806      !!      The three methods for stretching available are:
1807      !!
1808      !!           s_sh94 (Song and Haidvogel 1994)
1809      !!                a sinh/tanh function that allows sigma and stretched sigma
1810      !!
1811      !!           s_sf12 (Siddorn and Furner 2012?)
1812      !!                allows the maintenance of fixed surface and or
1813      !!                bottom cell resolutions (cf. geopotential coordinates)
1814      !!                within an analytically derived stretched S-coordinate framework.
1815      !!
1816      !!          s_tanh  (Madec et al 1996)
1817      !!                a cosh/tanh function that gives stretched coordinates       
1818      !!
1819      !!----------------------------------------------------------------------
1820      !
1821      INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
1822      INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
1823      INTEGER  ::   ios                      ! Local integer output status for namelist read
1824      REAL(wp) ::   zrmax, ztaper   ! temporary scalars
1825      REAL(wp) ::   zrfact
1826      !
1827      REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
1828      REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
1829
1830      NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
1831                           rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
1832     !!----------------------------------------------------------------------
1833      !
1834      IF( nn_timing == 1 )  CALL timing_start('zgr_sco')
1835      !
1836      CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
1837      !
1838      REWIND( numnam_ref )              ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
1839      READ  ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901)
1840901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist', lwp )
1841
1842      REWIND( numnam_cfg )              ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters
1843      READ  ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 )
1844902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist', lwp )
1845      IF(lwm .AND. nprint > 2) WRITE ( numond, namzgr_sco )
1846
1847      IF(lwp) THEN                           ! control print
1848         WRITE(numout,*)
1849         WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
1850         WRITE(numout,*) '~~~~~~~~~~~'
1851         WRITE(numout,*) '   Namelist namzgr_sco'
1852         WRITE(numout,*) '     stretching coeffs '
1853         WRITE(numout,*) '        maximum depth of s-bottom surface (>0)       rn_sbot_max   = ',rn_sbot_max
1854         WRITE(numout,*) '        minimum depth of s-bottom surface (>0)       rn_sbot_min   = ',rn_sbot_min
1855         WRITE(numout,*) '        Critical depth                               rn_hc         = ',rn_hc
1856         WRITE(numout,*) '        maximum cut-off r-value allowed              rn_rmax       = ',rn_rmax
1857         WRITE(numout,*) '     Song and Haidvogel 1994 stretching              ln_s_sh94     = ',ln_s_sh94
1858         WRITE(numout,*) '        Song and Haidvogel 1994 stretching coefficients'
1859         WRITE(numout,*) '        surface control parameter (0<=rn_theta<=20)  rn_theta      = ',rn_theta
1860         WRITE(numout,*) '        bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ',rn_thetb
1861         WRITE(numout,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb
1862         WRITE(numout,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12
1863         WRITE(numout,*) '        switching to sigma (T) or Z (F) at H<Hc      ln_sigcrit    = ',ln_sigcrit
1864         WRITE(numout,*) '        Siddorn and Furner 2012 stretching coefficients'
1865         WRITE(numout,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha
1866         WRITE(numout,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold
1867         WRITE(numout,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs
1868         WRITE(numout,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a
1869         WRITE(numout,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b
1870         WRITE(numout,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b'
1871         IF(lflush) CALL flush(numout)
1872      ENDIF
1873
1874      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate
1875      hifu(:,:) = rn_sbot_min
1876      hifv(:,:) = rn_sbot_min
1877      hiff(:,:) = rn_sbot_min
1878
1879      !                                        ! set maximum ocean depth
1880      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1881
1882      DO jj = 1, jpj
1883         DO ji = 1, jpi
1884           IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1885         END DO
1886      END DO
1887      !                                        ! =============================
1888      !                                        ! Define the envelop bathymetry   (hbatt)
1889      !                                        ! =============================
1890      ! use r-value to create hybrid coordinates
1891      zenv(:,:) = bathy(:,:)
1892      !
1893      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
1894      DO jj = 1, jpj
1895         DO ji = 1, jpi
1896           IF( bathy(ji,jj) == 0._wp ) THEN
1897             iip1 = MIN( ji+1, jpi )
1898             ijp1 = MIN( jj+1, jpj )
1899             iim1 = MAX( ji-1, 1 )
1900             ijm1 = MAX( jj-1, 1 )
1901             IF( ( + bathy(iim1,ijp1) + bathy(ji,ijp1) + bathy(iip1,ijp1)  &
1902                &  + bathy(iim1,jj  )                  + bathy(iip1,jj  )  &
1903                &  + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1)  ) > 0._wp ) THEN
1904                zenv(ji,jj) = rn_sbot_min
1905             ENDIF
1906           ENDIF
1907         END DO
1908      END DO
1909      ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1910      CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1911      !
1912      ! smooth the bathymetry (if required)
1913      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
1914      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
1915      !
1916      jl = 0
1917      zrmax = 1._wp
1918      !   
1919      !     
1920      ! set scaling factor used in reducing vertical gradients
1921      zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
1922      !
1923      ! initialise temporary evelope depth arrays
1924      ztmpi1(:,:) = zenv(:,:)
1925      ztmpi2(:,:) = zenv(:,:)
1926      ztmpj1(:,:) = zenv(:,:)
1927      ztmpj2(:,:) = zenv(:,:)
1928      !
1929      ! initialise temporary r-value arrays
1930      zri(:,:) = 1._wp
1931      zrj(:,:) = 1._wp
1932      !                                                            ! ================ !
1933      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) !  Iterative loop  !
1934         !                                                         ! ================ !
1935         jl = jl + 1
1936         zrmax = 0._wp
1937         ! we set zrmax from previous r-values (zri and zrj) first
1938         ! if set after current r-value calculation (as previously)
1939         ! we could exit DO WHILE prematurely before checking r-value
1940         ! of current zenv
1941         DO jj = 1, nlcj
1942            DO ji = 1, nlci
1943               zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
1944            END DO
1945         END DO
1946         zri(:,:) = 0._wp
1947         zrj(:,:) = 0._wp
1948         DO jj = 1, nlcj
1949            DO ji = 1, nlci
1950               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
1951               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
1952               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
1953                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
1954               END IF
1955               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
1956                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
1957               END IF
1958               IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact
1959               IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact
1960               IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact
1961               IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact
1962            END DO
1963         END DO
1964         IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
1965         !
1966         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax
1967         !
1968         DO jj = 1, nlcj
1969            DO ji = 1, nlci
1970               zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
1971            END DO
1972         END DO
1973         ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1974         CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1975         !                                                  ! ================ !
1976      END DO                                                !     End loop     !
1977      !                                                     ! ================ !
1978      DO jj = 1, jpj
1979         DO ji = 1, jpi
1980            zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
1981         END DO
1982      END DO
1983      !
1984      ! Envelope bathymetry saved in hbatt
1985      hbatt(:,:) = zenv(:,:) 
1986      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
1987         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1988         DO jj = 1, jpj
1989            DO ji = 1, jpi
1990               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
1991               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
1992            END DO
1993         END DO
1994      ENDIF
1995      !
1996      IF(lwp) THEN                             ! Control print
1997         WRITE(numout,*)
1998         WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
1999         WRITE(numout,*)
2000         CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )
2001         IF( nprint > 3 )   THEN       
2002            WRITE(numout,*) ' bathy  MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
2003            WRITE(numout,*) ' hbatt  MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
2004         ENDIF
2005      ENDIF
2006
2007      !                                        ! ==============================
2008      !                                        !   hbatu, hbatv, hbatf fields
2009      !                                        ! ==============================
2010      IF(lwp) THEN
2011         WRITE(numout,*)
2012         WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
2013      ENDIF
2014      hbatu(:,:) = rn_sbot_min
2015      hbatv(:,:) = rn_sbot_min
2016      hbatf(:,:) = rn_sbot_min
2017      DO jj = 1, jpjm1
2018        DO ji = 1, jpim1   ! NO vector opt.
2019           hbatu(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
2020           hbatv(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
2021           hbatf(ji,jj) = 0.25_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
2022              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
2023        END DO
2024      END DO
2025      !
2026      ! Apply lateral boundary condition
2027!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
2028      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk( hbatu, 'U', 1._wp )
2029      DO jj = 1, jpj
2030         DO ji = 1, jpi
2031            IF( hbatu(ji,jj) == 0._wp ) THEN
2032               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
2033               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
2034            ENDIF
2035         END DO
2036      END DO
2037      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk( hbatv, 'V', 1._wp )
2038      DO jj = 1, jpj
2039         DO ji = 1, jpi
2040            IF( hbatv(ji,jj) == 0._wp ) THEN
2041               IF( zhbat(ji,jj) == 0._wp )   hbatv(ji,jj) = rn_sbot_min
2042               IF( zhbat(ji,jj) /= 0._wp )   hbatv(ji,jj) = zhbat(ji,jj)
2043            ENDIF
2044         END DO
2045      END DO
2046      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk( hbatf, 'F', 1._wp )
2047      DO jj = 1, jpj
2048         DO ji = 1, jpi
2049            IF( hbatf(ji,jj) == 0._wp ) THEN
2050               IF( zhbat(ji,jj) == 0._wp )   hbatf(ji,jj) = rn_sbot_min
2051               IF( zhbat(ji,jj) /= 0._wp )   hbatf(ji,jj) = zhbat(ji,jj)
2052            ENDIF
2053         END DO
2054      END DO
2055
2056!!bug:  key_helsinki a verifer
2057      hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
2058      hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
2059      hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
2060      hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
2061
2062      IF( nprint > 3 .AND. lwp )   THEN
2063         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
2064            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
2065         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
2066            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
2067         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
2068            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
2069         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
2070            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
2071      ENDIF
2072!! helsinki
2073
2074      !                                            ! =======================
2075      !                                            !   s-ccordinate fields     (gdep., e3.)
2076      !                                            ! =======================
2077      !
2078      ! non-dimensional "sigma" for model level depth at w- and t-levels
2079
2080
2081!========================================================================
2082! Song and Haidvogel  1994 (ln_s_sh94=T)
2083! Siddorn and Furner 2012 (ln_sf12=T)
2084! or  tanh function       (both false)                   
2085!========================================================================
2086      IF      ( ln_s_sh94 ) THEN
2087                           CALL s_sh94()
2088      ELSE IF ( ln_s_sf12 ) THEN
2089                           CALL s_sf12()
2090      ELSE                 
2091                           CALL s_tanh()
2092      ENDIF
2093
2094      CALL lbc_lnk( e3t_0 , 'T', 1._wp )
2095      CALL lbc_lnk( e3u_0 , 'U', 1._wp )
2096      CALL lbc_lnk( e3v_0 , 'V', 1._wp )
2097      CALL lbc_lnk( e3f_0 , 'F', 1._wp )
2098      CALL lbc_lnk( e3w_0 , 'W', 1._wp )
2099      CALL lbc_lnk( e3uw_0, 'U', 1._wp )
2100      CALL lbc_lnk( e3vw_0, 'V', 1._wp )
2101
2102      fsdepw(:,:,:) = gdepw_0 (:,:,:)
2103      fsde3w(:,:,:) = gdep3w_0(:,:,:)
2104      !
2105      where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0
2106      where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0
2107      where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0
2108      where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0
2109      where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0
2110      where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0
2111      where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0
2112
2113#if defined key_agrif
2114      ! Ensure meaningful vertical scale factors in ghost lines/columns
2115      IF( .NOT. Agrif_Root() ) THEN
2116         
2117         IF((nbondi == -1).OR.(nbondi == 2)) THEN
2118            e3u_0(1,:,:) = e3u_0(2,:,:)
2119         ENDIF
2120         !
2121         IF((nbondi ==  1).OR.(nbondi == 2)) THEN
2122            e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:)
2123         ENDIF
2124         !
2125         IF((nbondj == -1).OR.(nbondj == 2)) THEN
2126            e3v_0(:,1,:) = e3v_0(:,2,:)
2127         ENDIF
2128         !
2129         IF((nbondj ==  1).OR.(nbondj == 2)) THEN
2130            e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:)
2131         ENDIF
2132         !
2133      ENDIF
2134#endif
2135
2136      fsdept(:,:,:) = gdept_0 (:,:,:)
2137      fsdepw(:,:,:) = gdepw_0 (:,:,:)
2138      fsde3w(:,:,:) = gdep3w_0(:,:,:)
2139      fse3t (:,:,:) = e3t_0   (:,:,:)
2140      fse3u (:,:,:) = e3u_0   (:,:,:)
2141      fse3v (:,:,:) = e3v_0   (:,:,:)
2142      fse3f (:,:,:) = e3f_0   (:,:,:)
2143      fse3w (:,:,:) = e3w_0   (:,:,:)
2144      fse3uw(:,:,:) = e3uw_0  (:,:,:)
2145      fse3vw(:,:,:) = e3vw_0  (:,:,:)
2146!!
2147      ! HYBRID :
2148      DO jj = 1, jpj
2149         DO ji = 1, jpi
2150            DO jk = 1, jpkm1
2151               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
2152            END DO
2153            IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0
2154         END DO
2155      END DO
2156      IF( nprint > 3 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
2157         &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
2158
2159      IF( nprint > 3  .AND. lwp )   THEN         ! min max values over the local domain
2160         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) )
2161         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
2162            &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w '  , MINVAL( gdep3w_0(:,:,:) )
2163         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0   (:,:,:) ),   &
2164            &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0   (:,:,:) ),   &
2165            &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0  (:,:,:) ),   &
2166            &                          ' w ', MINVAL( e3w_0  (:,:,:) )
2167
2168         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
2169            &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w '  , MAXVAL( gdep3w_0(:,:,:) )
2170         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0   (:,:,:) ),   &
2171            &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0   (:,:,:) ),   &
2172            &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0  (:,:,:) ),   &
2173            &                          ' w ', MAXVAL( e3w_0  (:,:,:) )
2174      ENDIF
2175      !  END DO
2176      IF(lwp) THEN                                  ! selected vertical profiles
2177         WRITE(numout,*)
2178         WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
2179         WRITE(numout,*) ' ~~~~~~  --------------------'
2180         WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2181         WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     &
2182            &                                 e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
2183         DO jj = mj0(20), mj1(20)
2184            DO ji = mi0(20), mi1(20)
2185               WRITE(numout,*)
2186               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
2187               WRITE(numout,*) ' ~~~~~~  --------------------'
2188               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2189               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
2190                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2191            END DO
2192         END DO
2193         DO jj = mj0(74), mj1(74)
2194            DO ji = mi0(100), mi1(100)
2195               WRITE(numout,*)
2196               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
2197               WRITE(numout,*) ' ~~~~~~  --------------------'
2198               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2199               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
2200                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2201            END DO
2202         END DO
2203      ENDIF
2204
2205!================================================================================
2206! check the coordinate makes sense
2207!================================================================================
2208      DO ji = 1, jpi
2209         DO jj = 1, jpj
2210
2211            IF( hbatt(ji,jj) > 0._wp) THEN
2212               DO jk = 1, mbathy(ji,jj)
2213                 ! check coordinate is monotonically increasing
2214                 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
2215                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
2216                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
2217                    WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
2218                    WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
2219                    CALL ctl_stop( ctmp1 )
2220                 ENDIF
2221                 ! and check it has never gone negative
2222                 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
2223                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
2224                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk
2225                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2226                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2227                    CALL ctl_stop( ctmp1 )
2228                 ENDIF
2229                 ! and check it never exceeds the total depth
2230                 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
2231                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
2232                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
2233                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2234                    CALL ctl_stop( ctmp1 )
2235                 ENDIF
2236               END DO
2237
2238               DO jk = 1, mbathy(ji,jj)-1
2239                 ! and check it never exceeds the total depth
2240                IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
2241                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
2242                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
2243                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2244                    CALL ctl_stop( ctmp1 )
2245                 ENDIF
2246               END DO
2247
2248            ENDIF
2249
2250         END DO
2251      END DO
2252      !
2253      IF(lwp .AND. lflush) CALL flush(numout)
2254      !
2255      CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
2256      !
2257      IF( nn_timing == 1 )  CALL timing_stop('zgr_sco')
2258      !
2259   END SUBROUTINE zgr_sco
2260
2261!!======================================================================
2262   SUBROUTINE s_sh94()
2263
2264      !!----------------------------------------------------------------------
2265      !!                  ***  ROUTINE s_sh94  ***
2266      !!                     
2267      !! ** Purpose :   stretch the s-coordinate system
2268      !!
2269      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994
2270      !!                mixed S/sigma coordinate
2271      !!
2272      !! Reference : Song and Haidvogel 1994.
2273      !!----------------------------------------------------------------------
2274      !
2275      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2276      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
2277      !
2278      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2279      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
2280
2281      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2282      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2283
2284      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
2285      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
2286      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
2287      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
2288
2289      DO ji = 1, jpi
2290         DO jj = 1, jpj
2291
2292            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
2293               DO jk = 1, jpk
2294                  z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
2295                  z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
2296               END DO
2297            ELSE ! shallow water, uniform sigma
2298               DO jk = 1, jpk
2299                  z_gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
2300                  z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
2301                  END DO
2302            ENDIF
2303            !
2304            DO jk = 1, jpkm1
2305               z_esigt3(ji,jj,jk  ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2306               z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2307            END DO
2308            z_esigw3(ji,jj,1  ) = 2._wp * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) )
2309            z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
2310            !
2311            ! Coefficients for vertical depth as the sum of e3w scale factors
2312            z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
2313            DO jk = 2, jpk
2314               z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2315            END DO
2316            !
2317            DO jk = 1, jpk
2318               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2319               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2320               gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
2321               gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
2322               gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
2323            END DO
2324           !
2325         END DO   ! for all jj's
2326      END DO    ! for all ji's
2327
2328      DO ji = 1, jpim1
2329         DO jj = 1, jpjm1
2330            DO jk = 1, jpk
2331               z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
2332                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2333               z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
2334                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2335               z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     &
2336                  &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   &
2337                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2338               z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
2339                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2340               z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
2341                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2342               !
2343               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2344               e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2345               e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2346               e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2347               !
2348               e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2349               e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2350               e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2351            END DO
2352        END DO
2353      END DO
2354
2355      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2356      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2357
2358   END SUBROUTINE s_sh94
2359
2360   SUBROUTINE s_sf12
2361
2362      !!----------------------------------------------------------------------
2363      !!                  ***  ROUTINE s_sf12 ***
2364      !!                     
2365      !! ** Purpose :   stretch the s-coordinate system
2366      !!
2367      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
2368      !!                mixed S/sigma/Z coordinate
2369      !!
2370      !!                This method allows the maintenance of fixed surface and or
2371      !!                bottom cell resolutions (cf. geopotential coordinates)
2372      !!                within an analytically derived stretched S-coordinate framework.
2373      !!
2374      !!
2375      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
2376      !!----------------------------------------------------------------------
2377      !
2378      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2379      REAL(wp) ::   zsmth               ! smoothing around critical depth
2380      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
2381      !
2382      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2383      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
2384
2385      !
2386      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2387      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2388
2389      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
2390      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
2391      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
2392      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
2393
2394      DO ji = 1, jpi
2395         DO jj = 1, jpj
2396
2397          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
2398             
2399              zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
2400                                                     ! could be changed by users but care must be taken to do so carefully
2401              zzb = 1.0_wp-(zzb/hbatt(ji,jj))
2402           
2403              zzs = rn_zs / hbatt(ji,jj) 
2404             
2405              IF (rn_efold /= 0.0_wp) THEN
2406                zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
2407              ELSE
2408                zsmth = 1.0_wp 
2409              ENDIF
2410               
2411              DO jk = 1, jpk
2412                z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
2413                z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
2414              ENDDO
2415              z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth  )
2416              z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth  )
2417 
2418          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
2419
2420            DO jk = 1, jpk
2421              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
2422              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
2423            END DO
2424
2425          ELSE  ! shallow water, z coordinates
2426
2427            DO jk = 1, jpk
2428              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2429              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2430            END DO
2431
2432          ENDIF
2433
2434          DO jk = 1, jpkm1
2435             z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2436             z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2437          END DO
2438          z_esigw3(ji,jj,1  ) = 2.0_wp * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ))
2439          z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
2440
2441          ! Coefficients for vertical depth as the sum of e3w scale factors
2442          z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
2443          DO jk = 2, jpk
2444             z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2445          END DO
2446
2447          DO jk = 1, jpk
2448             gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
2449             gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
2450             gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
2451          END DO
2452
2453        ENDDO   ! for all jj's
2454      ENDDO    ! for all ji's
2455
2456      DO ji=1,jpi-1
2457        DO jj=1,jpj-1
2458
2459          DO jk = 1, jpk
2460                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / &
2461                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2462                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / &
2463                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2464                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  &
2465                                      hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / &
2466                                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2467                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / &
2468                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2469                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / &
2470                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2471
2472             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
2473             e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
2474             e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
2475             e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
2476             !
2477             e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
2478             e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
2479             e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
2480          END DO
2481
2482        ENDDO
2483      ENDDO
2484      !
2485      CALL lbc_lnk(e3t_0 ,'T',1.) ; CALL lbc_lnk(e3u_0 ,'T',1.)
2486      CALL lbc_lnk(e3v_0 ,'T',1.) ; CALL lbc_lnk(e3f_0 ,'T',1.)
2487      CALL lbc_lnk(e3w_0 ,'T',1.)
2488      CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.)
2489      !
2490      !                                               ! =============
2491
2492      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2493      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2494
2495   END SUBROUTINE s_sf12
2496
2497   SUBROUTINE s_tanh()
2498
2499      !!----------------------------------------------------------------------
2500      !!                  ***  ROUTINE s_tanh***
2501      !!                     
2502      !! ** Purpose :   stretch the s-coordinate system
2503      !!
2504      !! ** Method  :   s-coordinate stretch
2505      !!
2506      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
2507      !!----------------------------------------------------------------------
2508
2509      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2510      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
2511
2512      REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
2513      REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw
2514
2515      CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      )
2516      CALL wrk_alloc( jpk, z_esigt, z_esigw                                               )
2517
2518      z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp   ;   z_gsi3w  = 0._wp
2519      z_esigt  = 0._wp   ;   z_esigw  = 0._wp 
2520
2521      DO jk = 1, jpk
2522        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
2523        z_gsigt(jk) = -fssig( REAL(jk,wp)        )
2524      END DO
2525      IF( nprint > 3 .AND. lwp )   WRITE(numout,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk)
2526      !
2527      ! Coefficients for vertical scale factors at w-, t- levels
2528!!gm bug :  define it from analytical function, not like juste bellow....
2529!!gm        or betteroffer the 2 possibilities....
2530      DO jk = 1, jpkm1
2531         z_esigt(jk  ) = z_gsigw(jk+1) - z_gsigw(jk)
2532         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
2533      END DO
2534      z_esigw( 1 ) = 2._wp * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
2535      z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
2536      !
2537      ! Coefficients for vertical depth as the sum of e3w scale factors
2538      z_gsi3w(1) = 0.5_wp * z_esigw(1)
2539      DO jk = 2, jpk
2540         z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
2541      END DO
2542!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
2543      DO jk = 1, jpk
2544         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2545         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2546         gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
2547         gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
2548         gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )
2549      END DO
2550!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
2551      DO jj = 1, jpj
2552         DO ji = 1, jpi
2553            DO jk = 1, jpk
2554              e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2555              e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2556              e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2557              e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
2558              !
2559              e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2560              e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2561              e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2562            END DO
2563         END DO
2564      END DO
2565
2566      IF(lwp .AND. lflush) CALL flush(numout)
2567
2568      CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      )
2569      CALL wrk_dealloc( jpk, z_esigt, z_esigw                                               )
2570
2571   END SUBROUTINE s_tanh
2572
2573   FUNCTION fssig( pk ) RESULT( pf )
2574      !!----------------------------------------------------------------------
2575      !!                 ***  ROUTINE fssig ***
2576      !!       
2577      !! ** Purpose :   provide the analytical function in s-coordinate
2578      !!         
2579      !! ** Method  :   the function provide the non-dimensional position of
2580      !!                T and W (i.e. between 0 and 1)
2581      !!                T-points at integer values (between 1 and jpk)
2582      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2583      !!----------------------------------------------------------------------
2584      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
2585      REAL(wp)             ::   pf   ! sigma value
2586      !!----------------------------------------------------------------------
2587      !
2588      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   &
2589         &     - TANH( rn_thetb * rn_theta                                )  )   &
2590         & * (   COSH( rn_theta                           )                      &
2591         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              &
2592         & / ( 2._wp * SINH( rn_theta ) )
2593      !
2594   END FUNCTION fssig
2595
2596
2597   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
2598      !!----------------------------------------------------------------------
2599      !!                 ***  ROUTINE fssig1 ***
2600      !!
2601      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
2602      !!
2603      !! ** Method  :   the function provides the non-dimensional position of
2604      !!                T and W (i.e. between 0 and 1)
2605      !!                T-points at integer values (between 1 and jpk)
2606      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2607      !!----------------------------------------------------------------------
2608      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
2609      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
2610      REAL(wp)             ::   pf1   ! sigma value
2611      !!----------------------------------------------------------------------
2612      !
2613      IF ( rn_theta == 0 ) then      ! uniform sigma
2614         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
2615      ELSE                        ! stretched sigma
2616         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
2617            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
2618            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
2619      ENDIF
2620      !
2621   END FUNCTION fssig1
2622
2623
2624   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
2625      !!----------------------------------------------------------------------
2626      !!                 ***  ROUTINE fgamma  ***
2627      !!
2628      !! ** Purpose :   provide analytical function for the s-coordinate
2629      !!
2630      !! ** Method  :   the function provides the non-dimensional position of
2631      !!                T and W (i.e. between 0 and 1)
2632      !!                T-points at integer values (between 1 and jpk)
2633      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2634      !!
2635      !!                This method allows the maintenance of fixed surface and or
2636      !!                bottom cell resolutions (cf. geopotential coordinates)
2637      !!                within an analytically derived stretched S-coordinate framework.
2638      !!
2639      !! Reference  :   Siddorn and Furner, in prep
2640      !!----------------------------------------------------------------------
2641      REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate
2642      REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate
2643      REAL(wp), INTENT(in   ) ::   pzb           ! Bottom box depth
2644      REAL(wp), INTENT(in   ) ::   pzs           ! surface box depth
2645      REAL(wp), INTENT(in   ) ::   psmth       ! Smoothing parameter
2646      REAL(wp)                ::   za1,za2,za3    ! local variables
2647      REAL(wp)                ::   zn1,zn2        ! local variables
2648      REAL(wp)                ::   za,zb,zx       ! local variables
2649      integer                 ::   jk
2650      !!----------------------------------------------------------------------
2651      !
2652
2653      zn1  =  1./(jpk-1.)
2654      zn2  =  1. -  zn1
2655
2656      za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 
2657      za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp)
2658      za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
2659     
2660      za = pzb - za3*(pzs-za1)-za2
2661      za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) )
2662      zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
2663      zx = 1.0_wp-za/2.0_wp-zb
2664 
2665      DO jk = 1, jpk
2666        p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp +  &
2667                    & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &
2668                    &      (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
2669        p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth)
2670      ENDDO 
2671
2672      !
2673   END FUNCTION fgamma
2674
2675   !!======================================================================
2676END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.