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

Last change on this file since 57 was 57, checked in by opalod, 20 years ago

CT : BUGFIX031 : Running error corrected: when duplicating e3X_ps() scale factors in mpp case

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