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

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

Branch 2016/dev_merge_2016. Wetting and drying updates. Includes a slightly modified MY_SRC copy of domain.F90 to enable a north-south cyclic test domain (jperio=7). Other changes include fixes in usrdef_zgr.F90 to correct metrics on the northern boundary and a new test case (#8) which uses the north-south cyclic option.
File size: 18.5 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, zwet
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, -2.0) 
108               END DO
109               zht(mi0(1):mi1(1),:) = -4._wp
110               zht(mi0(jpiglo):mi1(jpiglo),:) = -4._wp
111               zht(:,mj0(1):mj1(1)) = -4._wp
112               zht(:,mj0(jpjglo):mj1(jpjglo)) = -4._wp
113               !                                     ! ====================
114            CASE ( 2, 3, 8 )                         ! WAD 2 or 3  configuration
115               !                                     ! ====================
116               !
117               IF(lwp) WRITE(numout,*)
118               IF(lwp) WRITE(numout,*) 'usr_def_zgr (WAD) : Parobolic EW channel'
119               IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
120               !
121               DO ji = 1, jpi
122                 zi = MAX(1.0-((glamt(ji,1)-25.0)**2)/484.0, -0.3 )
123                 zht(ji,:) = MAX(zbathy*zi, -2.0)
124               END DO
125               zht(mi0(1):mi1(1),:) = -4._wp
126               zht(mi0(jpiglo):mi1(jpiglo),:) = -4._wp
127               IF( nn_cfg /= 8 ) THEN
128                  zht(:,mj0(1):mj1(1)) = -4._wp
129                  zht(:,mj0(jpjglo):mj1(jpjglo)) = -4._wp
130               ENDIF
131               !                                     ! ====================
132            CASE ( 4 )                               ! WAD 4 configuration
133               !                                     ! ====================
134               !
135               IF(lwp) WRITE(numout,*)
136               IF(lwp) WRITE(numout,*) 'usr_def_zgr (WAD) : Parobolic bowl'
137               IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
138               !
139               DO ji = 1, jpi
140                 zi = MAX(1.0-((glamt(ji,1)-25.0)**2)/484.0, -2.0 )
141               DO jj = 1, jpj
142                 zj = MAX(1.0-((gphit(1,jj)-17.0)**2)/196.0, -2.0 )
143                 zht(ji,jj) = MAX(zbathy*zi*zj, -2.0)
144               END DO
145               END DO
146               zht(mi0(1):mi1(1),:) = -4._wp
147               zht(mi0(jpiglo):mi1(jpiglo),:) = -4._wp
148               zht(:,mj0(1):mj1(1)) = -4._wp
149               zht(:,mj0(jpjglo):mj1(jpjglo)) = -4._wp
150               !                                    ! ===========================
151            CASE ( 5 )                              ! WAD 5 configuration
152               !                                    ! ====================
153               !
154               IF(lwp) WRITE(numout,*)
155               IF(lwp) WRITE(numout,*) 'usr_def_zgr (WAD) : Double slope with shelf'
156               IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
157               !
158               DO ji = 1, jpi
159                 zi = MIN(glamt(ji,1)/45.0, 1.0 )
160                 zht(ji,:) = MAX(zbathy*zi, -2.0)
161                 IF( glamt(ji,1) >= 46.0 ) THEN
162                   zht(ji,:) = 10.0
163                 ELSE IF( glamt(ji,1) >= 20.0 .AND. glamt(ji,1) < 46.0 ) THEN
164                   zi = 7.5/25.
165                   zht(ji,:) = MAX(10. - zi*(47.-glamt(ji,1)),2.5)
166                 ELSE IF( glamt(ji,1) >= 15.0 .AND. glamt(ji,1) < 20.0 ) THEN
167                   zht(ji,:) = 2.5
168                 ELSE IF( glamt(ji,1) >= 4.0 .AND. glamt(ji,1) < 15.0 ) THEN
169                   zi = 4.5/11.0
170                   zht(ji,:) = MAX(2.5 - zi*(16.0-glamt(ji,1)), -2.0)
171                 ELSE IF( glamt(ji,1) >= 0.0 .AND. glamt(ji,1) < 4.0 ) THEN
172                   zht(ji,:) = -2.0
173                 ENDIF
174               END DO
175               !                                    ! ===========================
176               zht(mi0(1):mi1(1),:) = -4._wp
177               zht(mi0(jpiglo):mi1(jpiglo),:) = -4._wp
178               zht(:,mj0(1):mj1(1)) = -4._wp
179               zht(:,mj0(jpjglo):mj1(jpjglo)) = -4._wp
180               !                                    ! ===========================
181            CASE ( 6 )                              ! WAD 6 configuration
182               !                                    ! ====================
183               !
184               IF(lwp) WRITE(numout,*)
185               IF(lwp) WRITE(numout,*) 'usr_def_zgr (WAD) : Parabolic channel with gaussian ridge'
186               IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
187               !
188               DO ji = 1, jpi
189                 zi = MAX(1.0-((glamt(ji,1)-25.0)**2)/484.0, -2.0 )
190                 zj = 1.075*MAX(EXP(-1.0*((glamt(ji,1)-25.0)**2)/32.0) , 0.0 )
191                 zht(ji,:) = MAX(zbathy*(zi-zj), -2.0)
192               END DO
193               zht(mi0(1):mi1(1),:) = -4._wp
194               zht(mi0(jpiglo):mi1(jpiglo),:) = -4._wp
195               zht(:,mj0(1):mj1(1)) = -4._wp
196               zht(:,mj0(jpjglo):mj1(jpjglo)) = -4._wp
197               !                                    ! ===========================
198            CASE ( 7 )                              ! WAD 7 configuration
199               !                                    ! ====================
200               !
201               IF(lwp) WRITE(numout,*)
202               IF(lwp) WRITE(numout,*) 'usr_def_zgr (WAD) : Double slope with open boundary'
203               IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
204               !
205               DO ji = 1, jpi
206                 zi = MIN(glamt(ji,1)/45.0, 1.0 )
207                 zht(ji,:) = MAX(zbathy*zi, -2.0)
208                 IF( glamt(ji,1) >= 46.0 ) THEN
209                   zht(ji,:) = 10.0
210                 ELSE IF( glamt(ji,1) >= 20.0 .AND. glamt(ji,1) < 46.0 ) THEN
211                   zi = 7.5/25.
212                   zht(ji,:) = MAX(10. - zi*(47.-glamt(ji,1)),2.5)
213                 ELSE IF( glamt(ji,1) >= 15.0 .AND. glamt(ji,1) < 20.0 ) THEN
214                   zht(ji,:) = 2.5
215                 ELSE IF( glamt(ji,1) >= 4.0 .AND. glamt(ji,1) < 15.0 ) THEN
216                   zi = 4.5/11.0
217                   zht(ji,:) = MAX(2.5 - zi*(16.0-glamt(ji,1)), -2.0)
218                 ELSE IF( glamt(ji,1) >= 0.0 .AND. glamt(ji,1) < 4.0 ) THEN
219                   zht(ji,:) = -2.0
220                 ENDIF
221               END DO
222               !                                    ! ===========================
223               zht(mi0(1):mi1(1),:) = -4._wp
224               zht(:,mj0(1):mj1(1)) = -4._wp
225               zht(:,mj0(jpjglo):mj1(jpjglo)) = -4._wp
226            CASE DEFAULT
227               !                                    ! ===========================
228               WRITE(ctmp1,*) 'WAD test with a ', nn_cfg,' option is not coded'
229               !
230               CALL ctl_stop( ctmp1 )
231               !
232         END SELECT
233      END IF
234               !
235      !
236      ! at u-point: averaging zht
237      DO ji = 1, jpim1
238         zhu(ji,:) = 0.5_wp * ( zht(ji,:) + zht(ji+1,:) )
239      END DO
240      CALL lbc_lnk( zhu, 'U', 1. )     ! boundary condition: this mask the surrounding grid-points
241      !                                ! ==>>>  set by hand non-zero value on first/last columns & rows
242      DO ji = mi0(1), mi1(1)              ! first row of global domain only
243         zhu(ji,:) = zht(1,:)
244      END DO
245       DO ji = mi0(jpiglo), mi1(jpiglo)   ! last  row of global domain only
246         zhu(ji,:) = zht(jpi,:)
247      END DO
248      ! at v-point: averaging zht
249      zhv = 0._wp
250      DO jj = 1, jpjm1
251         zhv(:,jj) = 0.5_wp * ( zht(:,jj) + zht(:,jj+1) )
252      END DO
253      CALL lbc_lnk( zhv, 'V', 1. )     ! boundary condition: this mask the surrounding grid-points
254      DO jj = mj0(1), mj1(1)   ! first  row of global domain only
255         zhv(:,jj) = zht(:,jj)
256      END DO
257      DO jj = mj0(jpjglo), mj1(jpjglo)   ! last  row of global domain only
258         zhv(:,jj) = zht(:,jj)
259      END DO
260      !     
261      CALL zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! Reference z-coordinate system
262      !
263      !
264      !                       !==  top masked level bathymetry  ==!  (all coordinates)
265      !
266      ! no ocean cavities : top ocean level is ONE, except over land
267      ! the ocean basin surrounnded by land (1 grid-point) set through lbc_lnk call as jperio=0
268      z2d(:,:) = 1._wp                    ! surface ocean is the 1st level
269      z2d(mi0(1):mi1(1),:) = 0._wp
270      z2d(mi0(jpiglo):mi1(jpiglo),:) = 0._wp
271      z2d(:,mj0(1):mj1(1)) = 0._wp
272      z2d(:,mj0(jpjglo):mj1(jpjglo)) = 0._wp
273
274
275
276
277
278      CALL lbc_lnk( z2d, 'T', 1. )        ! closed basin since jperio = 0 (see userdef_nam.F90)
279      k_top(:,:) = NINT( z2d(:,:) )
280      !
281      !                             
282      !
283      IF ( ln_sco ) THEN      !==  s-coordinate  ==!   (terrain-following coordinate)
284         !
285         ht_wd = zht
286         k_bot(:,:) = jpkm1 * k_top(:,:)  !* bottom ocean = jpk-1 (here use k_top as a land mask)
287         DO jj = 1, jpj
288            DO ji = 1, jpi
289              IF( zht(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN
290                k_bot(ji,jj) = 0
291                k_top(ji,jj) = 0
292              ENDIF
293           END DO
294         END DO
295         !
296         !                                !* terrain-following coordinate with e3.(k)=cst)
297         !                                !  OVERFLOW case : identical with j-index (T=V, U=F)
298         DO jj = 1, jpjm1
299            DO ji = 1, jpim1
300              z1_jpkm1 = 1._wp / REAL( k_bot(ji,jj) - k_top(ji,jj) + 1 , wp)
301              DO jk = 1, jpk
302                  zwet = MAX( zht(ji,jj), rn_wdmin1 )
303                  pdept(ji,jj,jk) = zwet * z1_jpkm1 * ( REAL( jk   , wp ) - 0.5_wp )
304                  pdepw(ji,jj,jk) = zwet * z1_jpkm1 * ( REAL( jk-1 , wp )          )
305                  pe3t (ji,jj,jk) = zwet * z1_jpkm1
306                  pe3w (ji,jj,jk) = zwet * z1_jpkm1
307                  zwet = MAX( zhu(ji,jj), rn_wdmin1 )
308                  pe3u (ji,jj,jk) = zwet * z1_jpkm1
309                  pe3uw(ji,jj,jk) = zwet * z1_jpkm1
310                  pe3f (ji,jj,jk) = zwet * z1_jpkm1
311                  zwet = MAX( zhv(ji,jj), rn_wdmin1 )
312                  pe3v (ji,jj,jk) = zwet * z1_jpkm1
313                  pe3vw(ji,jj,jk) = zwet * z1_jpkm1
314              END DO     
315           END DO     
316         END DO     
317         CALL lbc_lnk( pdept, 'T', 1. )
318         CALL lbc_lnk( pdepw, 'T', 1. )
319         CALL lbc_lnk( pe3t , 'T', 1. )
320         CALL lbc_lnk( pe3w , 'T', 1. )
321         CALL lbc_lnk( pe3u , 'U', 1. )
322         CALL lbc_lnk( pe3uw, 'U', 1. )
323         CALL lbc_lnk( pe3f , 'F', 1. )
324         CALL lbc_lnk( pe3v , 'V', 1. )
325         CALL lbc_lnk( pe3vw, 'V', 1. )
326         WHERE( pe3t (:,:,:) == 0._wp )   pe3t (:,:,:) = 1._wp
327         WHERE( pe3u (:,:,:) == 0._wp )   pe3u (:,:,:) = 1._wp
328         WHERE( pe3v (:,:,:) == 0._wp )   pe3v (:,:,:) = 1._wp
329         WHERE( pe3f (:,:,:) == 0._wp )   pe3f (:,:,:) = 1._wp
330         WHERE( pe3w (:,:,:) == 0._wp )   pe3w (:,:,:) = 1._wp
331         WHERE( pe3uw(:,:,:) == 0._wp )   pe3uw(:,:,:) = 1._wp
332         WHERE( pe3vw(:,:,:) == 0._wp )   pe3vw(:,:,:) = 1._wp
333      ENDIF
334      !
335      !
336      !
337   END SUBROUTINE usr_def_zgr
338
339
340   SUBROUTINE zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! 1D reference vertical coordinate
341      !!----------------------------------------------------------------------
342      !!                   ***  ROUTINE zgr_z1d  ***
343      !!
344      !! ** Purpose :   set the depth of model levels and the resulting
345      !!      vertical scale factors.
346      !!
347      !! ** Method  :   1D z-coordinate system (use in all type of coordinate)
348      !!       The depth of model levels is set from dep(k), an analytical function:
349      !!                   w-level: depw_1d  = dep(k)
350      !!                   t-level: dept_1d  = dep(k+0.5)
351      !!       The scale factors are the discrete derivative of the depth:
352      !!                   e3w_1d(jk) = dk[ dept_1d ]
353      !!                   e3t_1d(jk) = dk[ depw_1d ]
354      !!
355      !!            ===    Here constant vertical resolution   ===
356      !!
357      !! ** Action  : - pdept_1d, pdepw_1d : depth of T- and W-point (m)
358      !!              - pe3t_1d , pe3w_1d  : scale factors at T- and W-levels (m)
359      !!----------------------------------------------------------------------
360      REAL(wp), DIMENSION(:), INTENT(out) ::   pdept_1d, pdepw_1d   ! 1D grid-point depth        [m]
361      REAL(wp), DIMENSION(:), INTENT(out) ::   pe3t_1d , pe3w_1d    ! 1D vertical scale factors  [m]
362      !
363      INTEGER  ::   jk       ! dummy loop indices
364      REAL(wp) ::   zt, zw   ! local scalar
365      !!----------------------------------------------------------------------
366      !
367      IF(lwp) THEN                         ! Parameter print
368         WRITE(numout,*)
369         WRITE(numout,*) '    zgr_z1d : Reference vertical z-coordinates: uniform dz = ', rn_dz
370         WRITE(numout,*) '    ~~~~~~~'
371      ENDIF
372      !
373      ! Reference z-coordinate (depth - scale factor at T- and W-points)   ! Madec & Imbard 1996 function
374      ! ----------------------
375      DO jk = 1, jpk
376         zw = REAL( jk , wp )
377         zt = REAL( jk , wp ) + 0.5_wp
378         pdepw_1d(jk) =    rn_dz *   REAL( jk-1 , wp )
379         pdept_1d(jk) =    rn_dz * ( REAL( jk-1 , wp ) + 0.5_wp )
380         pe3w_1d (jk) =    rn_dz
381         pe3t_1d (jk) =    rn_dz
382      END DO
383      !
384      IF(lwp) THEN                        ! control print
385         WRITE(numout,*)
386         WRITE(numout,*) '              Reference 1D z-coordinate depth and scale factors:'
387         WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
388         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
389      ENDIF
390      !
391   END SUBROUTINE zgr_z1d
392   
393   !!======================================================================
394END MODULE usrdef_zgr
Note: See TracBrowser for help on using the repository browser.