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

Last change on this file since 8977 was 8977, checked in by deazer, 7 years ago

Addresses several issues in the review except for rn_ssh_ref in the TEST cases

Builds and runs ok , little extra bracket dealt with in stpctl also

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         ! ocean space and time domain
20   USE usrdef_nam     ! User defined : namelist variables
21   USE wet_dry ,  ONLY: rn_wdmin1, rn_wdmin2, rn_wdld, &
22                      & rn_ht_0, rn_ssh_ref            ! Wetting and drying
23   !
24   USE in_out_manager ! I/O manager
25   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
26   USE lib_mpp        ! distributed memory computing library
27   USE wrk_nemo       ! Memory allocation
28   USE timing         ! Timing
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   usr_def_zgr   ! called by domzgr.F90
34
35  !! * Substitutions
36#  include "vectopt_loop_substitute.h90"
37   !!----------------------------------------------------------------------
38   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
39   !! $Id$
40   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42CONTAINS             
43
44   SUBROUTINE usr_def_zgr( ld_zco  , ld_zps  , ld_sco  , ld_isfcav,    &   ! type of vertical coordinate
45      &                    pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d  ,    &   ! 1D reference vertical coordinate
46      &                    pdept , pdepw ,                             &   ! 3D t & w-points depth
47      &                    pe3t  , pe3u  , pe3v , pe3f ,               &   ! vertical scale factors
48      &                    pe3w  , pe3uw , pe3vw,                      &   !     -      -      -
49      &                    k_top  , k_bot    )                             ! top & bottom ocean level
50      !!---------------------------------------------------------------------
51      !!              ***  ROUTINE usr_def_zgr  ***
52      !!
53      !! ** Purpose :   User defined the vertical coordinates
54      !!
55      !!----------------------------------------------------------------------
56      LOGICAL                   , INTENT(out) ::   ld_zco, ld_zps, ld_sco      ! vertical coordinate flags
57      LOGICAL                   , INTENT(out) ::   ld_isfcav                   ! under iceshelf cavity flag
58      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth     [m]
59      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d           ! 1D grid-point depth     [m]
60      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pdept, pdepw                ! grid-point depth        [m]
61      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors  [m]
62      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3w , pe3uw, pe3vw         ! i-scale factors
63      INTEGER , DIMENSION(:,:)  , INTENT(out) ::   k_top, k_bot                ! first & last ocean level
64      !
65      INTEGER  ::   ji, jj, jk        ! dummy indices
66      INTEGER  ::   ik                ! local integers
67      REAL(wp) ::   zfact, z1_jpkm1   ! local scalar
68      REAL(wp) ::   ze3min            ! local scalar
69      REAL(wp) ::   zi, zj, zbathy    ! local scalar
70      REAL(wp) ::   ztmpu, ztmpv, ztmpf, ztmpu1, ztmpv1, ztmpf1, zwet
71      REAL(wp), DIMENSION(jpi,jpj) ::   zht, zhu, zhv, z2d   ! 2D workspace
72      !!----------------------------------------------------------------------
73      !
74      IF(lwp) WRITE(numout,*)
75      IF(lwp) WRITE(numout,*) 'usr_def_zgr : WAD_TEST_CASES configuration (s-coordinate closed box ocean without cavities)'
76      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
77      !
78      !
79      ! type of vertical coordinate
80      ! ---------------------------
81      ! set in usrdef_nam.F90 by reading the namusr_def namelist only ln_zco
82      ln_zco    = .FALSE.      ! z-partial-step coordinate
83      ln_zps    = .FALSE.      ! z-partial-step coordinate
84      ln_sco    = .TRUE.       ! s-coordinate
85      ld_isfcav = .FALSE.      ! ISF Ice Shelves Flag
86      !
87      !
88      ! Build the vertical coordinate system
89      ! ------------------------------------
90      !
91      !                       !==  UNmasked meter bathymetry  ==!
92      !
93!
94      zbathy=10.0
95      IF( cn_cfg == 'wad' ) THEN
96         SELECT CASE ( nn_cfg )
97            !                                        ! ====================
98            CASE ( 1 )                               ! WAD 1 configuration
99               !                                     ! ====================
100               !
101               IF(lwp) WRITE(numout,*)
102               IF(lwp) WRITE(numout,*) 'usr_def_zgr (WAD) : Closed box with EW linear bottom slope'
103               IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
104               !
105               zht = 1.5_wp
106               DO ji = 1, jpi
107                 zi = MIN((glamt(ji,1) - 10.0)/40.0, 1.0 )
108                 zht(ji,:) = MAX(zbathy*zi, -2.0) 
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, 8 )                         ! 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.3 )
124                 zht(ji,:) = MAX(zbathy*zi, -2.0)
125               END DO
126               zht(mi0(1):mi1(1),:) = -4._wp
127               zht(mi0(jpiglo):mi1(jpiglo),:) = -4._wp
128               IF( nn_cfg /= 8 ) THEN
129                  zht(:,mj0(1):mj1(1)) = -4._wp
130                  zht(:,mj0(jpjglo):mj1(jpjglo)) = -4._wp
131               ENDIF
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               END DO
147               zht(mi0(1):mi1(1),:) = -4._wp
148               zht(mi0(jpiglo):mi1(jpiglo),:) = -4._wp
149               zht(:,mj0(1):mj1(1)) = -4._wp
150               zht(:,mj0(jpjglo):mj1(jpjglo)) = -4._wp
151               !                                    ! ===========================
152            CASE ( 5 )                              ! WAD 5 configuration
153               !                                    ! ====================
154               !
155               IF(lwp) WRITE(numout,*)
156               IF(lwp) WRITE(numout,*) 'usr_def_zgr (WAD) : Double slope with shelf'
157               IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
158               !
159               DO ji = 1, jpi
160                 zi = MIN(glamt(ji,1)/45.0, 1.0 )
161                 zht(ji,:) = MAX(zbathy*zi, -2.0)
162                 IF( glamt(ji,1) >= 46.0 ) THEN
163                   zht(ji,:) = 10.0
164                 ELSE IF( glamt(ji,1) >= 20.0 .AND. glamt(ji,1) < 46.0 ) THEN
165                   zi = 7.5/25.
166                   zht(ji,:) = MAX(10. - zi*(47.-glamt(ji,1)),2.5)
167                 ELSE IF( glamt(ji,1) >= 15.0 .AND. glamt(ji,1) < 20.0 ) THEN
168                   zht(ji,:) = 2.5
169                 ELSE IF( glamt(ji,1) >= 4.0 .AND. glamt(ji,1) < 15.0 ) THEN
170                   zi = 4.5/11.0
171                   zht(ji,:) = MAX(2.5 - zi*(16.0-glamt(ji,1)), -2.0)
172                 ELSE IF( glamt(ji,1) >= 0.0 .AND. glamt(ji,1) < 4.0 ) THEN
173                   zht(ji,:) = -2.0
174                 ENDIF
175               END DO
176               !                                    ! ===========================
177               zht(mi0(1):mi1(1),:) = -4._wp
178               zht(mi0(jpiglo):mi1(jpiglo),:) = -4._wp
179               zht(:,mj0(1):mj1(1)) = -4._wp
180               zht(:,mj0(jpjglo):mj1(jpjglo)) = -4._wp
181               !                                    ! ===========================
182            CASE ( 6 )                              ! WAD 6 configuration
183               !                                    ! ====================
184               !
185               IF(lwp) WRITE(numout,*)
186               IF(lwp) WRITE(numout,*) 'usr_def_zgr (WAD) : Parabolic channel with gaussian ridge'
187               IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
188               !
189               DO ji = 1, jpi
190                 zi = MAX(1.0-((glamt(ji,1)-25.0)**2)/484.0, -2.0 )
191                 zj = 1.075*MAX(EXP(-1.0*((glamt(ji,1)-25.0)**2)/32.0) , 0.0 )
192                 zht(ji,:) = MAX(zbathy*(zi-zj), -2.0)
193               END DO
194               zht(mi0(1):mi1(1),:) = -4._wp
195               zht(mi0(jpiglo):mi1(jpiglo),:) = -4._wp
196               zht(:,mj0(1):mj1(1)) = -4._wp
197               zht(:,mj0(jpjglo):mj1(jpjglo)) = -4._wp
198               !                                    ! ===========================
199            CASE ( 7 )                              ! WAD 7 configuration
200               !                                    ! ====================
201               !
202               IF(lwp) WRITE(numout,*)
203               IF(lwp) WRITE(numout,*) 'usr_def_zgr (WAD) : Double slope with open boundary'
204               IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
205               !
206               DO ji = 1, jpi
207                 zi = MIN(glamt(ji,1)/45.0, 1.0 )
208                 zht(ji,:) = MAX(zbathy*zi, -2.0)
209                 IF( glamt(ji,1) >= 46.0 ) THEN
210                   zht(ji,:) = 10.0
211                 ELSE IF( glamt(ji,1) >= 20.0 .AND. glamt(ji,1) < 46.0 ) THEN
212                   zi = 7.5/25.
213                   zht(ji,:) = MAX(10. - zi*(47.-glamt(ji,1)),2.5)
214                 ELSE IF( glamt(ji,1) >= 15.0 .AND. glamt(ji,1) < 20.0 ) THEN
215                   zht(ji,:) = 2.5
216                 ELSE IF( glamt(ji,1) >= 4.0 .AND. glamt(ji,1) < 15.0 ) THEN
217                   zi = 4.5/11.0
218                   zht(ji,:) = MAX(2.5 - zi*(16.0-glamt(ji,1)), -2.0)
219                 ELSE IF( glamt(ji,1) >= 0.0 .AND. glamt(ji,1) < 4.0 ) THEN
220                   zht(ji,:) = -2.0
221                 ENDIF
222               END DO
223               !                                    ! ===========================
224               zht(mi0(1):mi1(1),:) = -4._wp
225               zht(:,mj0(1):mj1(1)) = -4._wp
226               zht(:,mj0(jpjglo):mj1(jpjglo)) = -4._wp
227            CASE DEFAULT
228               !                                    ! ===========================
229               WRITE(ctmp1,*) 'WAD test with a ', nn_cfg,' option is not coded'
230               !
231               CALL ctl_stop( ctmp1 )
232               !
233         END SELECT
234      END IF
235             
236! increase the depth of the bathymetry by rn_ssh_ref and 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.