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

source: trunk/NEMO/OPA_SRC/DOM/domzgr.F90 @ 900

Last change on this file since 900 was 900, checked in by rblod, 16 years ago

Update 1D configuration according to SBC and LIM3, see ticket #117

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 50.5 KB
Line 
1MODULE domzgr
2   !!==============================================================================
3   !!                       ***  MODULE domzgr   ***
4   !! Ocean initialization : domain initialization
5   !!==============================================================================
6
7   !!----------------------------------------------------------------------
8   !!   dom_zgr     : defined the ocean vertical coordinate system
9   !!       zgr_bat      : bathymetry fields (levels and meters)
10   !!       zgr_bat_zoom : modify the bathymetry field if zoom domain
11   !!       zgr_bat_ctl  : check the bathymetry files
12   !!       zgr_z        : reference z-coordinate
13   !!       zgr_zco      : z-coordinate
14   !!       zgr_zps      : z-coordinate with partial steps
15   !!       zgr_sco      : s-coordinate
16   !!---------------------------------------------------------------------
17   !! * Modules used
18   USE oce             ! ocean dynamics and tracers
19   USE dom_oce         ! ocean space and time domain
20   USE in_out_manager  ! I/O manager
21   USE lib_mpp         ! distributed memory computing library
22   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
23   USE closea
24   USE solisl
25   USE c1d             ! 1D configuration
26
27   IMPLICIT NONE
28   PRIVATE
29
30   !! * Routine accessibility
31   PUBLIC dom_zgr        ! called by dom_init.F90
32
33   !! * Module variables
34      REAL(wp) ::          & !!: Namelist nam_zgr_sco
35      sbot_min =  300.  ,  &  !: minimum depth of s-bottom surface (>0) (m)
36      sbot_max = 5250.  ,  &  !: maximum depth of s-bottom surface (= ocean depth) (>0) (m)
37      theta    =    6.0 ,  &  !: surface control parameter (0<=theta<=20)
38      thetb    =    0.75,  &  !: bottom control parameter  (0<=thetb<= 1)
39      r_max    =    0.15      !: maximum cut-off r-value allowed (0<r_max<1)
40
41
42 
43   !! * Substitutions
44#  include "domzgr_substitute.h90"
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !!   OPA 9.0 , LOCEAN-IPSL (2005)
48   !!----------------------------------------------------------------------
49
50CONTAINS       
51
52   SUBROUTINE dom_zgr
53      !!----------------------------------------------------------------------
54      !!                ***  ROUTINE dom_zgr  ***
55      !!                   
56      !! ** Purpose :  set the depth of model levels and the resulting
57      !!      vertical scale factors.
58      !!
59      !! ** Method  reference vertical coordinate
60      !!
61      !! ** Action :
62      !!
63      !! History :
64      !!   9.0  !  03-08  (G. Madec)  original code
65      !!   9.0  !  05-10  (A. Beckmann)  modifications for hybrid s-ccordinates
66      !!----------------------------------------------------------------------
67      INTEGER ::   ioptio = 0      ! temporary integer
68
69      NAMELIST/nam_zgr/ ln_zco, ln_zps, ln_sco
70      !!----------------------------------------------------------------------
71
72      ! Read Namelist nam_zgr : vertical coordinate'
73      ! ---------------------
74      REWIND ( numnam )
75      READ   ( numnam, nam_zgr )
76
77      ! Parameter control and print
78      ! ---------------------------
79      ! Control print
80      IF(lwp) THEN
81         WRITE(numout,*)
82         WRITE(numout,*) 'dom_zgr : vertical coordinate'
83         WRITE(numout,*) '~~~~~~~'
84         WRITE(numout,*) '          Namelist nam_zgr : set vertical coordinate'
85         WRITE(numout,*) '             z-coordinate - full steps      ln_zco = ', ln_zco
86         WRITE(numout,*) '             z-coordinate - partial steps   ln_zps = ', ln_zps
87         WRITE(numout,*) '             s- or hybrid z-s-coordinate    ln_sco = ', ln_sco
88      ENDIF
89
90      ! Check Vertical coordinate options
91      ioptio = 0
92      IF( ln_zco ) ioptio = ioptio + 1
93      IF( ln_zps ) ioptio = ioptio + 1
94      IF( ln_sco ) ioptio = ioptio + 1
95      IF ( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' )
96
97      IF( lk_zco ) THEN
98          IF(lwp) WRITE(numout,*) '          z-coordinate with reduced incore memory requirement'
99          IF( ln_zps .OR. ln_sco ) CALL ctl_stop( ' reduced memory with zps or sco option is impossible' )
100      ENDIF
101
102      ! Build the vertical coordinate system
103      ! ------------------------------------
104
105                     CALL zgr_z        ! Reference z-coordinate system (always called)
106
107                     CALL zgr_bat      ! Bathymetry fields (levels and meters)
108
109      IF( ln_zco )   CALL zgr_zco      ! z-coordinate
110
111      IF( ln_zps )   CALL zgr_zps      ! Partial step z-coordinate
112
113      IF( ln_sco )   CALL zgr_sco      ! s-coordinate or hybrid z-s coordinate
114
115!!bug gm control print:
116      IF( nprint == 1 .AND. lwp )   THEN
117         WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
118         WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ),   &
119            &                   ' w ',   MINVAL( fsdepw(:,:,:) ), '3w ', MINVAL( fsde3w(:,:,:) )
120         WRITE(numout,*) ' MIN val e3    t ', MINVAL( fse3t(:,:,:) ), ' f ', MINVAL( fse3f(:,:,:) ),  &
121            &                   ' u ',   MINVAL( fse3u(:,:,:) ), ' u ', MINVAL( fse3v(:,:,:) ),  &
122            &                   ' uw',   MINVAL( fse3uw(:,:,:)), ' vw', MINVAL( fse3vw(:,:,:)),   &
123            &                   ' w ',   MINVAL( fse3w(:,:,:) )
124
125         WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ),   &
126            &                   ' w ',   MAXVAL( fsdepw(:,:,:) ), '3w ', MAXVAL( fsde3w(:,:,:) )
127         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( fse3t(:,:,:) ), ' f ', MAXVAL( fse3f(:,:,:) ),  &
128            &                   ' u ',   MAXVAL( fse3u(:,:,:) ), ' u ', MAXVAL( fse3v(:,:,:) ),  &
129            &                   ' uw',   MAXVAL( fse3uw(:,:,:)), ' vw', MAXVAL( fse3vw(:,:,:)),   &
130            &                   ' w ',   MAXVAL( fse3w(:,:,:) )
131      ENDIF
132!!!bug gm
133
134   END SUBROUTINE dom_zgr
135
136
137   SUBROUTINE zgr_z
138      !!----------------------------------------------------------------------
139      !!                   ***  ROUTINE zgr_z  ***
140      !!                   
141      !! ** Purpose :   set the depth of model levels and the resulting
142      !!      vertical scale factors.
143      !!
144      !! ** Method  :   z-coordinate system (use in all type of coordinate)
145      !!        The depth of model levels is defined from an analytical
146      !!      function the derivative of which gives the scale factors.
147      !!        both depth and scale factors only depend on k (1d arrays).
148      !!              w-level: gdepw_0  = fsdep(k)
149      !!                       e3w_0(k) = dk(fsdep)(k)     = fse3(k)
150      !!              t-level: gdept_0  = fsdep(k+0.5)
151      !!                       e3t_0(k) = dk(fsdep)(k+0.5) = fse3(k+0.5)
152      !!
153      !! ** Action  : - gdept_0, gdepw_0 : depth of T- and W-point (m)
154      !!              -  e3t_0, e3w_0    : scale factors at T- and W-levels (m)
155      !!
156      !! Reference :
157      !!      Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
158      !!
159      !! History :
160      !!   9.0  !  03-08  (G. Madec)  F90: Free form and module
161      !!----------------------------------------------------------------------
162      !! * Local declarations
163      INTEGER  ::   jk                     ! dummy loop indices
164      REAL(wp) ::   zt, zw                 ! temporary scalars
165      REAL(wp) ::   &
166      zsur , za0, za1, zkth, zacr,      &  ! Values set from parameters in
167      zdzmin, zhmax                        ! par_CONFIG_Rxx.h90
168      !!----------------------------------------------------------------------
169
170      ! Set variables from parameters
171      ! ------------------------------
172       zkth = ppkth       ;   zacr = ppacr
173       zdzmin = ppdzmin   ;   zhmax = pphmax
174
175      ! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
176      !  za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
177      !
178       IF(  ppa1  == pp_to_be_computed  .AND.  &
179         &  ppa0  == pp_to_be_computed  .AND.  &
180         &  ppsur == pp_to_be_computed           ) THEN
181         za1 = ( ppdzmin - pphmax / FLOAT(jpk-1) )          &
182             / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) &
183             &                         *  (  LOG( COSH( (jpk - ppkth) / ppacr) )      &
184             &                             - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
185
186         za0  = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
187         zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
188
189       ELSE
190         za1 = ppa1 ;       za0 = ppa0 ;          zsur = ppsur
191       ENDIF
192
193
194      ! Parameter print
195      ! ---------------
196      IF(lwp) THEN
197         WRITE(numout,*)
198         WRITE(numout,*) '    zgr_z   : Reference vertical z-coordinates'
199         WRITE(numout,*) '    ~~~~~~~'
200         IF(  ppkth == 0. ) THEN             
201              WRITE(numout,*) '            Uniform grid with ',jpk-1,' layers'
202              WRITE(numout,*) '            Total depth    :', zhmax
203              WRITE(numout,*) '            Layer thickness:', zhmax/(jpk-1)
204         ELSE
205            IF( ppa1 == 0. .AND. ppa0 == 0. .AND. ppsur == 0. ) THEN
206               WRITE(numout,*) '         zsur, za0, za1 computed from '
207               WRITE(numout,*) '                 zdzmin = ', zdzmin
208               WRITE(numout,*) '                 zhmax  = ', zhmax
209            ENDIF
210            WRITE(numout,*) '           Value of coefficients for vertical mesh:'
211            WRITE(numout,*) '                 zsur = ', zsur
212            WRITE(numout,*) '                 za0  = ', za0
213            WRITE(numout,*) '                 za1  = ', za1
214            WRITE(numout,*) '                 zkth = ', zkth
215            WRITE(numout,*) '                 zacr = ', zacr
216         ENDIF
217      ENDIF
218
219
220      ! Reference z-coordinate (depth - scale factor at T- and W-points)
221      ! ======================
222      IF( ppkth == 0.e0 ) THEN            !  uniform vertical grid       
223
224         za1 = zhmax / FLOAT(jpk-1) 
225         DO jk = 1, jpk
226            zw = FLOAT( jk )
227            zt = FLOAT( jk ) + 0.5
228            gdepw_0(jk) = ( zw - 1 ) * za1
229            gdept_0(jk) = ( zt - 1 ) * za1
230            e3w_0  (jk) =  za1
231            e3t_0  (jk) =  za1
232         END DO
233
234      ELSE
235
236         DO jk = 1, jpk
237            zw = FLOAT( jk )
238            zt = FLOAT( jk ) + 0.5
239            gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG( COSH( (zw-zkth)/zacr ) )  )
240            gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG( COSH( (zt-zkth)/zacr ) )  )
241            e3w_0  (jk) =          za0      + za1        * TANH(      (zw-zkth)/zacr   )
242            e3t_0  (jk) =          za0      + za1        * TANH(      (zt-zkth)/zacr   )
243         END DO
244         gdepw_0(1) = 0.e0      ! force first w-level to be exactly at zero
245
246      ENDIF
247
248      ! Control and  print
249      ! ==================
250
251      IF(lwp) THEN
252         WRITE(numout,*)
253         WRITE(numout,*) '              Reference z-coordinate depth and scale factors:'
254         WRITE(numout, "(9x,' level   gdept    gdepw     e3t      e3w  ')" )
255         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_0(jk), gdepw_0(jk), e3t_0(jk), e3w_0(jk), jk = 1, jpk )
256      ENDIF
257
258      DO jk = 1, jpk
259         IF( e3w_0(jk)  <= 0. .OR. e3t_0(jk)  <= 0. ) CALL ctl_stop( ' e3w or e3t =< 0 ' )
260         IF( gdepw_0(jk) < 0. .OR. gdept_0(jk) < 0. ) CALL ctl_stop( ' gdepw or gdept < 0 ' )
261      END DO
262
263   END SUBROUTINE zgr_z
264
265
266   SUBROUTINE zgr_bat
267      !!----------------------------------------------------------------------
268      !!                    ***  ROUTINE zgr_bat  ***
269      !!
270      !! ** Purpose :   set bathymetry both in levels and meters
271      !!
272      !! ** Method  :   read or define mbathy and bathy arrays
273      !!       * level bathymetry:
274      !!      The ocean basin geometry is given by a two-dimensional array,
275      !!      mbathy, which is defined as follow :
276      !!            mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
277      !!                              at t-point (ji,jj).
278      !!                            = 0  over the continental t-point.
279      !!                            = -n over the nth island t-point.
280      !!      The array mbathy is checked to verified its consistency with
281      !!      model option. in particular:
282      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
283      !!                  along closed boundary.
284      !!            mbathy must be cyclic IF jperio=1.
285      !!            mbathy must be lower or equal to jpk-1.
286      !!            isolated ocean grid points are suppressed from mbathy
287      !!                  since they are only connected to remaining
288      !!                  ocean through vertical diffusion.
289      !!      ntopo=-1 :   rectangular channel or bassin with a bump
290      !!      ntopo= 0 :   flat rectangular channel or basin
291      !!      ntopo= 1 :   mbathy is read in 'bathy_level.nc' NetCDF file
292      !!                   bathy  is read in 'bathy_meter.nc' NetCDF file
293      !!      C A U T I O N : mbathy will be modified during the initializa-
294      !!      tion phase to become the number of non-zero w-levels of a water
295      !!      column, with a minimum value of 1.
296      !!
297      !! ** Action  : - mbathy: level bathymetry (in level index)
298      !!              - bathy : meter bathymetry (in meters)
299      !!
300      !! History :
301      !!   9.0  !  03-08  (G. Madec)  Original code
302      !!   9.0  !  05-10  (A. Beckmann)  modifications for s-ccordinates
303      !!----------------------------------------------------------------------
304      !! * Modules used
305      USE iom
306
307      !! * Local declarations
308      INTEGER ::   ji, jj, jl, jk       ! dummy loop indices
309      INTEGER ::   inum                 ! temporary logical unit
310      INTEGER  ::   &
311         ii_bump, ij_bump, ih           ! bump center position
312      INTEGER , DIMENSION(jpidta,jpjdta) ::   &
313         idta                           ! global domain integer data
314      REAL(wp) ::   &
315         r_bump, h_bump, h_oce,      &  ! bump characteristics
316         zi, zj, zh                     ! temporary scalars
317      REAL(wp), DIMENSION(jpidta,jpjdta) ::   &
318         zdta                           ! global domain scalar data
319      !!----------------------------------------------------------------------
320
321      IF(lwp) WRITE(numout,*)
322      IF(lwp) WRITE(numout,*) '    zgr_bat : defines level and meter bathymetry'
323      IF(lwp) WRITE(numout,*) '    ~~~~~~~'
324
325      ! ========================================
326      ! global domain level and meter bathymetry (idta,zdta)
327      ! ========================================
328      !                                               ! =============== !
329      IF( ntopo == 0 .OR. ntopo == -1 ) THEN          ! defined by hand !
330         !                                            ! =============== !
331
332         IF( ntopo == 0 ) THEN                        ! flat basin
333
334            IF(lwp) WRITE(numout,*)
335            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin'
336
337            idta(:,:) = jpkm1                            ! flat basin
338            zdta(:,:) = gdepw_0(jpk)
339            h_oce     = gdepw_0(jpk)
340
341         ELSE                                         ! bump
342            IF(lwp) WRITE(numout,*)
343            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin with a bump'
344
345            ii_bump = jpidta / 2           ! i-index of the bump center
346            ij_bump = jpjdta / 2           ! j-index of the bump center
347            r_bump  = 50000.e0             ! bump radius (meters)       
348            h_bump  = 2700.e0              ! bump height (meters)
349            h_oce   = gdepw_0(jpk)         ! background ocean depth (meters)
350            IF(lwp) WRITE(numout,*) '            bump characteristics: '
351            IF(lwp) WRITE(numout,*) '               bump center (i,j)   = ', ii_bump, ii_bump
352            IF(lwp) WRITE(numout,*) '               bump height         = ', h_bump , ' meters'
353            IF(lwp) WRITE(numout,*) '               bump radius         = ', r_bump , ' index'
354            IF(lwp) WRITE(numout,*) '            background ocean depth = ', h_oce  , ' meters'
355            ! zdta :
356            DO jj = 1, jpjdta
357               DO ji = 1, jpidta
358                  zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
359                  zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
360                  zdta(ji,jj) = h_oce - h_bump * EXP( -( zi*zi + zj*zj ) )
361               END DO
362            END DO
363
364            ! idta :
365            IF( ln_sco ) THEN      ! s-coordinate (zsc       ): idta()=jpk
366               idta(:,:) = jpkm1
367            ELSE                   ! z-coordinate (zco or zps): step-like topography
368               idta(:,:) = jpkm1
369               DO jk = 1, jpkm1
370                  DO jj = 1, jpjdta
371                     DO ji = 1, jpidta
372                        IF( gdept_0(jk) < zdta(ji,jj) .AND. zdta(ji,jj) <= gdept_0(jk+1) )   idta(ji,jj) = jk
373                     END DO
374                  END DO
375               END DO
376            ENDIF
377         ENDIF
378
379         ! set boundary conditions (caution, idta on the global domain: use of jperio, not nperio)
380         IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
381            idta( :    , 1    ) = -1                ;      zdta( :    , 1    ) = -1.e0
382            idta( :    ,jpjdta) =  0                ;      zdta( :    ,jpjdta) =  0.e0
383         ELSEIF( jperio == 2 ) THEN
384            idta( :    , 1    ) = idta( : ,  3  )   ;      zdta( :    , 1    ) = zdta( : ,  3  )
385            idta( :    ,jpjdta) = 0                 ;      zdta( :    ,jpjdta) =  0.e0
386            idta( 1    , :    ) = 0                 ;      zdta( 1    , :    ) =  0.e0
387            idta(jpidta, :    ) = 0                 ;      zdta(jpidta, :    ) =  0.e0
388         ELSE
389            ih = 0                                  ;      zh = 0.e0
390            IF( ln_sco )   ih = jpkm1               ;      IF( ln_sco )   zh = h_oce
391            idta( :    , 1    ) = ih                ;      zdta( :    , 1    ) =  zh
392            idta( :    ,jpjdta) = ih                ;      zdta( :    ,jpjdta) =  zh
393            idta( 1    , :    ) = ih                ;      zdta( 1    , :    ) =  zh
394            idta(jpidta, :    ) = ih                ;      zdta(jpidta, :    ) =  zh
395         ENDIF
396
397         ! =======================================
398         ! local domain level and meter bathymetry (mbathy,bathy)
399         ! =======================================
400         
401         mbathy(:,:) = 0                                 ! set to zero extra halo points
402         bathy (:,:) = 0.e0                              ! (require for mpp case)
403         
404         DO jj = 1, nlcj                                 ! interior values
405            DO ji = 1, nlci
406               mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
407               bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
408            END DO
409         END DO
410
411         !                                            ! =============== !
412      ELSEIF( ntopo == 1 ) THEN                       !   read in file  !
413         !                                            ! =============== !
414
415         IF( ln_zco )   THEN
416            CALL iom_open ( 'bathy_level.nc', inum )   ! Level bathymetry
417            CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy )
418            CALL iom_close (inum)
419            mbathy(:,:) = INT( bathy(:,:) )
420         ENDIF
421
422         IF( ln_zps .OR. ln_sco )   THEN
423            CALL iom_open ( 'bathy_meter.nc', inum )   ! meter bathymetry
424            CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy )
425            CALL iom_close (inum)
426         ENDIF
427         !                                            ! =============== !
428      ELSE                                            !      error      !
429         !                                            ! =============== !
430         WRITE(ctmp1,*) '          parameter , ntopo = ', ntopo
431         CALL ctl_stop( '    zgr_bat : '//trim(ctmp1) )
432      ENDIF
433
434      ! =======================
435      ! NO closed seas or lakes
436      ! =======================
437
438      IF( nclosea == 0 ) THEN
439         DO jl = 1, jpncs
440            DO jj = ncsj1(jl), ncsj2(jl)
441               DO ji = ncsi1(jl), ncsi2(jl)
442                  mbathy(ji,jj) = 0                   ! suppress closed seas
443                  bathy (ji,jj) = 0.e0                ! and lakes
444!!bug gm          bathy (ji,jj) = 300.e0                ! and lakes
445               END DO
446            END DO
447         END DO
448      ENDIF
449
450#if defined key_orca_lev10
451      ! 10 time the vertical resolution
452      mbathy(:,:) = 10 * mbathy(:,:)
453      IF(lwp) WRITE(numout,*) ' ATTENTION: 300 niveaux avec bathy levels "vraie?"'
454#endif
455      ! ===========
456      ! Zoom domain
457      ! ===========
458
459      IF( lzoom )   CALL zgr_bat_zoom
460
461      ! ================
462      ! Bathymetry check
463      ! ================
464
465      CALL zgr_bat_ctl
466
467   END SUBROUTINE zgr_bat
468
469
470   SUBROUTINE zgr_bat_zoom
471      !!----------------------------------------------------------------------
472      !!                    ***  ROUTINE zgr_bat_zoom  ***
473      !!
474      !! ** Purpose : - Close zoom domain boundary if necessary
475      !!              - Suppress Med Sea from ORCA R2 and R05 arctic zoom
476      !!
477      !! ** Method  :
478      !!
479      !! ** Action  : - update mbathy: level bathymetry (in level index)
480      !!
481      !! History :
482      !!   9.0  !  03-08  (G. Madec)  Original code
483      !!----------------------------------------------------------------------
484      !! * Local variables
485      INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integers
486      !!----------------------------------------------------------------------
487
488      IF(lwp) WRITE(numout,*)
489      IF(lwp) WRITE(numout,*) '    zgr_bat_zoom : modify the level bathymetry for zoom domain'
490      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~'
491
492      ! Zoom domain
493      ! ===========
494
495      ! Forced closed boundary if required
496      IF( lzoom_w )   mbathy( mi0(jpizoom):mi1(jpizoom) , :  ) = 0
497      IF( lzoom_s )   mbathy(  : , mj0(jpjzoom):mj1(jpjzoom) ) = 0
498      IF( lzoom_e )   mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , :  ) = 0
499      IF( lzoom_n )   mbathy(  : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0
500
501      ! Configuration specific domain modifications
502      ! (here, ORCA arctic configuration: suppress Med Sea)
503      IF( cp_cfg == "orca" .AND. lzoom_arct ) THEN
504         SELECT CASE ( jp_cfg )
505         !                                        ! =======================
506         CASE ( 2 )                               !  ORCA_R2 configuration
507            !                                     ! =======================
508            IF(lwp) WRITE(numout,*) '                   ORCA R2 arctic zoom: suppress the Med Sea'
509
510            ii0 = 141   ;   ii1 = 162      ! Sea box i,j indices
511            ij0 =  98   ;   ij1 = 110
512            !                                     ! =======================
513         CASE ( 05 )                              !  ORCA_R05 configuration
514            !                                     ! =======================
515            IF(lwp) WRITE(numout,*) '                   ORCA R05 arctic zoom: suppress the Med Sea'
516 
517            ii0 = 563   ;   ii1 = 642      ! zero over the Med Sea boxe
518            ij0 = 314   ;   ij1 = 370 
519
520         END SELECT
521         !
522         mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0   ! zero over the Med Sea boxe
523         !
524      ENDIF
525
526   END SUBROUTINE zgr_bat_zoom
527
528
529   SUBROUTINE zgr_bat_ctl
530      !!----------------------------------------------------------------------
531      !!                    ***  ROUTINE zgr_bat_ctl  ***
532      !!
533      !! ** Purpose :   check the bathymetry in levels
534      !!
535      !! ** Method  :   The array mbathy is checked to verified its consistency
536      !!      with the model options. in particular:
537      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
538      !!                  along closed boundary.
539      !!            mbathy must be cyclic IF jperio=1.
540      !!            mbathy must be lower or equal to jpk-1.
541      !!            isolated ocean grid points are suppressed from mbathy
542      !!                  since they are only connected to remaining
543      !!                  ocean through vertical diffusion.
544      !!      C A U T I O N : mbathy will be modified during the initializa-
545      !!      tion phase to become the number of non-zero w-levels of a water
546      !!      column, with a minimum value of 1.
547      !!
548      !! ** Action  : - update mbathy: level bathymetry (in level index)
549      !!              - update bathy : meter bathymetry (in meters)
550      !!
551      !! History :
552      !!   9.0  !  03-08  (G. Madec)  Original code
553      !!   9.0  !  05-10  (A. Beckmann)  modifications for s-ccordinates
554      !!----------------------------------------------------------------------
555      !! * Local declarations
556      INTEGER ::   ji, jj, jl           ! dummy loop indices
557      INTEGER ::   &
558         icompt, ibtest, ikmax          ! temporary integers
559      REAL(wp), DIMENSION(jpi,jpj) ::   &
560         zbathy                         ! temporary workspace
561      !!----------------------------------------------------------------------
562
563      IF(lwp) WRITE(numout,*)
564      IF(lwp) WRITE(numout,*) '    zgr_bat_ctl : check the bathymetry'
565      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
566
567      ! ================
568      ! Bathymetry check
569      ! ================
570
571      IF( .NOT. lk_c1d )   THEN
572
573         ! Suppress isolated ocean grid points
574
575         IF(lwp) WRITE(numout,*)
576         IF(lwp) WRITE(numout,*)'                   suppress isolated ocean grid points'
577         IF(lwp) WRITE(numout,*)'                   -----------------------------------'
578
579         icompt = 0
580
581         DO jl = 1, 2
582
583            IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
584               mbathy( 1 ,:) = mbathy(jpim1,:)
585               mbathy(jpi,:) = mbathy(  2  ,:)
586            ENDIF
587            DO jj = 2, jpjm1
588               DO ji = 2, jpim1
589                  ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj),   &
590                     &          mbathy(ji,jj-1),mbathy(ji,jj+1) )
591                  IF( ibtest < mbathy(ji,jj) ) THEN
592                     IF(lwp) WRITE(numout,*) ' the number of ocean level at ',   &
593                        &   'grid-point (i,j) =  ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
594                     mbathy(ji,jj) = ibtest
595                     icompt = icompt + 1
596                  ENDIF
597               END DO
598            END DO
599
600         END DO
601         IF( icompt == 0 ) THEN
602            IF(lwp) WRITE(numout,*)'     no isolated ocean grid points'
603         ELSE
604            IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed'
605         ENDIF
606         IF( lk_mpp ) THEN
607            zbathy(:,:) = FLOAT( mbathy(:,:) )
608            CALL lbc_lnk( zbathy, 'T', 1. )
609            mbathy(:,:) = INT( zbathy(:,:) )
610         ENDIF
611
612         ! 3.2 East-west cyclic boundary conditions
613
614         IF( nperio == 0 ) THEN
615            IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west',   &
616               ' boundary: nperio = ', nperio
617            IF( lk_mpp ) THEN
618               IF( nbondi == -1 .OR. nbondi == 2 ) THEN
619                  IF( jperio /= 1 )   mbathy(1,:) = 0
620               ENDIF
621               IF( nbondi == 1 .OR. nbondi == 2 ) THEN
622                  IF( jperio /= 1 )   mbathy(nlci,:) = 0
623               ENDIF
624            ELSE
625               IF( ln_zco .OR. ln_zps ) THEN
626                  mbathy( 1 ,:) = 0
627                  mbathy(jpi,:) = 0
628               ELSE
629                  mbathy( 1 ,:) = jpkm1
630                  mbathy(jpi,:) = jpkm1
631               ENDIF
632            ENDIF
633         ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio ==  6 ) THEN
634            IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions',   &
635               ' on mbathy: nperio = ', nperio
636            mbathy( 1 ,:) = mbathy(jpim1,:)
637            mbathy(jpi,:) = mbathy(  2  ,:)
638         ELSEIF( nperio == 2 ) THEN
639            IF(lwp) WRITE(numout,*) '   equatorial boundary conditions',   &
640               ' on mbathy: nperio = ', nperio
641         ELSE
642            IF(lwp) WRITE(numout,*) '    e r r o r'
643            IF(lwp) WRITE(numout,*) '    parameter , nperio = ', nperio
644            !         STOP 'dom_mba'
645         ENDIF
646
647         ! Set to zero mbathy over islands if necessary  (lk_isl=F)
648         IF( .NOT. lk_isl ) THEN    ! No island
649            IF(lwp) WRITE(numout,*)
650            IF(lwp) WRITE(numout,*) '         mbathy set to 0 over islands'
651            IF(lwp) WRITE(numout,*) '         ----------------------------'
652
653            mbathy(:,:) = MAX( 0, mbathy(:,:) )
654
655            !  Boundary condition on mbathy
656            IF( .NOT.lk_mpp ) THEN 
657               !!bug ???  y reflechir!
658               !   ... mono- or macro-tasking: T-point, >0, 2D array, no slab
659               zbathy(:,:) = FLOAT( mbathy(:,:) )
660               CALL lbc_lnk( zbathy, 'T', 1. )
661               mbathy(:,:) = INT( zbathy(:,:) )
662            ENDIF
663
664         ENDIF
665
666      ENDIF
667
668      ! Number of ocean level inferior or equal to jpkm1
669
670      ikmax = 0
671      DO jj = 1, jpj
672         DO ji = 1, jpi
673            ikmax = MAX( ikmax, mbathy(ji,jj) )
674         END DO
675      END DO
676      !!! test a faire:   ikmax = MAX( mbathy(:,:) )   ???
677
678      IF( ikmax > jpkm1 ) THEN
679         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' >  jpk-1'
680         IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
681      ELSE IF( ikmax < jpkm1 ) THEN
682         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1' 
683         IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
684      ENDIF
685
686      IF( lwp .AND. nprint == 1 ) THEN
687         WRITE(numout,*)
688         WRITE(numout,*) ' bathymetric field '
689         WRITE(numout,*) ' ------------------'
690         WRITE(numout,*) ' number of non-zero T-levels '
691         CALL prihin( mbathy, jpi, jpj, 1, jpi,   &
692                      1     , 1  , jpj, 1, 3  ,   &
693                      numout )
694         WRITE(numout,*)
695      ENDIF
696
697   END SUBROUTINE zgr_bat_ctl
698
699
700   SUBROUTINE zgr_zco
701      !!----------------------------------------------------------------------
702      !!                  ***  ROUTINE zgr_zco  ***
703      !!
704      !! ** Purpose :   define the z-coordinate system
705      !!
706      !! ** Method  :   set 3D coord. arrays to reference 1D array if lk_zco=T
707      !!
708      !! History :
709      !!        !  06-04  (R. Benshila, G. Madec)  Original code
710      !!----------------------------------------------------------------------
711      INTEGER  ::   jk
712      !!----------------------------------------------------------------------
713
714      IF( .NOT.lk_zco ) THEN
715         DO jk = 1, jpk
716            fsdept(:,:,jk) = gdept_0(jk)
717            fsdepw(:,:,jk) = gdepw_0(jk)
718            fsde3w(:,:,jk) = gdepw_0(jk)
719            fse3t (:,:,jk) = e3t_0(jk)
720            fse3u (:,:,jk) = e3t_0(jk)
721            fse3v (:,:,jk) = e3t_0(jk)
722            fse3f (:,:,jk) = e3t_0(jk)
723            fse3w (:,:,jk) = e3w_0(jk)
724            fse3uw(:,:,jk) = e3w_0(jk)
725            fse3vw(:,:,jk) = e3w_0(jk)
726         END DO
727      ENDIF
728
729   END SUBROUTINE zgr_zco
730
731
732
733#  include "domzgr_zps.h90"
734
735
736#if ! defined key_zco
737   !!----------------------------------------------------------------------
738   !!   NOT 'key_zco' :                                    3D gdep. and e3.
739   !!----------------------------------------------------------------------
740
741   SUBROUTINE zgr_sco
742      !!----------------------------------------------------------------------
743      !!                  ***  ROUTINE zgr_sco  ***
744      !!                     
745      !! ** Purpose :   define the s-coordinate system
746      !!
747      !! ** Method  :   s-coordinate
748      !!         The depth of model levels is defined as the product of an
749      !!      analytical function by the local bathymetry, while the vertical
750      !!      scale factors are defined as the product of the first derivative
751      !!      of the analytical function by the bathymetry.
752      !!      (this solution save memory as depth and scale factors are not
753      !!      3d fields)
754      !!          - Read bathymetry (in meters) at t-point and compute the
755      !!         bathymetry at u-, v-, and f-points.
756      !!            hbatu = mi( hbatt )
757      !!            hbatv = mj( hbatt )
758      !!            hbatf = mi( mj( hbatt ) )
759      !!          - Compute gsigt, gsigw, esigt, esigw from an analytical
760      !!         function and its derivative given as statement function.
761      !!            gsigt(k) = fssig (k+0.5)
762      !!            gsigw(k) = fssig (k    )
763      !!            esigt(k) = fsdsig(k+0.5)
764      !!            esigw(k) = fsdsig(k    )
765      !!      This routine is given as an example, it must be modified
766      !!      following the user s desiderata. nevertheless, the output as
767      !!      well as the way to compute the model levels and scale factors
768      !!      must be respected in order to insure second order a!!uracy
769      !!      schemes.
770      !!
771      !! Reference :
772      !!      Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
773      !!
774      !! History :
775      !!        !  95-12  (G. Madec)  Original code : s vertical coordinate
776      !!        !  97-07  (G. Madec)  lbc_lnk call
777      !!        !  97-04  (J.-O. Beismann)
778      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
779      !!   9.0  !  05-10  (A. Beckmann)  new stretching function
780      !!----------------------------------------------------------------------
781      !! * Local declarations
782      INTEGER  ::   ji, jj, jk, jl
783      REAL(wp) ::   fssig, fsdsig, pfkk
784
785      INTEGER  ::   iip1, ijp1, iim1, ijm1
786      REAL(wp) ::   &
787         fssigt, fssigw, zcoeft, zcoefw,   &
788         zrmax, ztaper
789
790      REAL(wp), DIMENSION(jpi,jpj) ::   &
791         zenv, ztmp, zmsk, zri, zrj, zhbat
792
793      NAMELIST/nam_zgr_sco/ sbot_max, sbot_min, theta, thetb, r_max
794      !!----------------------------------------------------------------------
795
796      ! a. Sigma stretching coefficient
797       fssigt(pfkk) = (tanh(theta*(-(pfkk-0.5)/FLOAT(jpkm1)+thetb))-tanh(thetb*theta)) &
798                    * (cosh(theta)+cosh(theta*(2.*thetb-1.)))/(2.*sinh(theta))
799       fssigw(pfkk) = (tanh(theta*(-(pfkk-1.0)/FLOAT(jpkm1)+thetb))-tanh(thetb*theta)) &
800                    * (cosh(theta)+cosh(theta*(2.*thetb-1.)))/(2.*sinh(theta))
801
802      ! b. Vertical derivative of sigma stretching coefficient
803 
804!!bug fsdsig(pfkk)= put here the analytical derivative of fssigt and w ...
805      !!----------------------------------------------------------------------
806
807      ! Read Namelist nam_zgr_sco : sigma-stretching parameters
808      ! -------------------------
809      REWIND( numnam )
810      READ  ( numnam, nam_zgr_sco )
811
812      IF(lwp) THEN
813         WRITE(numout,*)
814         WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate'
815         WRITE(numout,*) '~~~~~~~~~~~'
816         WRITE(numout,*) '         Namelist nam_zgr_sco'
817         WRITE(numout,*) '            sigma-stretching coeffs '
818         WRITE(numout,*) '            maximum depth of s-bottom surface (>0)    sbot_max = ' ,sbot_max
819         WRITE(numout,*) '            minimum depth of s-bottom surface (>0)    sbot_min = ' ,sbot_min
820         WRITE(numout,*) '            surface control parameter (0<=theta<=20)  theta    = ', theta
821         WRITE(numout,*) '            bottom control parameter  (0<=thetb<= 1)  thetb    = ', thetb
822         WRITE(numout,*) '            maximum cut-off r-value allowed           r_max    = ' , r_max
823      ENDIF
824
825      ! ???
826      ! ------------------
827      hift(:,:) = sbot_min
828      hifu(:,:) = sbot_min
829      hifv(:,:) = sbot_min
830      hiff(:,:) = sbot_min
831
832
833      ! set maximum ocean depth
834      ! -----------------------
835       DO jj = 1, jpj
836          DO ji = 1, jpi
837             bathy(ji,jj) = MIN( sbot_max, bathy(ji,jj) )
838          END DO
839       END DO
840
841      ! Define the envelop bathymetry
842      ! =============================
843      ! Smooth the bathymetry (if required)
844
845      scosrf(:,:) = 0.e0              ! ocean surface depth (here zero: no under ice-shelf sea)
846      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
847
848 
849      ! use r-value to create hybrid coordinates
850 
851      DO jj = 1, jpj
852         DO ji = 1, jpi
853            zenv(ji,jj) = MAX( bathy(ji,jj), sbot_min )
854         END DO
855      END DO
856
857      jl = 0
858      zrmax = 1.e0
859      !                                                     ! ================ !
860      DO WHILE ( jl <= 10000 .AND. zrmax > r_max )          !  Iterative loop  !
861         !                                                  ! ================ !
862         jl = jl + 1
863
864         zrmax = 0.e0
865         zmsk(:,:) = 0.e0
866 
867         DO jj = 1, nlcj
868            DO ji = 1, nlci
869               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
870               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
871               zri(ji,jj) = ABS( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
872               zrj(ji,jj) = ABS( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
873               zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) )
874               IF( zri(ji,jj) > r_max )   zmsk(ji  ,jj  ) = 1.0
875               IF( zri(ji,jj) > r_max )   zmsk(iip1,jj  ) = 1.0
876               IF( zrj(ji,jj) > r_max )   zmsk(ji  ,jj  ) = 1.0
877               IF( zrj(ji,jj) > r_max )   zmsk(ji  ,ijp1) = 1.0
878            END DO
879         END DO
880         IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
881         ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)
882         ztmp(:,:) = zmsk(:,:)
883         CALL lbc_lnk( zmsk, 'T', 1. )
884         DO jj = 1, nlcj
885            DO ji = 1, nlci
886                zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )   
887            END DO
888         END DO
889 
890         WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) )
891
892         DO jj = 1, nlcj
893            DO ji = 1, nlci
894               iip1 = MIN( ji+1, nlci )     ! last  line (ji=nlci)
895               ijp1 = MIN( jj+1, nlcj )     ! last  raw  (jj=nlcj)
896               iim1 = MAX( ji-1,  1  )      ! first line (ji=nlci)
897               ijm1 = MAX( jj-1,  1  )      ! first raw  (jj=nlcj)
898               IF( zmsk(ji,jj) == 1.0 ) THEN
899                  ztmp(ji,jj) =   (                                                                                   &
900             &      zenv(iim1,ijp1)*zmsk(iim1,ijp1) + zenv(ji,ijp1)*zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1)   &
901             &    + zenv(iim1,jj  )*zmsk(iim1,jj  ) + zenv(ji,jj  )*    2.e0      + zenv(iip1,jj  )*zmsk(iip1,jj  )   &
902             &    + zenv(iim1,ijm1)*zmsk(iim1,ijm1) + zenv(ji,ijm1)*zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1)   &
903             &                    ) / (                                                                               &
904             &                      zmsk(iim1,ijp1) +               zmsk(ji,ijp1) +                 zmsk(iip1,ijp1)   &
905             &    +                 zmsk(iim1,jj  ) +                   2.e0      +                 zmsk(iip1,jj  )   &
906             &    +                 zmsk(iim1,ijm1) +               zmsk(ji,ijm1) +                 zmsk(iip1,ijm1)   &
907             &                        )
908               ENDIF
909            END DO
910         END DO
911 
912         DO jj = 1, nlcj
913            DO ji = 1, nlci
914               IF( zmsk(ji,jj) == 1.0 )   zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )
915            END DO
916         END DO
917
918         ! Apply lateral boundary condition   CAUTION: kept the value when the lbc field is zero
919         ztmp(:,:) = zenv(:,:)
920         CALL lbc_lnk( zenv, 'T', 1. )
921         DO jj = 1, nlcj
922            DO ji = 1, nlci
923               IF( zenv(ji,jj) == 0.e0 )   zenv(ji,jj) = ztmp(ji,jj)
924            END DO
925         END DO
926         !                                                  ! ================ !
927      END DO                                                !     End loop     !
928      !                                                     ! ================ !
929
930      ! save envelop bathymetry in hbatt
931      hbatt(:,:) = zenv(:,:) 
932
933!!bug gm  :   CAUTION:  tapering hard coded !!!!  orca2 only
934!!bug gm                NOT valid in mpp ===> taper have been changed
935
936       DO jj = 1, jpj
937          DO ji = 1, jpi
938!!bug mpp    ztaper = EXP( -(FLOAT(jj-74)/10.)**2 )
939             ztaper = EXP( -(gphit(ji,jj)/8)**2 )
940             hbatt(ji,jj) = sbot_max * ztaper + hbatt(ji,jj) * (1.-ztaper)
941          END DO
942       END DO
943
944      DO jj = 1, jpj
945         zrmax  = EXP( -(gphit(10,jj)/8)**2 )
946         ztaper = EXP( -(FLOAT(jj-74)/10.)**2 )
947         IF( nprint == 1 .AND. lwp ) WRITE(numout,*) jj, (FLOAT(jj-74)/10.), ztaper,(gphit(10,jj)/8), zrmax
948      END DO
949
950      IF( nprint == 1 .AND. lwp )   THEN
951         WRITE(numout,*) ' bathy  MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
952         WRITE(numout,*) ' hbatt  MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
953      ENDIF
954
955      ! Control print
956      IF(lwp) THEN
957         WRITE(numout,*)
958         WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
959         WRITE(numout,*)
960         CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0., numout )
961      ENDIF
962
963      ! 1. computation of hbatu, hbatv, hbatf fields
964      ! --------------------------------------------
965
966      hbatu(:,:) = sbot_min
967      hbatv(:,:) = sbot_min
968      hbatf(:,:) = sbot_min
969
970      IF(lwp) THEN
971         WRITE(numout,*)
972         WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', sbot_min
973         WRITE(numout,*)
974      ENDIF
975
976      DO jj = 1, jpjm1
977        DO ji = 1, jpim1
978           hbatu(ji,jj)= 0.5 *( hbatt(ji  ,jj)+hbatt(ji+1,jj  ) )
979           hbatv(ji,jj)= 0.5 *( hbatt(ji  ,jj)+hbatt(ji  ,jj+1) )
980           hbatf(ji,jj)= 0.25*( hbatt(ji  ,jj)+hbatt(ji  ,jj+1)   &
981                               +hbatt(ji+1,jj)+hbatt(ji+1,jj+1) )
982        END DO
983      END DO
984 
985      ! Apply lateral boundary condition
986      ! CAUTION: retain non zero value in the initial file
987      !          this should be OK for orca cfg, not for EEL
988      zhbat(:,:) = hbatu(:,:)
989      CALL lbc_lnk( hbatu, 'U', 1. )
990      DO jj = 1, jpj
991         DO ji = 1, jpi
992            IF( hbatu(ji,jj) == 0.e0 ) THEN
993               IF( zhbat(ji,jj) == 0.e0 )   hbatu(ji,jj) = sbot_min
994               IF( zhbat(ji,jj) /= 0.e0 )   hbatu(ji,jj) = zhbat(ji,jj)
995            ENDIF
996         END DO
997      END DO
998
999      zhbat(:,:) = hbatv(:,:)
1000      CALL lbc_lnk( hbatv, 'V', 1. )
1001      DO jj = 1, jpj
1002         DO ji = 1, jpi
1003            IF( hbatv(ji,jj) == 0.e0 ) THEN
1004               IF( zhbat(ji,jj) == 0.e0 )   hbatv(ji,jj) = sbot_min
1005               IF( zhbat(ji,jj) /= 0.e0 )   hbatv(ji,jj) = zhbat(ji,jj)
1006            ENDIF
1007         END DO
1008      END DO
1009
1010      zhbat(:,:) = hbatf(:,:)
1011      CALL lbc_lnk( hbatf, 'F', 1. )
1012      DO jj = 1, jpj
1013         DO ji = 1, jpi
1014            IF( hbatf(ji,jj) == 0.e0 ) THEN
1015               IF( zhbat(ji,jj) == 0.e0 )   hbatf(ji,jj) = sbot_min
1016               IF( zhbat(ji,jj) /= 0.e0 )   hbatf(ji,jj) = zhbat(ji,jj)
1017            ENDIF
1018         END DO
1019      END DO
1020
1021
1022!!bug:  key_helsinki a verifer
1023      hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1024      hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1025      hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1026      hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1027
1028      IF( nprint == 1 .AND. lwp )   THEN
1029         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift(:,:) ), ' f ', MAXVAL( hiff(:,:) ),  &
1030            &                        ' u ',   MAXVAL( hifu(:,:) ), ' v ', MAXVAL( hifv(:,:) )
1031         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift(:,:) ), ' f ', MINVAL( hiff(:,:) ),  &
1032            &                        ' u ',   MINVAL( hifu(:,:) ), ' v ', MINVAL( hifv(:,:) )
1033         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
1034            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1035         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
1036            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1037      ENDIF
1038!! helsinki
1039
1040      ! 2. Computation of gsig and esig fields
1041      ! --------------------------------------
1042
1043      ! 2.1 Coefficients for model level depth at w- and t-levels
1044
1045!!bug gm : change the def of statment function....
1046      DO jk = 1, jpk
1047        gsigw(jk) = -fssigw(FLOAT(jk))
1048        gsigt(jk) = -fssigt(FLOAT(jk))
1049      END DO
1050
1051      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk)
1052
1053!!org DO jk = 1, jpk
1054!!org    gsigw(jk) = -fssig ( FLOAT(jk)    )
1055!!org    gsigt(jk) = -fssig ( FLOAT(jk)+0.5)
1056!!org END DO
1057
1058      ! 2.2 Coefficients for vertical scale factors at w-, t- levels
1059
1060!!bug gm :  define it from analytical function, not like juste bellow....
1061!!          or betteroffer the 2 possibilities....
1062      DO jk=2,jpk
1063        esigw(jk)=gsigt(jk)-gsigt(jk-1)
1064      END DO
1065      DO jk=1,jpkm1
1066        esigt(jk)=gsigw(jk+1)-gsigw(jk)
1067      END DO
1068        esigw(1)=esigw(2)
1069        esigt(jpk)=esigt(jpkm1)
1070
1071!!org DO jk = 1, jpk
1072!!org    esigw(jk)=fsdsig( FLOAT(jk))
1073!!org    esigt(jk)=fsdsig( FLOAT(jk)+0.5)
1074!!org END DO
1075
1076      ! 2.3 Coefficients for vertical depth as the sum of e3w scale factors
1077
1078      gsi3w(1) = 0.5 * esigw(1)
1079      DO jk = 2, jpk
1080        gsi3w(jk) = gsi3w(jk-1)+ esigw(jk)
1081      END DO
1082
1083!!bug gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
1084      DO jk = 1, jpk
1085!! change the solver stat....
1086!!       zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1)
1087!!       gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigt(jk)+hift(:,:)*zcoeft)
1088         gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigt(jk)+hift(:,:)*(FLOAT(jk)-0.5)/FLOAT(jpkm1))
1089!!       zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1)
1090!!       gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigw(jk)+hift(:,:)*zcoefw)
1091!!       gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsi3w(jk)+hift(:,:)*zcoefw)
1092         gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigw(jk)+hift(:,:)*FLOAT(jk-1)/FLOAT(jpkm1))
1093         gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsi3w(jk)+hift(:,:)*FLOAT(jk-1)/FLOAT(jpkm1))
1094      END DO
1095
1096!!bug gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
1097      DO jj = 1, jpj
1098         DO ji = 1, jpi
1099            DO jk = 1, jpk
1100              e3t(ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/FLOAT(jpkm1))
1101              e3u(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/FLOAT(jpkm1))
1102              e3v(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/FLOAT(jpkm1))
1103              e3f(ji,jj,jk)=((hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/FLOAT(jpkm1))
1104
1105              e3w (ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/FLOAT(jpkm1))
1106              e3uw(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/FLOAT(jpkm1))
1107              e3vw(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/FLOAT(jpkm1))
1108            END DO
1109         END DO
1110      END DO
1111
1112! HYBRID
1113      DO jj = 1, jpj
1114         DO ji = 1, jpi
1115            DO jk = 1, jpkm1
1116               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
1117               IF( scobot(ji,jj) == 0.e0             )   mbathy(ji,jj) = 0
1118            END DO
1119         END DO
1120      END DO
1121
1122      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
1123
1124
1125      ! ===========
1126      ! Zoom domain
1127      ! ===========
1128
1129      IF( lzoom )   CALL zgr_bat_zoom
1130
1131      ! 2.4 Control print
1132
1133      IF(lwp) THEN
1134         WRITE(numout,*) 
1135         WRITE(numout,*) ' domzgr: vertical coefficients for model level'
1136         WRITE(numout,9400)
1137         WRITE(numout,9410) (jk,gsigt(jk),gsigw(jk),esigt(jk),esigw(jk),gsi3w(jk),jk=1,jpk)
1138      ENDIF
1139 9400 FORMAT(9x,'  level    gsigt      gsigw      esigt      esigw      gsi3w')
1140 9410 FORMAT(10x,i4,5f11.4)
1141
1142      IF(lwp) THEN
1143         WRITE(numout,*)
1144         WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
1145         WRITE(numout,*) ' ~~~~~~  --------------------'
1146         WRITE(numout,9420)
1147         WRITE(numout,9430) (jk,fsdept(1,1,jk),fsdepw(1,1,jk),     &
1148                             fse3t (1,1,jk),fse3w (1,1,jk),jk=1,jpk)
1149         DO jj = mj0(20), mj1(20)
1150            DO ji = mi0(20), mi1(20)
1151               WRITE(numout,*)
1152               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
1153               WRITE(numout,*) ' ~~~~~~  --------------------'
1154               WRITE(numout,9420)
1155               WRITE(numout,9430) (jk,fsdept(ji,jj,jk),fsdepw(ji,jj,jk),     &
1156                    &                 fse3t (ji,jj,jk),fse3w (ji,jj,jk),jk=1,jpk)
1157            END DO
1158         END DO
1159         DO jj = mj0(74), mj1(74)
1160            DO ji = mi0(100), mi1(100)
1161               WRITE(numout,*)
1162               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
1163               WRITE(numout,*) ' ~~~~~~  --------------------'
1164               WRITE(numout,9420)
1165               WRITE(numout,9430) (jk,fsdept(ji,jj,jk),fsdepw(ji,jj,jk),     &
1166                    &                 fse3t (ji,jj,jk),fse3w (ji,jj,jk),jk=1,jpk)
1167            END DO
1168         END DO
1169      ENDIF
1170
1171 9420 FORMAT(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')
1172 9430 FORMAT(10x,i4,4f9.2)
1173
1174!!bug gm   no more necessary?  if ! defined key_helsinki
1175      DO jk = 1, jpk
1176         DO jj = 1, jpj
1177            DO ji = 1, jpi
1178               IF( fse3w(ji,jj,jk) <= 0. .OR. fse3t(ji,jj,jk) <= 0. ) THEN
1179                  IF(lwp) THEN
1180                     WRITE(numout,*)
1181                     WRITE(numout,*) ' e r r o r :         e3w or e3t =< 0 '
1182                     WRITE(numout,*) ' =========           --------------- '
1183                     WRITE(numout,*)
1184                     WRITE(numout,*) '             point (i,j,k)= ',ji,jj,jk
1185                     WRITE(numout,*)
1186                  ENDIF
1187                  STOP 'domzgr'
1188               ENDIF
1189               IF( fsdepw(ji,jj,jk) < 0. .OR. fsdept(ji,jj,jk) < 0. ) THEN
1190                  IF(lwp) THEN
1191                     WRITE(numout,*)
1192                     WRITE(numout,*) ' e r r o r :      gdepw or gdept < 0 '
1193                     WRITE(numout,*) ' =========        ------------------ '
1194                     WRITE(numout,*)
1195                     WRITE(numout,*) '             point (i,j,k)= ', ji, jj, jk
1196                     WRITE(numout,*)
1197                  ENDIF
1198                  STOP 'domzgr'
1199               ENDIF
1200            END DO
1201         END DO
1202      END DO
1203!!bug gm    #endif
1204
1205   END SUBROUTINE zgr_sco
1206
1207#else
1208   !!----------------------------------------------------------------------
1209   !!   Default option :                                      Empty routine
1210   !!----------------------------------------------------------------------
1211   SUBROUTINE zgr_sco
1212   END SUBROUTINE zgr_sco
1213#endif
1214
1215
1216   !!======================================================================
1217END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.