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/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/tests/WAD/MY_SRC – NEMO

source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/tests/WAD/MY_SRC/usrdef_zgr.F90 @ 12939

Last change on this file since 12939 was 12939, checked in by smasson, 4 years ago

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