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

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

source: branches/UKMO/dev_3.6_FVPS/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 8476

Last change on this file since 8476 was 8476, checked in by davestorkey, 7 years ago

UKMO/dev_3.6_FVPS : implement FVPS HPG scheme.

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