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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 4468

Last change on this file since 4468 was 4468, checked in by trackstand2, 10 years ago

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