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/2017/dev_merge_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_SRC – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_SRC/usrdef_zgr.F90 @ 9125

Last change on this file since 9125 was 9125, checked in by timgraham, 6 years ago

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