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

source: branches/UKMO/ROMS_WAD_7832/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_SRC/usrdef_zgr.F90 @ 8865

Last change on this file since 8865 was 8865, checked in by deazer, 6 years ago

Moving Changes from CS15mini config into NEMO main src
Updating TEST configs to run wit this version of the code
all sette tests pass at this revision other than AGRIF
Includes updates to dynnxt and tranxt required for 3D rives in WAD case to be conservative.

Next commit will update naming conventions and tidy the code.
File size: 18.7 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, ht_0   ! ocean space and time domain
20   USE usrdef_nam     ! User defined : namelist variables
21   USE wet_dry ,  ONLY: rn_wdmin1, rn_wdmin2, rn_wdld,  rn_ht_0,rn_ssh_ref
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! increase the depth of the bathymetry by rn_ssh_ref and rn_ht_0   
236      !zht(:,:) = zht(:,:) + rn_ssh_ref + rn_ht_0 
237      !zht(:,:) = zht(:,:) + rn_ssh_ref + rn_ht_0 
238
239      ! at u-point: averaging zht
240      DO ji = 1, jpim1
241         zhu(ji,:) = 0.5_wp * ( zht(ji,:) + zht(ji+1,:) )
242      END DO
243      CALL lbc_lnk( zhu, 'U', 1. )     ! boundary condition: this mask the surrounding grid-points
244      !                                ! ==>>>  set by hand non-zero value on first/last columns & rows
245      DO ji = mi0(1), mi1(1)              ! first row of global domain only
246         zhu(ji,:) = zht(1,:)
247      END DO
248       DO ji = mi0(jpiglo), mi1(jpiglo)   ! last  row of global domain only
249         zhu(ji,:) = zht(jpi,:)
250      END DO
251      ! at v-point: averaging zht
252      zhv = 0._wp
253      DO jj = 1, jpjm1
254         zhv(:,jj) = 0.5_wp * ( zht(:,jj) + zht(:,jj+1) )
255      END DO
256      CALL lbc_lnk( zhv, 'V', 1. )     ! boundary condition: this mask the surrounding grid-points
257      DO jj = mj0(1), mj1(1)   ! first  row of global domain only
258         zhv(:,jj) = zht(:,jj)
259      END DO
260      DO jj = mj0(jpjglo), mj1(jpjglo)   ! last  row of global domain only
261         zhv(:,jj) = zht(:,jj)
262      END DO
263      !     
264      CALL zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! Reference z-coordinate system
265      !
266      !
267      !                       !==  top masked level bathymetry  ==!  (all coordinates)
268      !
269      ! no ocean cavities : top ocean level is ONE, except over land
270      ! the ocean basin surrounnded by land (1 grid-point) set through lbc_lnk call as jperio=0
271      z2d(:,:) = 1._wp                    ! surface ocean is the 1st level
272      z2d(mi0(1):mi1(1),:) = 0._wp
273      z2d(mi0(jpiglo):mi1(jpiglo),:) = 0._wp
274      z2d(:,mj0(1):mj1(1)) = 0._wp
275      z2d(:,mj0(jpjglo):mj1(jpjglo)) = 0._wp
276
277
278
279
280
281      CALL lbc_lnk( z2d, 'T', 1. )        ! closed basin since jperio = 0 (see userdef_nam.F90)
282      k_top(:,:) = NINT( z2d(:,:) )
283      !
284      !                             
285      !
286      IF ( ln_sco ) THEN      !==  s-coordinate  ==!   (terrain-following coordinate)
287         !
288         ht_0 = zht
289         k_bot(:,:) = jpkm1 * k_top(:,:)  !* bottom ocean = jpk-1 (here use k_top as a land mask)
290         DO jj = 1, jpj
291            DO ji = 1, jpi
292              IF( zht(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN
293                k_bot(ji,jj) = 0
294                k_top(ji,jj) = 0
295              ENDIF
296           END DO
297         END DO
298         !
299         !                                !* terrain-following coordinate with e3.(k)=cst)
300         !                                !  OVERFLOW case : identical with j-index (T=V, U=F)
301         DO jj = 1, jpjm1
302            DO ji = 1, jpim1
303              z1_jpkm1 = 1._wp / REAL( k_bot(ji,jj) - k_top(ji,jj) + 1 , wp)
304              DO jk = 1, jpk
305                  zwet = MAX( zht(ji,jj), rn_wdmin1 )
306                  pdept(ji,jj,jk) = zwet * z1_jpkm1 * ( REAL( jk   , wp ) - 0.5_wp )
307                  pdepw(ji,jj,jk) = zwet * z1_jpkm1 * ( REAL( jk-1 , wp )          )
308                  pe3t (ji,jj,jk) = zwet * z1_jpkm1
309                  pe3w (ji,jj,jk) = zwet * z1_jpkm1
310                  zwet = MAX( zhu(ji,jj), rn_wdmin1 )
311                  pe3u (ji,jj,jk) = zwet * z1_jpkm1
312                  pe3uw(ji,jj,jk) = zwet * z1_jpkm1
313                  pe3f (ji,jj,jk) = zwet * z1_jpkm1
314                  zwet = MAX( zhv(ji,jj), rn_wdmin1 )
315                  pe3v (ji,jj,jk) = zwet * z1_jpkm1
316                  pe3vw(ji,jj,jk) = zwet * z1_jpkm1
317              END DO     
318           END DO     
319         END DO     
320         CALL lbc_lnk( pdept, 'T', 1. )
321         CALL lbc_lnk( pdepw, 'T', 1. )
322         CALL lbc_lnk( pe3t , 'T', 1. )
323         CALL lbc_lnk( pe3w , 'T', 1. )
324         CALL lbc_lnk( pe3u , 'U', 1. )
325         CALL lbc_lnk( pe3uw, 'U', 1. )
326         CALL lbc_lnk( pe3f , 'F', 1. )
327         CALL lbc_lnk( pe3v , 'V', 1. )
328         CALL lbc_lnk( pe3vw, 'V', 1. )
329         WHERE( pe3t (:,:,:) == 0._wp )   pe3t (:,:,:) = 1._wp
330         WHERE( pe3u (:,:,:) == 0._wp )   pe3u (:,:,:) = 1._wp
331         WHERE( pe3v (:,:,:) == 0._wp )   pe3v (:,:,:) = 1._wp
332         WHERE( pe3f (:,:,:) == 0._wp )   pe3f (:,:,:) = 1._wp
333         WHERE( pe3w (:,:,:) == 0._wp )   pe3w (:,:,:) = 1._wp
334         WHERE( pe3uw(:,:,:) == 0._wp )   pe3uw(:,:,:) = 1._wp
335         WHERE( pe3vw(:,:,:) == 0._wp )   pe3vw(:,:,:) = 1._wp
336      ENDIF
337      !
338      !
339      !
340   END SUBROUTINE usr_def_zgr
341
342
343   SUBROUTINE zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! 1D reference vertical coordinate
344      !!----------------------------------------------------------------------
345      !!                   ***  ROUTINE zgr_z1d  ***
346      !!
347      !! ** Purpose :   set the depth of model levels and the resulting
348      !!      vertical scale factors.
349      !!
350      !! ** Method  :   1D z-coordinate system (use in all type of coordinate)
351      !!       The depth of model levels is set from dep(k), an analytical function:
352      !!                   w-level: depw_1d  = dep(k)
353      !!                   t-level: dept_1d  = dep(k+0.5)
354      !!       The scale factors are the discrete derivative of the depth:
355      !!                   e3w_1d(jk) = dk[ dept_1d ]
356      !!                   e3t_1d(jk) = dk[ depw_1d ]
357      !!
358      !!            ===    Here constant vertical resolution   ===
359      !!
360      !! ** Action  : - pdept_1d, pdepw_1d : depth of T- and W-point (m)
361      !!              - pe3t_1d , pe3w_1d  : scale factors at T- and W-levels (m)
362      !!----------------------------------------------------------------------
363      REAL(wp), DIMENSION(:), INTENT(out) ::   pdept_1d, pdepw_1d   ! 1D grid-point depth        [m]
364      REAL(wp), DIMENSION(:), INTENT(out) ::   pe3t_1d , pe3w_1d    ! 1D vertical scale factors  [m]
365      !
366      INTEGER  ::   jk       ! dummy loop indices
367      REAL(wp) ::   zt, zw   ! local scalar
368      !!----------------------------------------------------------------------
369      !
370      IF(lwp) THEN                         ! Parameter print
371         WRITE(numout,*)
372         WRITE(numout,*) '    zgr_z1d : Reference vertical z-coordinates: uniform dz = ', rn_dz
373         WRITE(numout,*) '    ~~~~~~~'
374      ENDIF
375      !
376      ! Reference z-coordinate (depth - scale factor at T- and W-points)   ! Madec & Imbard 1996 function
377      ! ----------------------
378      DO jk = 1, jpk
379         zw = REAL( jk , wp )
380         zt = REAL( jk , wp ) + 0.5_wp
381         pdepw_1d(jk) =    rn_dz *   REAL( jk-1 , wp )
382         pdept_1d(jk) =    rn_dz * ( REAL( jk-1 , wp ) + 0.5_wp )
383         pe3w_1d (jk) =    rn_dz
384         pe3t_1d (jk) =    rn_dz
385      END DO
386      !
387      IF(lwp) THEN                        ! control print
388         WRITE(numout,*)
389         WRITE(numout,*) '              Reference 1D z-coordinate depth and scale factors:'
390         WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
391         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
392      ENDIF
393      !
394   END SUBROUTINE zgr_z1d
395   
396   !!======================================================================
397END MODULE usrdef_zgr
Note: See TracBrowser for help on using the repository browser.