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 NEMO/trunk/tests/WAD/MY_SRC – NEMO

source: NEMO/trunk/tests/WAD/MY_SRC/usrdef_zgr.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

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