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_zps.h90 in trunk/NEMO/OPA_SRC/DOM – NEMO

source: trunk/NEMO/OPA_SRC/DOM/domzgr_zps.h90 @ 719

Last change on this file since 719 was 719, checked in by ctlod, 17 years ago

get back to the nemo_v2_3 version for trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.4 KB
Line 
1   !!----------------------------------------------------------------------
2   !!                      ***  domzgr_zps.h90  ***
3   !!----------------------------------------------------------------------
4   !!----------------------------------------------------------------------
5   !!  OPA 9.0 , LOCEAN-IPSL (2005)
6   !! $Header$
7   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
8   !!----------------------------------------------------------------------
9
10#if ! defined key_zco
11   !!----------------------------------------------------------------------
12   !!   NOT 'key_zco' :                                    3D gdep. and e3.
13   !!----------------------------------------------------------------------
14
15   SUBROUTINE zgr_zps
16      !!----------------------------------------------------------------------
17      !!                  ***  ROUTINE zgr_zps  ***
18      !!                     
19      !! ** Purpose :   the depth and vertical scale factor in partial step
20      !!      z-coordinate case
21      !!
22      !! ** Method  :   Partial steps : computes the 3D vertical scale factors
23      !!      of T-, U-, V-, W-, UW-, VW and F-points that are associated with
24      !!      a partial step representation of bottom topography.
25      !!
26      !!        The reference depth of model levels is defined from an analytical
27      !!      function the derivative of which gives the reference vertical
28      !!      scale factors.
29      !!        From  depth and scale factors reference, we compute there new value
30      !!      with partial steps  on 3d arrays ( i, j, k ).
31      !!
32      !!              w-level: gdepw(i,j,k)  = fsdep(k)
33      !!                       e3w(i,j,k) = dk(fsdep)(k)     = fse3(i,j,k)
34      !!              t-level: gdept(i,j,k)  = fsdep(k+0.5)
35      !!                       e3t(i,j,k) = dk(fsdep)(k+0.5) = fse3(i,j,k+0.5)
36      !!
37      !!        With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
38      !!      we find the mbathy index of the depth at each grid point.
39      !!      This leads us to three cases:
40      !!
41      !!              - bathy = 0 => mbathy = 0
42      !!              - 1 < mbathy < jpkm1   
43      !!              - bathy > gdepw(jpk) => mbathy = jpkm1 
44      !!
45      !!        Then, for each case, we find the new depth at t- and w- levels
46      !!      and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
47      !!      and f-points.
48      !!
49      !!        This routine is given as an example, it must be modified
50      !!      following the user s desiderata. nevertheless, the output as
51      !!      well as the way to compute the model levels and scale factors
52      !!      must be respected in order to insure second order accuracy
53      !!      schemes.
54      !!
55      !!         c a u t i o n : gdept_0, gdepw_0 and e3._0 are positives
56      !!         - - - - - - -   gdept, gdepw and e3. are positives
57      !!     
58      !!  Reference :
59      !!     Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
60      !!
61      !!  History :
62      !!    8.5  !  02-09 (A. Bozec, G. Madec)  F90: Free form and module
63      !!    9.0  !  02-09 (A. de Miranda)  rigid-lid + islands
64      !!----------------------------------------------------------------------
65      !! * Local declarations
66      INTEGER  ::   ji, jj, jk    ! dummy loop indices
67      INTEGER  ::   &
68         ik, it            ! temporary integers
69     
70      REAL(wp) ::   & 
71         ze3tp, ze3wp,    &  ! Last ocean level thickness at T- and W-points
72         zdepwp,          &  ! Ajusted ocean depth to avoid too small e3t
73         zdepth,          &  !    "         "
74         zmax, zmin,      &  ! Maximum and minimum depth
75         zdiff               ! temporary scalar
76
77      REAL(wp), DIMENSION(jpi,jpj) ::   &
78         zprt  !    "           "
79
80      LOGICAL ::  ll_print                          ! Allow  control print for debugging
81
82      !!---------------------------------------------------------------------
83      !! OPA8.5, LODYC-IPSL (2002)
84      !!---------------------------------------------------------------------
85     
86      ! Local variable for debugging
87      ll_print=.FALSE.
88!!!   ll_print=.TRUE.
89     
90      ! Initialization of constant
91      zmax = gdepw_0(jpk) + e3t_0(jpk)
92      zmin = gdepw_0(4)
93     
94      ! Ocean depth
95      IF(lwp .AND. ll_print) THEN
96         WRITE(numout,*)
97         WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)'
98         CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
99      ENDIF
100
101      IF(lwp) WRITE(numout,*)
102      IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps'
103      IF(lwp) WRITE(numout,*) '    ~~~~~~~ '
104      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used'
105
106
107      ! bathymetry in level (from bathy_meter)
108      ! ===================
109
110      ! initialize mbathy to the maximum ocean level available
111      mbathy(:,:) = jpkm1
112
113      ! storage of land and island's number (zera and negative values) in mbathy
114      DO jj = 1, jpj
115         DO ji= 1, jpi
116            IF( bathy(ji,jj) <= 0. )   mbathy(ji,jj) = INT( bathy(ji,jj) )
117         END DO
118      END DO
119
120      ! bounded value of bathy
121      ! minimum depth == 3 levels
122      ! maximum depth == gdepw_0(jpk)+e3t_0(jpk)
123      ! i.e. the last ocean level thickness cannot exceed e3t_0(jpkm1)+e3t_0(jpk)
124      DO jj = 1, jpj
125         DO ji= 1, jpi
126            IF( bathy(ji,jj) <= 0. ) THEN
127               bathy(ji,jj) = 0.e0
128            ELSE
129               bathy(ji,jj) = MAX( bathy(ji,jj), zmin )
130               bathy(ji,jj) = MIN( bathy(ji,jj), zmax )
131            ENDIF
132         END DO
133      END DO
134
135      ! Compute mbathy for ocean points (i.e. the number of ocean levels)
136      ! find the number of ocean levels such that the last level thickness
137      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_0 (where
138      ! e3t_0 is the reference level thickness
139      DO jk = jpkm1, 1, -1
140         zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat )
141         DO jj = 1, jpj
142            DO ji = 1, jpi
143               IF( 0. < bathy(ji,jj) .AND. bathy(ji,jj) <= zdepth )   mbathy(ji,jj) = jk-1
144            END DO
145         END DO
146      END DO
147
148     ! Scale factors and depth at T- and W-points
149     
150     ! intitialization to the reference z-coordinate
151     DO jk = 1, jpk
152        gdept(:,:,jk) = gdept_0(jk)
153        gdepw(:,:,jk) = gdepw_0(jk)
154        e3t  (:,:,jk) = e3t_0  (jk)
155        e3w  (:,:,jk) = e3w_0  (jk)
156     END DO
157     hdept(:,:) = gdept(:,:,2 )
158     hdepw(:,:) = gdepw(:,:,3 )     
159     
160     !
161     DO jj = 1, jpj
162        DO ji = 1, jpi
163           ik = mbathy(ji,jj)
164           ! ocean point only
165           IF( ik > 0 ) THEN
166              ! max ocean level case
167              IF( ik == jpkm1 ) THEN
168                 zdepwp = bathy(ji,jj)
169                 ze3tp  = bathy(ji,jj) - gdepw_0(ik)
170                 ze3wp = 0.5 * e3w_0(ik) * ( 1. + ( ze3tp/e3t_0(ik) ) )
171                 e3t(ji,jj,ik  ) = ze3tp
172                 e3t(ji,jj,ik+1) = ze3tp
173                 e3w(ji,jj,ik  ) = ze3wp
174                 e3w(ji,jj,ik+1) = ze3tp
175                 gdepw(ji,jj,ik+1) = zdepwp
176                 gdept(ji,jj,ik  ) = gdept_0(ik-1) + ze3wp
177                 gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp
178                 ! standard case
179              ELSE
180!!alex
181                 IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN
182                    gdepw(ji,jj,ik+1) = bathy(ji,jj)
183                 ELSE
184!!alex ctl          write(*,*) 'zps',ji,jj,'bathy', bathy(ji,jj), 'depw ',gdepw_0(ik+1)
185                    gdepw(ji,jj,ik+1) = gdepw_0(ik+1)
186                 ENDIF
187!!Alex
188!!Alex           gdepw(ji,jj,ik+1) = bathy(ji,jj)
189!!bug gm  verifier les gdepw_0
190                 gdept(ji,jj,ik  ) =  gdepw_0(ik) + ( gdepw(ji,jj,ik+1) - gdepw_0(ik) )   &
191                    &                             * ((gdept_0(      ik  ) - gdepw_0(ik) )   &
192                    &                             / ( gdepw_0(      ik+1) - gdepw_0(ik) ))
193                 e3t(ji,jj,ik) = e3t_0(ik) * ( gdepw  (ji,jj,ik+1) - gdepw_0(ik) )   &
194                    &                      / ( gdepw_0(      ik+1) - gdepw_0(ik) )
195                 e3w(ji,jj,ik) = 0.5 * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2.*gdepw_0(ik) )   &
196                    &                * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) )
197                 !       ... on ik+1
198                 e3w(ji,jj,ik+1) = e3t(ji,jj,ik)
199                 e3t(ji,jj,ik+1) = e3t(ji,jj,ik)
200                 gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)
201              ENDIF
202           ENDIF
203        END DO
204     END DO
205
206     it = 0
207     DO jj = 1, jpj
208        DO ji = 1, jpi
209           ik = mbathy(ji,jj)
210           ! ocean point only
211           IF( ik > 0 ) THEN
212              ! bathymetry output
213              hdept(ji,jj) = gdept(ji,jj,ik  )
214              hdepw(ji,jj) = gdepw(ji,jj,ik+1)
215              e3tp (ji,jj) = e3t(ji,jj,ik  )
216              e3wp (ji,jj) = e3w(ji,jj,ik  )
217              ! test
218              zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik  )
219              IF( zdiff <= 0. .AND. lwp ) THEN
220                 it=it+1
221                 WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj
222                 WRITE(numout,*) ' bathy = ', bathy(ji,jj)
223                 WRITE(numout,*) ' gdept= ', gdept(ji,jj,ik), ' gdepw= ', gdepw(ji,jj,ik+1),   &
224                                 ' zdiff   = ', zdiff
225                 WRITE(numout,*) ' e3tp    = ', e3t(ji,jj,ik  ), ' e3wp    = ', e3w(ji,jj,ik  )
226              ENDIF
227           ENDIF
228        END DO
229      END DO
230
231      ! Scale factors and depth at U-, V-, UW and VW-points
232
233      ! initialisation to z-scale factors
234      DO jk = 1, jpk
235         e3u (:,:,jk)  = e3t_0(jk)
236         e3v (:,:,jk)  = e3t_0(jk)
237         e3uw(:,:,jk)  = e3w_0(jk)
238         e3vw(:,:,jk)  = e3w_0(jk)
239      END DO
240
241     ! Computed as the minimum of neighbooring scale factors
242     DO jk = 1,jpk
243        DO jj = 1, jpjm1
244           DO ji = 1, fs_jpim1   ! vector opt.
245              e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk))
246              e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk))
247              e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) )
248              e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) )
249           END DO
250        END DO
251     END DO
252     
253     ! Boundary conditions
254     CALL lbc_lnk( e3u , 'U', 1. )   ;   CALL lbc_lnk( e3uw, 'U', 1. )
255     CALL lbc_lnk( e3v , 'V', 1. )   ;   CALL lbc_lnk( e3vw, 'V', 1. )
256     
257     ! set to z-scale factor if zero (i.e. along closed boundaries)
258     DO jk = 1, jpk
259        DO jj = 1, jpj
260           DO ji = 1, jpi
261              IF( e3u (ji,jj,jk) == 0.e0 ) e3u (ji,jj,jk) = e3t_0(jk)
262              IF( e3v (ji,jj,jk) == 0.e0 ) e3v (ji,jj,jk) = e3t_0(jk)
263              IF( e3uw(ji,jj,jk) == 0.e0 ) e3uw(ji,jj,jk) = e3w_0(jk)
264              IF( e3vw(ji,jj,jk) == 0.e0 ) e3vw(ji,jj,jk) = e3w_0(jk)
265           END DO
266        END DO
267     END DO
268     
269     ! Scale factor at F-point
270     
271     ! initialisation to z-scale factors
272     DO jk = 1, jpk
273        e3f (:,:,jk) = e3t_0(jk)
274     END DO
275     
276     ! Computed as the minimum of neighbooring V-scale factors
277     DO jk = 1, jpk
278        DO jj = 1, jpjm1
279           DO ji = 1, fs_jpim1   ! vector opt.
280              e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) )
281           END DO
282        END DO
283     END DO
284     ! Boundary conditions
285     CALL lbc_lnk( e3f, 'F', 1. )
286     
287     ! set to z-scale factor if zero (i.e. along closed boundaries)
288     DO jk = 1, jpk
289        DO jj = 1, jpj
290           DO ji = 1, jpi
291              IF( e3f(ji,jj,jk) == 0.e0 ) e3f(ji,jj,jk) = e3t_0(jk)
292           END DO
293        END DO
294     END DO
295!!bug ? gm:  must be a do loop with mj0,mj1
296
297     ! we duplicate factor scales for jj = 1 and jj = 2
298     e3t(:,mj0(1),:) = e3t(:,mj0(2),:)
299     e3w(:,mj0(1),:) = e3w(:,mj0(2),:)
300     e3u(:,mj0(1),:) = e3u(:,mj0(2),:)
301     e3v(:,mj0(1),:) = e3v(:,mj0(2),:)
302     e3f(:,mj0(1),:) = e3f(:,mj0(2),:)
303
304
305     
306     ! Compute gdep3w (vertical sum of e3w)
307     
308     gdep3w   (:,:,1) = 0.5 * e3w (:,:,1)
309     DO jk = 2, jpk
310        gdep3w   (:,:,jk) = gdep3w   (:,:,jk-1) + e3w (:,:,jk)
311     END DO
312         
313     ! Control print
314 9600 FORMAT(9x,' level   gdept    gdepw     e3t      e3w   ')
315 9610 FORMAT(10x,i4,4f9.2)
316      IF(lwp .AND. ll_print) THEN
317         DO jj = 1,jpj
318            DO ji = 1, jpi
319               ik = MAX(mbathy(ji,jj),1)
320               zprt(ji,jj) = e3t(ji,jj,ik)
321            END DO
322         END DO
323         WRITE(numout,*)
324         WRITE(numout,*) 'domzgr e3t(mbathy)'
325         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
326         DO jj = 1,jpj
327            DO ji = 1, jpi
328               ik = MAX(mbathy(ji,jj),1)
329               zprt(ji,jj) = e3w(ji,jj,ik)
330            END DO
331         END DO
332         WRITE(numout,*)
333         WRITE(numout,*) 'domzgr e3w(mbathy)'
334         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
335         DO jj = 1,jpj
336            DO ji = 1, jpi
337               ik = MAX(mbathy(ji,jj),1)
338               zprt(ji,jj) = e3u(ji,jj,ik)
339            END DO
340         END DO
341
342         WRITE(numout,*)
343         WRITE(numout,*) 'domzgr e3u(mbathy)'
344         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
345         DO jj = 1,jpj
346            DO ji = 1, jpi
347               ik = MAX(mbathy(ji,jj),1)
348               zprt(ji,jj) = e3v(ji,jj,ik)
349            END DO
350         END DO
351         WRITE(numout,*)
352         WRITE(numout,*) 'domzgr e3v(mbathy)'
353         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
354         DO jj = 1,jpj
355            DO ji = 1, jpi
356               ik = MAX(mbathy(ji,jj),1)
357               zprt(ji,jj) = e3f(ji,jj,ik)
358            END DO
359         END DO
360
361         WRITE(numout,*)
362         WRITE(numout,*) 'domzgr e3f(mbathy)'
363         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
364         DO jj = 1,jpj
365            DO ji = 1, jpi
366               ik =  MAX(mbathy(ji,jj),1)
367               zprt(ji,jj) = gdep3w(ji,jj,ik)
368            END DO
369         END DO
370         WRITE(numout,*)
371         WRITE(numout,*) 'domzgr gdep3w(mbathy)'
372         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
373
374      ENDIF 
375
376
377      DO jk = 1,jpk
378         DO jj = 1, jpj
379            DO ji = 1, jpi
380               IF( e3w(ji,jj,jk) <= 0. .or. e3t(ji,jj,jk) <= 0. ) THEN
381                  IF(lwp) THEN
382                     WRITE(numout,*) ' e r r o r :         e3w or e3t =< 0 '
383                     WRITE(numout,*) ' =========           --------------- '
384                     WRITE(numout,*)
385                  ENDIF
386                  STOP 'domzgr.psteps'
387               ENDIF
388               IF( gdepw(ji,jj,jk) < 0. .or. gdept(ji,jj,jk) < 0. ) THEN
389                  IF(lwp) THEN
390                     WRITE(numout,*) ' e r r o r :      gdepw or gdept < 0 '
391                     WRITE(numout,*) ' =========        ------------------ '
392                     WRITE(numout,*)
393                  ENDIF
394                  STOP 'domzgr.psteps'
395               ENDIF
396            END DO
397         END DO
398      END DO 
399
400   IF(lwp) THEN
401      WRITE(numout,*) ' e3t lev 21 '
402      CALL prihre(e3t(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
403      WRITE(numout,*) ' e3w lev 21  '
404      CALL prihre(e3w(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
405      WRITE(numout,*) ' e3u lev 21  '
406      CALL prihre(e3u(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
407      WRITE(numout,*) ' e3v lev 21  '
408      CALL prihre(e3v(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
409      WRITE(numout,*) ' e3f lev 21  '
410      CALL prihre(e3f(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
411      WRITE(numout,*) ' e3t lev 22 '
412      CALL prihre(e3t(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
413      WRITE(numout,*) ' e3w lev 22  '
414      CALL prihre(e3w(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
415      WRITE(numout,*) ' e3u lev 22  '
416      CALL prihre(e3u(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
417      WRITE(numout,*) ' e3v lev 22  '
418      CALL prihre(e3v(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
419      WRITE(numout,*) ' e3f lev 22  '
420      CALL prihre(e3f(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
421   ENDIF
422
423      ! ===========
424      ! Zoom domain
425      ! ===========
426
427      IF( lzoom )   CALL zgr_bat_zoom
428
429      ! ================
430      ! Bathymetry check
431      ! ================
432
433      CALL zgr_bat_ctl
434
435   END SUBROUTINE zgr_zps
436
437#else
438   !!----------------------------------------------------------------------
439   !!   Default option :                                      Empty routine
440   !!----------------------------------------------------------------------
441   SUBROUTINE zgr_zps
442   END SUBROUTINE zgr_zps
443#endif
Note: See TracBrowser for help on using the repository browser.