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

Last change on this file since 9124 was 9124, checked in by gm, 6 years ago

dev_merge_2017: ln_timing instead of nn_timing + restricted timing to nemo_init and routine called by step in OPA_SRC

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