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

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

Initial revision

  • 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(:,1,:) = e3t_ps(:,2,:)
291     e3w_ps(:,1,:) = e3w_ps(:,2,:)
292     e3u_ps(:,1,:) = e3u_ps(:,2,:)
293     e3v_ps(:,1,:) = e3v_ps(:,2,:)
294     e3f_ps(:,1,:) = e3f_ps(:,2,:)
295
296     
297     ! Compute gdep3w (vertical sum of e3w)
298     
299     gdep3w   (:,:,1) = 0.5 * e3w_ps (:,:,1)
300     DO jk = 2, jpk
301        gdep3w   (:,:,jk) = gdep3w   (:,:,jk-1) + e3w_ps (:,:,jk)
302     END DO
303         
304     ! Control print
305 9600 FORMAT(9x,' level   gdept    gdepw     e3t      e3w   ')
306 9610 FORMAT(10x,i4,4f9.2)
307      IF(lwp .AND. ll_print) THEN
308         DO jj = 1,jpj
309            DO ji = 1, jpi
310               ik = MAX(mbathy(ji,jj),1)
311               zprt(ji,jj) = e3t_ps(ji,jj,ik)
312            END DO
313         END DO
314         WRITE(numout,*)
315         WRITE(numout,*) 'domzgr e3t(mbathy)'
316         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
317         DO jj = 1,jpj
318            DO ji = 1, jpi
319               ik = MAX(mbathy(ji,jj),1)
320               zprt(ji,jj) = e3w_ps(ji,jj,ik)
321            END DO
322         END DO
323         WRITE(numout,*)
324         WRITE(numout,*) 'domzgr e3w(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) = e3u_ps(ji,jj,ik)
330            END DO
331         END DO
332
333         WRITE(numout,*)
334         WRITE(numout,*) 'domzgr e3u(mbathy)'
335         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
336         DO jj = 1,jpj
337            DO ji = 1, jpi
338               ik = MAX(mbathy(ji,jj),1)
339               zprt(ji,jj) = e3v_ps(ji,jj,ik)
340            END DO
341         END DO
342         WRITE(numout,*)
343         WRITE(numout,*) 'domzgr e3v(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) = e3f_ps(ji,jj,ik)
349            END DO
350         END DO
351
352         WRITE(numout,*)
353         WRITE(numout,*) 'domzgr e3f(mbathy)'
354         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
355         DO jj = 1,jpj
356            DO ji = 1, jpi
357               ik =  MAX(mbathy(ji,jj),1)
358               zprt(ji,jj) = gdep3w(ji,jj,ik)
359            END DO
360         END DO
361         WRITE(numout,*)
362         WRITE(numout,*) 'domzgr gdep3w(mbathy)'
363         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
364
365      ENDIF 
366
367
368      DO jk = 1,jpk
369         DO jj = 1, jpj
370            DO ji = 1, jpi
371               IF( e3w_ps(ji,jj,jk) <= 0. .or. e3t_ps(ji,jj,jk) <= 0. ) THEN
372                  IF(lwp) THEN
373                     WRITE(numout,*) ' e r r o r :         e3w or e3t =< 0 '
374                     WRITE(numout,*) ' =========           --------------- '
375                     WRITE(numout,*)
376                  ENDIF
377                  STOP 'domzgr.psteps'
378               ENDIF
379               IF( gdepw_ps(ji,jj,jk) < 0. .or. gdept_ps(ji,jj,jk) < 0. ) THEN
380                  IF(lwp) THEN
381                     WRITE(numout,*) ' e r r o r :      gdepw or gdept < 0 '
382                     WRITE(numout,*) ' =========        ------------------ '
383                     WRITE(numout,*)
384                  ENDIF
385                  STOP 'domzgr.psteps'
386               ENDIF
387            END DO
388         END DO
389      END DO 
390
391   IF(lwp) THEN
392      WRITE(numout,*) ' e3t lev 21 '
393      CALL prihre(e3t_ps(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
394      WRITE(numout,*) ' e3w lev 21  '
395      CALL prihre(e3w_ps(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
396      WRITE(numout,*) ' e3u lev 21  '
397      CALL prihre(e3u_ps(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
398      WRITE(numout,*) ' e3v lev 21  '
399      CALL prihre(e3v_ps(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
400      WRITE(numout,*) ' e3f lev 21  '
401      CALL prihre(e3f_ps(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
402      WRITE(numout,*) ' e3t lev 22 '
403      CALL prihre(e3t_ps(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
404      WRITE(numout,*) ' e3w lev 22  '
405      CALL prihre(e3w_ps(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
406      WRITE(numout,*) ' e3u lev 22  '
407      CALL prihre(e3u_ps(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
408      WRITE(numout,*) ' e3v lev 22  '
409      CALL prihre(e3v_ps(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
410      WRITE(numout,*) ' e3f lev 22  '
411      CALL prihre(e3f_ps(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
412   ENDIF
413
414      ! ================
415      ! Bathymetry check
416      ! ================
417
418      CALL zgr_bat_ctl
419
420
421   END SUBROUTINE zgr_zps
422
423#else
424   !!----------------------------------------------------------------------
425   !!   Default option :                                      Empty routine
426   !!----------------------------------------------------------------------
427   SUBROUTINE zgr_zps
428   END SUBROUTINE zgr_zps
429#endif
Note: See TracBrowser for help on using the repository browser.