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 @ 389

Last change on this file since 389 was 389, checked in by opalod, 18 years ago

RB:nemo_v1_update_038: first integration of Agrif :

  • configuration parameters are just integer when agrif is used
  • add call to agrif routines with key_agrif
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.5 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_partial_steps
11   !!----------------------------------------------------------------------
12   !!   'key_partial_steps' :               z-coordinate with partial steps
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_ps(i,j,k)  = fsdep(k)
33      !!                       e3w_ps(i,j,k) = dk(fsdep)(k)     = fse3(i,j,k)
34      !!              t-level: gdept_ps(i,j,k)  = fsdep(k+0.5)
35      !!                       e3t_ps(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, gdepw and e3 are positives
56      !!         - - - - - - -   gdept_ps, gdepw_ps and e3_ps 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(jpk) + e3t(jpk)
92      zmin = gdepw(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(jpk)+e3t(jpk)
123      ! i.e. the last ocean level thickness cannot exceed e3t(jpkm1)+e3t(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 (where
138      ! e3t is the reference level thickness
139      DO jk = jpkm1, 1, -1
140         zdepth = gdepw(jk) + MIN( e3zps_min, e3t(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_ps(:,:,jk) = gdept(jk)
153        gdepw_ps(:,:,jk) = gdepw(jk)
154        e3t_ps(:,:,jk) = e3t(jk)
155        e3w_ps(:,:,jk) = e3w(jk)
156     END DO
157     hdept(:,:) = gdept_ps(:,:,2 )
158     hdepw(:,:) = gdepw_ps(:,:,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(ik)
170                 ze3wp = 0.5 * e3w(ik) * ( 1. + ( ze3tp/e3t(ik) ) )
171                 e3t_ps(ji,jj,ik  ) = ze3tp
172                 e3t_ps(ji,jj,ik+1) = ze3tp
173                 e3w_ps(ji,jj,ik  ) = ze3wp
174                 e3w_ps(ji,jj,ik+1) = ze3tp
175                 gdepw_ps(ji,jj,ik+1) = zdepwp
176                 gdept_ps(ji,jj,ik  ) = gdept(ik-1) + ze3wp
177                 gdept_ps(ji,jj,ik+1) = gdept_ps(ji,jj,ik) + ze3tp
178                 ! standard case
179              ELSE
180!!alex
181                 IF( bathy(ji,jj) <= gdepw(ik+1) ) THEN
182                    gdepw_ps(ji,jj,ik+1) = bathy(ji,jj)
183                 ELSE
184!!alex ctl          write(*,*) 'zps',ji,jj,'bathy', bathy(ji,jj), 'depw_ps ',gdepw(ik+1)
185                    gdepw_ps(ji,jj,ik+1) = gdepw(ik+1)
186                 ENDIF
187!!Alex
188!!Alex           gdepw_ps(ji,jj,ik+1) = bathy(ji,jj)
189                 gdept_ps(ji,jj,ik  ) =  gdepw(ik) + ( gdepw_ps(ji,jj,ik+1) - gdepw(ik))   &
190                                      * ((gdept   (      ik  ) - gdepw(ik))   &
191                                      / ( gdepw   (      ik+1) - gdepw(ik)))
192                 e3t_ps(ji,jj,ik) = e3t(ik) * ( gdepw_ps(ji,jj,ik+1) - gdepw(ik))   &
193                                  /( gdepw   (      ik+1) - gdepw(ik))
194                 e3w_ps(ji,jj,ik) = 0.5 *( gdepw_ps(ji,jj,ik+1) + gdepw(ik+1)-2.*gdepw(ik))   &
195                                  *( e3w(ik) / ( gdepw(ik+1) - gdepw(ik)))
196                 !       ... on ik+1
197                 e3w_ps(ji,jj,ik+1) = e3t_ps(ji,jj,ik)
198                 e3t_ps(ji,jj,ik+1) = e3t_ps(ji,jj,ik)
199                 gdept_ps(ji,jj,ik+1) = gdept_ps(ji,jj,ik) + e3t_ps  (ji,jj,ik)
200              ENDIF
201           ENDIF
202        END DO
203     END DO
204
205     it = 0
206     DO jj = 1, jpj
207        DO ji = 1, jpi
208           ik = mbathy(ji,jj)
209           ! ocean point only
210           IF( ik > 0 ) THEN
211              ! bathymetry output
212              hdept(ji,jj) = gdept_ps(ji,jj,ik  )
213              hdepw(ji,jj) = gdepw_ps(ji,jj,ik+1)
214              e3tp (ji,jj) = e3t_ps(ji,jj,ik  )
215              e3wp (ji,jj) = e3w_ps(ji,jj,ik  )
216              ! test
217              zdiff= gdepw_ps(ji,jj,ik+1) - gdept_ps(ji,jj,ik  )
218              IF( zdiff <= 0. .AND. lwp ) THEN
219                 it=it+1
220                 WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj
221                 WRITE(numout,*) ' bathy = ', bathy(ji,jj)
222                 WRITE(numout,*) ' gdept_ps= ', gdept_ps(ji,jj,ik), ' gdepw_ps= ', gdepw_ps(ji,jj,ik+1),   &
223                                 ' zdiff   = ', zdiff
224                 WRITE(numout,*) ' e3tp    = ', e3t_ps(ji,jj,ik  ), ' e3wp    = ', e3w_ps(ji,jj,ik  )
225              ENDIF
226           ENDIF
227        END DO
228      END DO
229
230      ! Scale factors and depth at U-, V-, UW and VW-points
231
232      ! initialisation to z-scale factors
233      DO jk = 1, jpk
234         e3u_ps (:,:,jk)  = e3t(jk)
235         e3v_ps (:,:,jk)  = e3t(jk)
236         e3uw_ps(:,:,jk)  = e3w(jk)
237         e3vw_ps(:,:,jk)  = e3w(jk)
238      END DO
239
240     ! Computed as the minimum of neighbooring scale factors
241     DO jk = 1,jpk
242        DO jj = 1, jpjm1
243           DO ji = 1, fs_jpim1   ! vector opt.
244              e3u_ps (ji,jj,jk) = MIN( e3t_ps(ji,jj,jk), e3t_ps(ji+1,jj,jk))
245              e3v_ps (ji,jj,jk) = MIN( e3t_ps(ji,jj,jk), e3t_ps(ji,jj+1,jk))
246              e3uw_ps(ji,jj,jk) = MIN( e3w_ps(ji,jj,jk), e3w_ps(ji+1,jj,jk) )
247              e3vw_ps(ji,jj,jk) = MIN( e3w_ps(ji,jj,jk), e3w_ps(ji,jj+1,jk) )
248           END DO
249        END DO
250     END DO
251     
252     ! Boundary conditions
253     CALL lbc_lnk( e3u_ps , 'U', 1. )   ;   CALL lbc_lnk( e3uw_ps, 'U', 1. )
254     CALL lbc_lnk( e3v_ps , 'V', 1. )   ;   CALL lbc_lnk( e3vw_ps, 'V', 1. )
255     
256     ! set to z-scale factor if zero (i.e. along closed boundaries)
257     DO jk = 1, jpk
258        DO jj = 1, jpj
259           DO ji = 1, jpi
260              IF( e3u_ps (ji,jj,jk) == 0.e0 ) e3u_ps (ji,jj,jk) = e3t(jk)
261              IF( e3v_ps (ji,jj,jk) == 0.e0 ) e3v_ps (ji,jj,jk) = e3t(jk)
262              IF( e3uw_ps(ji,jj,jk) == 0.e0 ) e3uw_ps(ji,jj,jk) = e3w(jk)
263              IF( e3vw_ps(ji,jj,jk) == 0.e0 ) e3vw_ps(ji,jj,jk) = e3w(jk)
264           END DO
265        END DO
266     END DO
267     
268     ! Scale factor at F-point
269     
270     ! initialisation to z-scale factors
271     DO jk = 1, jpk
272        e3f_ps (:,:,jk) = e3t(jk)
273     END DO
274     
275     ! Computed as the minimum of neighbooring V-scale factors
276     DO jk = 1, jpk
277        DO jj = 1, jpjm1
278           DO ji = 1, fs_jpim1   ! vector opt.
279              e3f_ps(ji,jj,jk) = MIN( e3v_ps(ji,jj,jk), e3v_ps(ji+1,jj,jk) )
280           END DO
281        END DO
282     END DO
283     ! Boundary conditions
284     CALL lbc_lnk( e3f_ps, 'F', 1. )
285     
286     ! set to z-scale factor if zero (i.e. along closed boundaries)
287     DO jk = 1, jpk
288        DO jj = 1, jpj
289           DO ji = 1, jpi
290              IF( e3f_ps(ji,jj,jk) == 0.e0 ) e3f_ps(ji,jj,jk) = e3t(jk)
291           END DO
292        END DO
293     END DO
294     ! we duplicate factor scales for jj = 1 and jj = 2
295     e3t_ps(:,mj0(1),:) = e3t_ps(:,mj0(2),:)
296     e3w_ps(:,mj0(1),:) = e3w_ps(:,mj0(2),:)
297     e3u_ps(:,mj0(1),:) = e3u_ps(:,mj0(2),:)
298     e3v_ps(:,mj0(1),:) = e3v_ps(:,mj0(2),:)
299     e3f_ps(:,mj0(1),:) = e3f_ps(:,mj0(2),:)
300
301
302     
303     ! Compute gdep3w (vertical sum of e3w)
304     
305     gdep3w   (:,:,1) = 0.5 * e3w_ps (:,:,1)
306     DO jk = 2, jpk
307        gdep3w   (:,:,jk) = gdep3w   (:,:,jk-1) + e3w_ps (:,:,jk)
308     END DO
309         
310     ! Control print
311 9600 FORMAT(9x,' level   gdept    gdepw     e3t      e3w   ')
312 9610 FORMAT(10x,i4,4f9.2)
313      IF(lwp .AND. ll_print) THEN
314         DO jj = 1,jpj
315            DO ji = 1, jpi
316               ik = MAX(mbathy(ji,jj),1)
317               zprt(ji,jj) = e3t_ps(ji,jj,ik)
318            END DO
319         END DO
320         WRITE(numout,*)
321         WRITE(numout,*) 'domzgr e3t(mbathy)'
322         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
323         DO jj = 1,jpj
324            DO ji = 1, jpi
325               ik = MAX(mbathy(ji,jj),1)
326               zprt(ji,jj) = e3w_ps(ji,jj,ik)
327            END DO
328         END DO
329         WRITE(numout,*)
330         WRITE(numout,*) 'domzgr e3w(mbathy)'
331         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
332         DO jj = 1,jpj
333            DO ji = 1, jpi
334               ik = MAX(mbathy(ji,jj),1)
335               zprt(ji,jj) = e3u_ps(ji,jj,ik)
336            END DO
337         END DO
338
339         WRITE(numout,*)
340         WRITE(numout,*) 'domzgr e3u(mbathy)'
341         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
342         DO jj = 1,jpj
343            DO ji = 1, jpi
344               ik = MAX(mbathy(ji,jj),1)
345               zprt(ji,jj) = e3v_ps(ji,jj,ik)
346            END DO
347         END DO
348         WRITE(numout,*)
349         WRITE(numout,*) 'domzgr e3v(mbathy)'
350         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
351         DO jj = 1,jpj
352            DO ji = 1, jpi
353               ik = MAX(mbathy(ji,jj),1)
354               zprt(ji,jj) = e3f_ps(ji,jj,ik)
355            END DO
356         END DO
357
358         WRITE(numout,*)
359         WRITE(numout,*) 'domzgr e3f(mbathy)'
360         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
361         DO jj = 1,jpj
362            DO ji = 1, jpi
363               ik =  MAX(mbathy(ji,jj),1)
364               zprt(ji,jj) = gdep3w(ji,jj,ik)
365            END DO
366         END DO
367         WRITE(numout,*)
368         WRITE(numout,*) 'domzgr gdep3w(mbathy)'
369         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
370
371      ENDIF 
372
373
374      DO jk = 1,jpk
375         DO jj = 1, jpj
376            DO ji = 1, jpi
377               IF( e3w_ps(ji,jj,jk) <= 0. .or. e3t_ps(ji,jj,jk) <= 0. ) THEN
378                  IF(lwp) THEN
379                     WRITE(numout,*) ' e r r o r :         e3w or e3t =< 0 '
380                     WRITE(numout,*) ' =========           --------------- '
381                     WRITE(numout,*)
382                  ENDIF
383                  STOP 'domzgr.psteps'
384               ENDIF
385               IF( gdepw_ps(ji,jj,jk) < 0. .or. gdept_ps(ji,jj,jk) < 0. ) THEN
386                  IF(lwp) THEN
387                     WRITE(numout,*) ' e r r o r :      gdepw or gdept < 0 '
388                     WRITE(numout,*) ' =========        ------------------ '
389                     WRITE(numout,*)
390                  ENDIF
391                  STOP 'domzgr.psteps'
392               ENDIF
393            END DO
394         END DO
395      END DO 
396
397   IF(lwp) THEN
398      WRITE(numout,*) ' e3t lev 21 '
399      CALL prihre(e3t_ps(:,:,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
400      WRITE(numout,*) ' e3w lev 21  '
401      CALL prihre(e3w_ps(:,:,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
402      WRITE(numout,*) ' e3u lev 21  '
403      CALL prihre(e3u_ps(:,:,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
404      WRITE(numout,*) ' e3v lev 21  '
405      CALL prihre(e3v_ps(:,:,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
406      WRITE(numout,*) ' e3f lev 21  '
407      CALL prihre(e3f_ps(:,:,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
408      WRITE(numout,*) ' e3t lev 22 '
409      CALL prihre(e3t_ps(:,:,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
410      WRITE(numout,*) ' e3w lev 22  '
411      CALL prihre(e3w_ps(:,:,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
412      WRITE(numout,*) ' e3u lev 22  '
413      CALL prihre(e3u_ps(:,:,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
414      WRITE(numout,*) ' e3v lev 22  '
415      CALL prihre(e3v_ps(:,:,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
416      WRITE(numout,*) ' e3f lev 22  '
417      CALL prihre(e3f_ps(:,:,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
418   ENDIF
419
420      ! ===========
421      ! Zoom domain
422      ! ===========
423
424      IF( lzoom )   CALL zgr_bat_zoom
425
426      ! ================
427      ! Bathymetry check
428      ! ================
429
430      CALL zgr_bat_ctl
431
432
433   END SUBROUTINE zgr_zps
434
435#else
436   !!----------------------------------------------------------------------
437   !!   Default option :                                      Empty routine
438   !!----------------------------------------------------------------------
439   SUBROUTINE zgr_zps
440   END SUBROUTINE zgr_zps
441#endif
Note: See TracBrowser for help on using the repository browser.