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.
usrdef_zgr.F90 in branches/2016/dev_merge_2016/NEMOGCM/CONFIG/WAD_TEST_CASES/MY_SRC – NEMO

source: branches/2016/dev_merge_2016/NEMOGCM/CONFIG/WAD_TEST_CASES/MY_SRC/usrdef_zgr.F90 @ 7467

Last change on this file since 7467 was 7467, checked in by acc, 8 years ago

Changes to WAD_TEST_CASES/MY_SRC files and namelist_cfg to successfully run with the merge branch. Some changes will be moved into the main source once verified
File size: 21.3 KB
Line 
1MODULE usrdef_zgr
2   !!======================================================================
3   !!                   ***  MODULE  usrdef_zgr  ***
4   !!
5   !!                   ===  WAD_TEST_CASES case  ===
6   !!
7   !! Ocean domain : user defined vertical coordinate system
8   !!======================================================================
9   !! History :  4.0  ! 2016-08  (G. Madec, S. Flavoni)  Original code
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   usr_def_zgr   : user defined vertical coordinate system (required)
14   !!       zgr_z1d   : reference 1D z-coordinate
15   !!---------------------------------------------------------------------
16   USE oce            ! ocean variables
17   USE dom_oce ,  ONLY: ln_zco, ln_zps, ln_sco   ! ocean space and time domain
18   USE dom_oce ,  ONLY: mi0, mi1, nimpp, njmpp,  &
19                      & mj0, mj1, glamt, gphit   ! ocean space and time domain
20   USE usrdef_nam     ! User defined : namelist variables
21   USE wet_dry ,  ONLY: rn_wdmin1, rn_wdmin2, rn_wdld, ht_wd
22   !
23   USE in_out_manager ! I/O manager
24   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
25   USE lib_mpp        ! distributed memory computing library
26   USE wrk_nemo       ! Memory allocation
27   USE timing         ! Timing
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   usr_def_zgr   ! called by domzgr.F90
33
34  !! * Substitutions
35#  include "vectopt_loop_substitute.h90"
36   !!----------------------------------------------------------------------
37   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
38   !! $Id$
39   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
40   !!----------------------------------------------------------------------
41CONTAINS             
42
43   SUBROUTINE usr_def_zgr( ld_zco  , ld_zps  , ld_sco  , ld_isfcav,    &   ! type of vertical coordinate
44      &                    pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d  ,    &   ! 1D reference vertical coordinate
45      &                    pdept , pdepw ,                             &   ! 3D t & w-points depth
46      &                    pe3t  , pe3u  , pe3v , pe3f ,               &   ! vertical scale factors
47      &                    pe3w  , pe3uw , pe3vw,                      &   !     -      -      -
48      &                    k_top  , k_bot    )                             ! top & bottom ocean level
49      !!---------------------------------------------------------------------
50      !!              ***  ROUTINE usr_def_zgr  ***
51      !!
52      !! ** Purpose :   User defined the vertical coordinates
53      !!
54      !!----------------------------------------------------------------------
55      LOGICAL                   , INTENT(out) ::   ld_zco, ld_zps, ld_sco      ! vertical coordinate flags
56      LOGICAL                   , INTENT(out) ::   ld_isfcav                   ! under iceshelf cavity flag
57      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth     [m]
58      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d           ! 1D grid-point depth     [m]
59      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pdept, pdepw                ! grid-point depth        [m]
60      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors  [m]
61      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3w , pe3uw, pe3vw         ! i-scale factors
62      INTEGER , DIMENSION(:,:)  , INTENT(out) ::   k_top, k_bot                ! first & last ocean level
63      !
64      INTEGER  ::   ji, jj, jk        ! dummy indices
65      INTEGER  ::   ik                ! local integers
66      REAL(wp) ::   zfact, z1_jpkm1   ! local scalar
67      REAL(wp) ::   ze3min            ! local scalar
68      REAL(wp) ::   zi, zj, zbathy    ! local scalar
69      REAL(wp) ::   ztmpu, ztmpv, ztmpf, ztmpu1, ztmpv1, ztmpf1
70      REAL(wp), DIMENSION(jpi,jpj) ::   zht, zhu, zhv, z2d   ! 2D workspace
71      !!----------------------------------------------------------------------
72      !
73      IF(lwp) WRITE(numout,*)
74      IF(lwp) WRITE(numout,*) 'usr_def_zgr : WAD_TEST_CASES configuration (s-coordinate closed box ocean without cavities)'
75      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
76      !
77      !
78      ! type of vertical coordinate
79      ! ---------------------------
80      ! set in usrdef_nam.F90 by reading the namusr_def namelist only ln_zco
81      ln_zco    = .FALSE.      ! z-partial-step coordinate
82      ln_zps    = .FALSE.      ! z-partial-step coordinate
83      ln_sco    = .TRUE.       ! s-coordinate
84      ld_isfcav = .FALSE.      ! ISF Ice Shelves Flag
85      !
86      !
87      ! Build the vertical coordinate system
88      ! ------------------------------------
89      !
90      !                       !==  UNmasked meter bathymetry  ==!
91      !
92!
93      zbathy=10.0
94      IF( cn_cfg == 'wad' ) THEN
95         SELECT CASE ( nn_cfg )
96            !                                        ! ====================
97            CASE ( 1 )                               ! WAD 1 configuration
98               !                                     ! ====================
99               !
100               IF(lwp) WRITE(numout,*)
101               IF(lwp) WRITE(numout,*) 'usr_def_zgr (WAD) : Closed box with EW linear bottom slope'
102               IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
103               !
104               zht = 1.5_wp
105               DO ji = 1, jpi
106                 zi = MIN((glamt(ji,1) - 10.0)/40.0, 1.0 )
107                 zht(ji,:) = MAX(zbathy*zi, 1.5) 
108                 IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zht(ji,1)
109               END DO
110               zht(mi0(1):mi1(1),:) = -4._wp
111               zht(mi0(jpiglo):mi1(jpiglo),:) = -4._wp
112               zht(:,mj0(1):mj1(1)) = -4._wp
113               zht(:,mj0(jpjglo):mj1(jpjglo)) = -4._wp
114               !                                     ! ====================
115            CASE ( 2, 3 )                            ! WAD 2 or 3  configuration
116               !                                     ! ====================
117               !
118               IF(lwp) WRITE(numout,*)
119               IF(lwp) WRITE(numout,*) 'usr_def_zgr (WAD) : Parobolic EW channel'
120               IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
121               !
122               DO ji = 1, jpi
123                 zi = MAX(1.0-((glamt(ji,1)-25.0)**2)/484.0, 0.0 )
124                 zi = MAX(1.0-((glamt(ji,1)-25.0)**2)/484.0, -2.0 )
125                 zht(ji,:) = MAX(zbathy*zi, -20.0)
126                 IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zht(ji,1)
127               END DO
128               zht(mi0(1):mi1(2),:) = -4._wp
129               zht(mi0(jpiglo-1):mi1(jpiglo),:) = -4._wp
130               zht(:,mj0(1):mj1(1)) = -4._wp
131               zht(:,mj0(jpjglo):mj1(jpjglo)) = -4._wp
132               !                                     ! ====================
133            CASE ( 4 )                               ! WAD 4 configuration
134               !                                     ! ====================
135               !
136               IF(lwp) WRITE(numout,*)
137               IF(lwp) WRITE(numout,*) 'usr_def_zgr (WAD) : Parobolic bowl'
138               IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
139               !
140               DO ji = 1, jpi
141                 zi = MAX(1.0-((glamt(ji,1)-25.0)**2)/484.0, -2.0 )
142               DO jj = 1, jpj
143                 zj = MAX(1.0-((gphit(1,jj)-17.0)**2)/196.0, -2.0 )
144                 zht(ji,jj) = MAX(zbathy*zi*zj, -2.0)
145               END DO
146                 IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zht(ji,1)
147               END DO
148               zht(mi0(1):mi1(1),:) = -4._wp
149               zht(mi0(jpiglo-1):mi1(jpiglo),:) = -4._wp
150               zht(:,mj0(1):mj1(1)) = -4._wp
151               zht(:,mj0(jpjglo):mj1(jpjglo)) = -4._wp
152               !                                    ! ===========================
153            CASE ( 5 )                              ! WAD 5 configuration
154               !                                    ! ====================
155               !
156               IF(lwp) WRITE(numout,*)
157               IF(lwp) WRITE(numout,*) 'usr_def_zgr (WAD) : Double slope with shelf'
158               IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
159               !
160               DO ji = 1, jpi
161                 zi = MIN(glamt(ji,1)/45.0, 1.0 )
162                 zht(ji,:) = MAX(zbathy*zi, 0.5)
163                 IF( glamt(ji,1) >= 46.0 ) THEN
164                   zht(ji,:) = 10.0
165                 ELSE IF( glamt(ji,1) >= 20.0 .AND. glamt(ji,1) < 46.0 ) THEN
166                   zi = 7.5/25.
167                   zht(ji,:) = MAX(10. - zi*(47.-glamt(ji,1)),2.5)
168                 ELSE IF( glamt(ji,1) >= 15.0 .AND. glamt(ji,1) < 20.0 ) THEN
169                   zht(ji,:) = 2.5
170                 ELSE IF( glamt(ji,1) >= 4.0 .AND. glamt(ji,1) < 15.0 ) THEN
171                   zi = 2.0/11.0
172                   zht(ji,:) = MAX(2.5 - zi*(16.0-glamt(ji,1)), 0.5)
173                 ELSE IF( glamt(ji,1) >= 1.0 .AND. glamt(ji,1) < 4.0 ) THEN
174                   zht(ji,:) = 0.5
175                 ENDIF
176                 IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zht(ji,1)
177               END DO
178               !                                    ! ===========================
179               zht(mi0(1):mi1(1),:) = -4._wp
180               zht(mi0(jpiglo):mi1(jpiglo),:) = -4._wp
181               zht(:,mj0(1):mj1(1)) = -4._wp
182               zht(:,mj0(jpjglo):mj1(jpjglo)) = -4._wp
183               !                                    ! ===========================
184            CASE ( 6 )                              ! WAD 6 configuration
185               !                                    ! ====================
186               !
187               IF(lwp) WRITE(numout,*)
188               IF(lwp) WRITE(numout,*) 'usr_def_zgr (WAD) : Parabolic channel with gaussian ridge'
189               IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
190               !
191               DO ji = 1, jpi
192                 zi = MAX(1.0-((glamt(ji,1)-25.0)**2)/484.0, -2.0 )
193                 zj = 0.95*MAX(EXP(-1.0*((glamt(ji,1)-25.0)**2)/32.0) , 0.0 )
194                 zht(ji,:) = MAX(zbathy*(zi-zj), -2.0)
195                 IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zht(ji,1)
196               END DO
197               zht(mi0(1):mi1(1),:) = -4._wp
198               zht(mi0(jpiglo):mi1(jpiglo),:) = -4._wp
199               zht(:,mj0(1):mj1(1)) = -4._wp
200               zht(:,mj0(jpjglo):mj1(jpjglo)) = -4._wp
201               !                                    ! ===========================
202            CASE DEFAULT
203               !                                    ! ===========================
204               WRITE(ctmp1,*) 'WAD test with a ', nn_cfg,' option is not coded'
205               !
206               CALL ctl_stop( ctmp1 )
207               !
208         END SELECT
209      END IF
210               !
211      !
212      ! at u-point: averaging zht
213      DO ji = 1, jpim1
214         zhu(ji,:) = 0.5_wp * ( zht(ji,:) + zht(ji+1,:) )
215      END DO
216      IF(lwp) write(numout,*) 'E3a zu ', MINVAL(zht(:,:)), MINVAL(zhu(:,:)), MINVAL(zhv(:,:))
217      CALL lbc_lnk( zhu, 'U', 1. )     ! boundary condition: this mask the surrounding grid-points
218      !                                ! ==>>>  set by hand non-zero value on first/last columns & rows
219      IF(lwp) write(numout,*) 'E3b zu ', MINVAL(zht(:,:)), MINVAL(zhu(:,:)), MINVAL(zhv(:,:))
220      DO ji = mi0(1), mi1(1)              ! first row of global domain only
221         zhu(ji,:) = zht(1,:)
222      END DO
223       DO ji = mi0(jpiglo), mi1(jpiglo)   ! last  row of global domain only
224         zhu(ji,:) = zht(jpi,:)
225      END DO
226      IF(lwp) write(numout,*) 'E3c zu ', MINVAL(zht(:,:)), MINVAL(zhu(:,:)), MINVAL(zhv(:,:))
227      ! at v-point: averaging zht
228      zhv = 0._wp
229      DO jj = 1, jpjm1
230         zhv(:,jj) = 0.5_wp * ( zht(:,jj) + zht(:,jj+1) )
231      END DO
232      CALL lbc_lnk( zhv, 'V', 1. )     ! boundary condition: this mask the surrounding grid-points
233      IF(lwp) write(numout,*) 'E3d zu ', MINVAL(zht(:,:)), MINVAL(zhu(:,:)), MINVAL(zhv(:,:))
234      !avoid the zero depth on T- (U-,V-,F-) points
235        DO jj = 1, jpj
236          DO ji = 1, jpi
237            IF(ABS(zht(ji,jj)) < rn_wdmin1) &
238              & zht(ji,jj) = SIGN(1._wp, zht(ji,jj)) * rn_wdmin1
239            IF(ABS(zhu(ji,jj)) < rn_wdmin1) &
240              & zhu(ji,jj) = SIGN(1._wp, zhu(ji,jj)) * rn_wdmin1
241            IF(ABS(zhv(ji,jj)) < rn_wdmin1) &
242              & zhv(ji,jj) = SIGN(1._wp, zhv(ji,jj)) * rn_wdmin1
243          END DO
244        END DO
245      IF(lwp) write(numout,*) 'E3e zu ', MINVAL(zht(:,:)), MINVAL(zhu(:,:)), MINVAL(zhv(:,:))
246      !     
247      CALL zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! Reference z-coordinate system
248      !
249      !
250      !                       !==  top masked level bathymetry  ==!  (all coordinates)
251      !
252      ! no ocean cavities : top ocean level is ONE, except over land
253      ! the ocean basin surrounnded by land (1 grid-point) set through lbc_lnk call as jperio=0
254      z2d(:,:) = 1._wp                    ! surface ocean is the 1st level
255      z2d(mi0(1):mi1(1),:) = 0._wp
256      z2d(mi0(jpiglo):mi1(jpiglo),:) = 0._wp
257      z2d(:,mj0(1):mj1(1)) = 0._wp
258      z2d(:,mj0(jpjglo):mj1(jpjglo)) = 0._wp
259
260
261
262
263
264      CALL lbc_lnk( z2d, 'T', 1. )        ! closed basin since jperio = 0 (see userdef_nam.F90)
265      k_top(:,:) = NINT( z2d(:,:) )
266      !
267      !                             
268      !
269      IF ( ln_sco ) THEN      !==  s-coordinate  ==!   (terrain-following coordinate)
270         !
271         ht_wd = zht
272         k_bot(:,:) = jpkm1 * k_top(:,:)  !* bottom ocean = jpk-1 (here use k_top as a land mask)
273         DO jj = 1, jpj
274            DO ji = 1, jpi
275              IF( zht(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN
276                !WRITE(numout,*) 'KBOT ',ji,jj,zht(ji,jj),rn_wdld,rn_wdmin2
277                k_bot(ji,jj) = 0
278                zht(ji,jj) = rn_wdmin1 * REAL( jpkm1 , wp) * 0.5_wp
279              ELSEIF(zht(ji,jj) <= rn_wdmin1) THEN
280                k_bot(ji,jj) = 2
281                zht(ji,jj) = rn_wdmin1 * REAL( jpkm1 , wp) * 0.5_wp
282                !WRITE(numout,*) 'KBOT1 ',ji,jj,zht(ji,jj),rn_wdmin1
283              ENDIF
284           END DO
285         END DO
286      DO ji = 1, jpim1
287         zhu(ji,:) = 0.5_wp * ( zht(ji,:) + zht(ji+1,:) )
288      END DO
289      IF(lwp) write(numout,*) 'E3f zu ', MINVAL(zht(:,:)), MINVAL(zhu(:,:)), MINVAL(zhv(:,:))
290      !                                ! ==>>>  set by hand non-zero value on
291      !                                first/last columns & rows
292      IF(lwp) write(numout,*) 'E3f Zu ', MINVAL(zht(:,:)), MINVAL(zhu(:,:)), MINVAL(zhv(:,:))
293      CALL lbc_lnk( zhu, 'U', 1. )     ! boundary condition: this mask the surrounding grid-points
294      IF(lwp) write(numout,*) 'E3F zu ', MINVAL(zht(:,:)), MINVAL(zhu(:,:)), MINVAL(zhv(:,:))
295      DO ji = mi0(1), mi1(1)              ! first row of global domain only
296         zhu(ji,:) = zht(ji,:)
297      END DO
298      DO ji = mi0(jpiglo), mi1(jpiglo)   ! last  row of global domain only
299         zhu(ji,:) = zht(ji,:)
300      END DO
301      ! at v-point: averaging zht
302      IF(lwp) write(numout,*) 'E3g zu ', MINVAL(zht(:,:)), MINVAL(zhu(:,:)), MINVAL(zhv(:,:))
303      zhv = 0._wp
304      DO jj = 1, jpjm1
305         zhv(:,jj) = 0.5_wp * ( zht(:,jj) + zht(:,jj+1) )
306      END DO
307      IF(lwp) write(numout,*) 'E3h zu ', MINVAL(zht(:,:)), MINVAL(zhu(:,:)), MINVAL(zhv(:,:))
308      CALL lbc_lnk( zhv, 'V', 1. )
309      IF(lwp) write(numout,*) 'E3H zu ', MINVAL(zht(:,:)), MINVAL(zhu(:,:)), MINVAL(zhv(:,:))
310      DO jj = mj0(1), mj1(1)   ! first  row of global domain only
311         zhv(:,jj) = zht(:,jj)
312      END DO
313      DO jj = mj0(jpjglo), mj1(jpjglo)   ! last  row of global domain only
314         zhv(:,jj) = zht(:,jj)
315      END DO
316         !
317         !                                !* terrain-following coordinate with e3.(k)=cst)
318         !                                !  OVERFLOW case : identical with j-index (T=V, U=F)
319         z1_jpkm1 = 1._wp / REAL( jpkm1 , wp)
320         DO jk = 1, jpk
321                  pdept(:,:,jk) = zht(:,:) * z1_jpkm1 * ( REAL( jk   , wp ) - 0.5_wp )
322                  pdepw(:,:,jk) = zht(:,:) * z1_jpkm1 * ( REAL( jk-1 , wp )          )
323                  pe3t (:,:,jk) = zht(:,:) * z1_jpkm1
324                  pe3u (:,:,jk) = zhu(:,:) * z1_jpkm1
325                  pe3v (:,:,jk) = zhv(:,:) * z1_jpkm1
326                  pe3f (:,:,jk) = zhu(:,:) * z1_jpkm1
327                  pe3w (:,:,jk) = zht(:,:) * z1_jpkm1
328                  pe3uw(:,:,jk) = zhu(:,:) * z1_jpkm1
329                  pe3vw(:,:,jk) = zhv(:,:) * z1_jpkm1
330         END DO     
331         DO jj = 1, jpjm1
332            DO ji = 1, jpim1
333               ztmpu  = zht(ji,jj)+zht(ji+1,jj)
334               ztmpv  = zht(ji,jj)+zht(ji,jj+1)
335               ztmpf  = zht(ji,jj)+zht(ji+1,jj)+zht(ji,jj+1)+zht(ji+1,jj+1)
336               ztmpu1 = zht(ji,jj)*zht(ji+1,jj)
337               ztmpv1 = zht(ji,jj)*zht(ji,jj+1)
338               ztmpf1 = MIN(zht(ji,jj), zht(ji+1,jj), zht(ji,jj+1), zht(ji+1,jj+1)) * &
339                      & MAX(zht(ji,jj), zht(ji+1,jj), zht(ji,jj+1), zht(ji+1,jj+1))
340               IF(  (ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1) ) THEN
341                 pe3u (ji,jj,:) = 0.5_wp* ( pe3t(ji,jj,:) + pe3t(ji+1,jj,:) )
342                 pe3uw(ji,jj,:) = 0.5_wp* ( pe3w(ji,jj,:) + pe3w(ji+1,jj,:) )
343               ENDIF
344               IF(  (ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1) ) THEN
345                 pe3v (ji,jj,:) = 0.5_wp* ( pe3t(ji,jj,:) + pe3t(ji,jj+1,:) )
346                 pe3vw(ji,jj,:) = 0.5_wp* ( pe3w(ji,jj,:) + pe3w(ji,jj+1,:) )
347               ENDIF
348               IF(  (ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1) ) THEN
349                 pe3f (ji,jj,:) = 0.25_wp* ( pe3t(ji,jj,:)   + pe3t(ji+1,jj,:) &
350                                &          + pe3t(ji,jj+1,:) + pe3t(ji+1,jj+1,:) )
351               ENDIF
352                 pe3f (ji,jj,:) = 0.25_wp* ( pe3t(ji,jj,:)   + pe3t(ji+1,jj,:) &
353                                &          + pe3t(ji,jj+1,:) + pe3t(ji+1,jj+1,:) )
354           END DO
355         END DO
356         CALL lbc_lnk( pe3u , 'U', 1. )
357         CALL lbc_lnk( pe3uw, 'U', 1. )
358         CALL lbc_lnk( pe3v , 'V', 1. )
359         CALL lbc_lnk( pe3vw, 'V', 1. )
360         CALL lbc_lnk( pe3f , 'F', 1. )
361      WHERE( pe3t (:,:,:) == 0._wp )   pe3t (:,:,:) = 1._wp
362      WHERE( pe3u (:,:,:) == 0._wp )   pe3u (:,:,:) = 1._wp
363      WHERE( pe3v (:,:,:) == 0._wp )   pe3v (:,:,:) = 1._wp
364      WHERE( pe3f (:,:,:) == 0._wp )   pe3f (:,:,:) = 1._wp
365      WHERE( pe3w (:,:,:) == 0._wp )   pe3w (:,:,:) = 1._wp
366      WHERE( pe3uw(:,:,:) == 0._wp )   pe3uw(:,:,:) = 1._wp
367      WHERE( pe3vw(:,:,:) == 0._wp )   pe3vw(:,:,:) = 1._wp
368      IF(lwp) write(numout,*) 'E3 TU ',SUM(0._wp/pe3t(:,:,:)), MINVAL(pe3t(:,:,:)), rn_wdmin1 ; CALL FLUSH(numout)
369      IF(lwp) write(numout,*) 'E3 tu ',SUM(0._wp/pe3u(:,:,:)), MINVAL(pe3u(:,:,:))
370      IF(lwp) write(numout,*) 'E3 tv ',SUM(0._wp/pe3v(:,:,:)), MINVAL(pe3v(:,:,:))
371      IF(lwp) write(numout,*) 'E3 tf ',SUM(0._wp/pe3f(:,:,:)), MINVAL(pe3f(:,:,:))
372      IF(lwp) write(numout,*) 'E3 zu ', MINVAL(zht(:,:)), MINVAL(zhu(:,:)), MINVAL(zhv(:,:))
373      ENDIF
374      !
375      !
376      !
377   END SUBROUTINE usr_def_zgr
378
379
380   SUBROUTINE zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! 1D reference vertical coordinate
381      !!----------------------------------------------------------------------
382      !!                   ***  ROUTINE zgr_z1d  ***
383      !!
384      !! ** Purpose :   set the depth of model levels and the resulting
385      !!      vertical scale factors.
386      !!
387      !! ** Method  :   1D z-coordinate system (use in all type of coordinate)
388      !!       The depth of model levels is set from dep(k), an analytical function:
389      !!                   w-level: depw_1d  = dep(k)
390      !!                   t-level: dept_1d  = dep(k+0.5)
391      !!       The scale factors are the discrete derivative of the depth:
392      !!                   e3w_1d(jk) = dk[ dept_1d ]
393      !!                   e3t_1d(jk) = dk[ depw_1d ]
394      !!
395      !!            ===    Here constant vertical resolution   ===
396      !!
397      !! ** Action  : - pdept_1d, pdepw_1d : depth of T- and W-point (m)
398      !!              - pe3t_1d , pe3w_1d  : scale factors at T- and W-levels (m)
399      !!----------------------------------------------------------------------
400      REAL(wp), DIMENSION(:), INTENT(out) ::   pdept_1d, pdepw_1d   ! 1D grid-point depth        [m]
401      REAL(wp), DIMENSION(:), INTENT(out) ::   pe3t_1d , pe3w_1d    ! 1D vertical scale factors  [m]
402      !
403      INTEGER  ::   jk       ! dummy loop indices
404      REAL(wp) ::   zt, zw   ! local scalar
405      !!----------------------------------------------------------------------
406      !
407      IF(lwp) THEN                         ! Parameter print
408         WRITE(numout,*)
409         WRITE(numout,*) '    zgr_z1d : Reference vertical z-coordinates: uniform dz = ', rn_dz
410         WRITE(numout,*) '    ~~~~~~~'
411      ENDIF
412      !
413      ! Reference z-coordinate (depth - scale factor at T- and W-points)   ! Madec & Imbard 1996 function
414      ! ----------------------
415      DO jk = 1, jpk
416         zw = REAL( jk , wp )
417         zt = REAL( jk , wp ) + 0.5_wp
418         pdepw_1d(jk) =    rn_dz *   REAL( jk-1 , wp )
419         pdept_1d(jk) =    rn_dz * ( REAL( jk-1 , wp ) + 0.5_wp )
420         pe3w_1d (jk) =    rn_dz
421         pe3t_1d (jk) =    rn_dz
422      END DO
423      !
424      IF(lwp) THEN                        ! control print
425         WRITE(numout,*)
426         WRITE(numout,*) '              Reference 1D z-coordinate depth and scale factors:'
427         WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
428         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
429      ENDIF
430      !
431   END SUBROUTINE zgr_z1d
432   
433   !!======================================================================
434END MODULE usrdef_zgr
Note: See TracBrowser for help on using the repository browser.