Newer
Older
MODULE usrdef_zgr
!!======================================================================
!! *** MODULE usrdef_zgr ***
!!
!! === WAD_TEST_CASES case ===
!!
!! Ocean domain : user defined vertical coordinate system
!!======================================================================
!! History : 4.0 ! 2016-06 (G. Madec) Original code
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! usr_def_zgr : user defined vertical coordinate system (required)
!! zgr_z : reference 1D z-coordinate
!!---------------------------------------------------------------------
USE oce ! ocean variables
USE dom_oce ! ocean space and time domain
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
USE usrdef_nam ! User defined : namelist variables
USE wet_dry , ONLY: rn_wdmin1, rn_wdmin2, rn_wdld ! Wetting and drying
!
USE in_out_manager ! I/O manager
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
USE lib_mpp ! distributed memory computing library
IMPLICIT NONE
PRIVATE
PUBLIC usr_def_zgr ! called by domzgr.F90
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: usrdef_zgr.F90 15033 2021-06-21 10:24:45Z smasson $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE usr_def_zgr( ld_zco , ld_zps , ld_sco , ld_isfcav, & ! type of vertical coordinate
& pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d , & ! 1D reference vertical coordinate
& pdept , pdepw , & ! 3D t & w-points depth
& pe3t , pe3u , pe3v , pe3f , & ! vertical scale factors
& pe3w , pe3uw , pe3vw , & ! - - -
& k_top , k_bot ) ! top & bottom ocean level
!!---------------------------------------------------------------------
!! *** ROUTINE usr_def_zgr ***
!!
!! ** Purpose : User defined the vertical coordinates
!!
!!----------------------------------------------------------------------
LOGICAL , INTENT(out) :: ld_zco, ld_zps, ld_sco ! vertical coordinate flags
LOGICAL , INTENT(out) :: ld_isfcav ! under iceshelf cavity flag
REAL(wp), DIMENSION(:) , INTENT(out) :: pdept_1d, pdepw_1d ! 1D grid-point depth [m]
REAL(wp), DIMENSION(:) , INTENT(out) :: pe3t_1d , pe3w_1d ! 1D grid-point depth [m]
REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pdept, pdepw ! grid-point depth [m]
REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pe3t , pe3u , pe3v , pe3f ! vertical scale factors [m]
REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pe3w , pe3uw, pe3vw ! i-scale factors
INTEGER , DIMENSION(:,:) , INTENT(out) :: k_top, k_bot ! first & last ocean level
!
INTEGER :: ji, jj, jk ! dummy indices
INTEGER :: ik ! local integers
REAL(wp) :: zfact, z1_jpkm1 ! local scalar
REAL(wp) :: ze3min ! local scalar
REAL(wp) :: zi, zj, zbathy ! local scalar
REAL(wp) :: ztmpu, ztmpv, ztmpf, ztmpu1, ztmpv1, ztmpf1, zwet
REAL(wp), DIMENSION(jpi,jpj) :: zht, zhu, zhv, z2d ! 2D workspace
!!----------------------------------------------------------------------
!
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'usr_def_zgr : WAD_TEST_CASES configuration (s-coordinate closed box ocean without cavities)'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
!
!
! type of vertical coordinate ==>>> here WAD_TEST_CASES : s-coordinate always
! ---------------------------
ld_zco = .FALSE. ! z-partial-step coordinate
ld_zps = .FALSE. ! z-partial-step coordinate
ld_sco = .TRUE. ! s-coordinate
ld_isfcav = .FALSE. ! ISF Ice Shelves Flag
!
!
! Build the vertical coordinate system
! ------------------------------------
!
! !== UNmasked meter bathymetry ==!
!
!
zbathy=10.0
IF( cn_cfg == 'wad' ) THEN
SELECT CASE ( nn_cfg )
! ! ====================
CASE ( 1 ) ! WAD 1 configuration
! ! ====================
!
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'usr_def_zgr (WAD) : Closed box with EW linear bottom slope'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
!
zht = 1.5_wp
DO ji = 1, jpi
zi = MIN((glamt(ji,1) - 10.0)/40.0, 1.0 )
zht(ji,:) = MAX(zbathy*zi, -2.0)
END DO
zht(mi0( 1,nn_hls):mi1( 1,nn_hls),:) = -4._wp
zht(mi0(jpiglo,nn_hls):mi1(jpiglo,nn_hls),:) = -4._wp
zht(:,mj0( 1,nn_hls):mj1( 1,nn_hls)) = -4._wp
zht(:,mj0(jpjglo,nn_hls):mj1(jpjglo,nn_hls)) = -4._wp
! ! ====================
CASE ( 2, 3, 8 ) ! WAD 2 or 3 configuration
! ! ====================
!
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'usr_def_zgr (WAD) : Parobolic EW channel'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
!
DO ji = 1, jpi
zi = MAX(1.0-((glamt(ji,1)-25.0)**2)/484.0, -0.3 )
zht(ji,:) = MAX(zbathy*zi, -2.0)
END DO
zht(mi0( 1,nn_hls):mi1( 1,nn_hls),:) = -4._wp
zht(mi0(jpiglo,nn_hls):mi1(jpiglo,nn_hls),:) = -4._wp
zht(:,mj0( 1,nn_hls):mj1( 1,nn_hls)) = -4._wp
zht(:,mj0(jpjglo,nn_hls):mj1(jpjglo,nn_hls)) = -4._wp
ENDIF
! ! ====================
CASE ( 4 ) ! WAD 4 configuration
! ! ====================
!
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'usr_def_zgr (WAD) : Parobolic bowl'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
!
DO ji = 1, jpi
zi = MAX(1.0-((glamt(ji,1)-25.0)**2)/484.0, -2.0 )
DO jj = 1, jpj
zj = MAX(1.0-((gphit(1,jj)-17.0)**2)/196.0, -2.0 )
zht(ji,jj) = MAX(zbathy*zi*zj, -2.0)
END DO
END DO
zht(mi0(1 ,nn_hls):mi1( 1,nn_hls),:) = -4._wp
zht(mi0(jpiglo,nn_hls):mi1(jpiglo,nn_hls),:) = -4._wp
zht(:,mj0( 1,nn_hls):mj1( 1,nn_hls)) = -4._wp
zht(:,mj0(jpjglo,nn_hls):mj1(jpjglo,nn_hls)) = -4._wp
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
! ! ===========================
CASE ( 5 ) ! WAD 5 configuration
! ! ====================
!
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'usr_def_zgr (WAD) : Double slope with shelf'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
!
DO ji = 1, jpi
zi = MIN(glamt(ji,1)/45.0, 1.0 )
zht(ji,:) = MAX(zbathy*zi, -2.0)
IF( glamt(ji,1) >= 46.0 ) THEN
zht(ji,:) = 10.0
ELSE IF( glamt(ji,1) >= 20.0 .AND. glamt(ji,1) < 46.0 ) THEN
zi = 7.5/25.
zht(ji,:) = MAX(10. - zi*(47.-glamt(ji,1)),2.5)
ELSE IF( glamt(ji,1) >= 15.0 .AND. glamt(ji,1) < 20.0 ) THEN
zht(ji,:) = 2.5
ELSE IF( glamt(ji,1) >= 4.0 .AND. glamt(ji,1) < 15.0 ) THEN
zi = 4.5/11.0
zht(ji,:) = MAX(2.5 - zi*(16.0-glamt(ji,1)), -2.0)
ELSE IF( glamt(ji,1) >= 0.0 .AND. glamt(ji,1) < 4.0 ) THEN
zht(ji,:) = -2.0
ENDIF
END DO
! ! ===========================
zht(mi0( 1,nn_hls):mi1( 1,nn_hls),:) = -4._wp
zht(mi0(jpiglo,nn_hls):mi1(jpiglo,nn_hls),:) = -4._wp
zht(:,mj0( 1,nn_hls):mj1( 1,nn_hls)) = -4._wp
zht(:,mj0(jpjglo,nn_hls):mj1(jpjglo,nn_hls)) = -4._wp
! ! ===========================
CASE ( 6 ) ! WAD 6 configuration
! ! ====================
!
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'usr_def_zgr (WAD) : Parabolic channel with gaussian ridge'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
!
DO ji = 1, jpi
zi = MAX(1.0-((glamt(ji,1)-25.0)**2)/484.0, -2.0 )
zj = 1.075*MAX(EXP(-1.0*((glamt(ji,1)-25.0)**2)/32.0) , 0.0 )
zht(ji,:) = MAX(zbathy*(zi-zj), -2.0)
END DO
zht(mi0( 1,nn_hls):mi1( 1,nn_hls),:) = -4._wp
zht(mi0(jpiglo,nn_hls):mi1(jpiglo,nn_hls),:) = -4._wp
zht(:,mj0( 1,nn_hls):mj1( 1,nn_hls)) = -4._wp
zht(:,mj0(jpjglo,nn_hls):mj1(jpjglo,nn_hls)) = -4._wp
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
! ! ===========================
CASE ( 7 ) ! WAD 7 configuration
! ! ====================
!
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'usr_def_zgr (WAD) : Double slope with open boundary'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
!
DO ji = 1, jpi
zi = MIN(glamt(ji,1)/45.0, 1.0 )
zht(ji,:) = MAX(zbathy*zi, -2.0)
IF( glamt(ji,1) >= 46.0 ) THEN
zht(ji,:) = 10.0
ELSE IF( glamt(ji,1) >= 20.0 .AND. glamt(ji,1) < 46.0 ) THEN
zi = 7.5/25.
zht(ji,:) = MAX(10. - zi*(47.-glamt(ji,1)),2.5)
ELSE IF( glamt(ji,1) >= 15.0 .AND. glamt(ji,1) < 20.0 ) THEN
zht(ji,:) = 2.5
ELSE IF( glamt(ji,1) >= 4.0 .AND. glamt(ji,1) < 15.0 ) THEN
zi = 4.5/11.0
zht(ji,:) = MAX(2.5 - zi*(16.0-glamt(ji,1)), -2.0)
ELSE IF( glamt(ji,1) >= 0.0 .AND. glamt(ji,1) < 4.0 ) THEN
zht(ji,:) = -2.0
ENDIF
END DO
! ! ===========================
zht(mi0( 1,nn_hls):mi1( 1,nn_hls),:) = -4._wp
zht(:,mj0( 1,nn_hls):mj1( 1,nn_hls)) = -4._wp
zht(:,mj0(jpjglo,nn_hls):mj1(jpjglo,nn_hls)) = -4._wp
CASE DEFAULT
! ! ===========================
WRITE(ctmp1,*) 'WAD test with a ', nn_cfg,' option is not coded'
!
CALL ctl_stop( ctmp1 )
!
END SELECT
END IF
! at u-point: averaging zht
DO_2D( 0, 0, 0, 0 )
zhu(ji,jj) = 0.5_wp * ( zht(ji,jj) + zht(ji+1,jj) )
END_2D
CALL lbc_lnk( 'usrdef_zgr', zhu, 'U', 1. ) ! boundary condition: this mask the surrounding grid-points
! ! ==>>> set by hand non-zero value on first/last columns & rows
DO ji = mi0( 1,nn_hls), mi1( 1,nn_hls) ! first row of global domain only
DO ji = mi0(jpiglo,nn_hls), mi1(jpiglo,nn_hls) ! last row of global domain only
zhu(ji,:) = zht(jpi,:)
END DO
! at v-point: averaging zht
zhv = 0._wp
DO_2D( 0, 0, 0, 0 )
zhv(ji,jj) = 0.5_wp * ( zht(ji,jj) + zht(ji,jj+1) )
END_2D
CALL lbc_lnk( 'usrdef_zgr', zhv, 'V', 1. ) ! boundary condition: this mask the surrounding grid-points
DO jj = mj0( 1,nn_hls), mj1( 1,nn_hls) ! first row of global domain only
DO jj = mj0(jpjglo,nn_hls), mj1(jpjglo,nn_hls) ! last row of global domain only
zhv(:,jj) = zht(:,jj)
END DO
!
CALL zgr_z( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d ) ! Reference z-coordinate system
!
!
! !== top masked level bathymetry ==! (all coordinates)
!
! no ocean cavities : top ocean level is ONE, except over land
! the ocean basin surrounnded by land (1+nn_hls grid-points) set through lbc_lnk call
z2d(:,:) = 1._wp ! surface ocean is the 1st level
z2d(mi0( 1,nn_hls):mi1( 1,nn_hls),:) = 0._wp
z2d(mi0(jpiglo,nn_hls):mi1(jpiglo,nn_hls),:) = 0._wp
z2d(:,mj0( 1,nn_hls):mj1( 1,nn_hls)) = 0._wp
z2d(:,mj0(jpjglo,nn_hls):mj1(jpjglo,nn_hls)) = 0._wp
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. ) ! closed basin, see userdef_nam.F90
k_top(:,:) = NINT( z2d(:,:) )
!
!
!
IF ( ld_sco ) THEN !== s-coordinate ==! (terrain-following coordinate)
!
ht_0 = zht
k_bot(:,:) = jpkm1 * k_top(:,:) !* bottom ocean = jpk-1 (here use k_top as a land mask)
DO_2D( 1, 1, 1, 1 )
IF( zht(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN
k_bot(ji,jj) = 0
k_top(ji,jj) = 0
ENDIF
END_2D
!
! !* terrain-following coordinate with e3.(k)=cst)
! ! OVERFLOW case : identical with j-index (T=V, U=F)
DO_2D( 0, 0, 0, 0 )
z1_jpkm1 = 1._wp / REAL( k_bot(ji,jj) - k_top(ji,jj) + 1 , wp)
DO jk = 1, jpk
zwet = MAX( zht(ji,jj), rn_wdmin1 )
pdept(ji,jj,jk) = zwet * z1_jpkm1 * ( REAL( jk , wp ) - 0.5_wp )
pdepw(ji,jj,jk) = zwet * z1_jpkm1 * ( REAL( jk-1 , wp ) )
pe3t (ji,jj,jk) = zwet * z1_jpkm1
pe3w (ji,jj,jk) = zwet * z1_jpkm1
zwet = MAX( zhu(ji,jj), rn_wdmin1 )
pe3u (ji,jj,jk) = zwet * z1_jpkm1
pe3uw(ji,jj,jk) = zwet * z1_jpkm1
pe3f (ji,jj,jk) = zwet * z1_jpkm1
zwet = MAX( zhv(ji,jj), rn_wdmin1 )
pe3v (ji,jj,jk) = zwet * z1_jpkm1
pe3vw(ji,jj,jk) = zwet * z1_jpkm1
END DO
END_2D
CALL lbc_lnk( 'usrdef_zgr', pdept, 'T', 1. )
CALL lbc_lnk( 'usrdef_zgr', pdepw, 'T', 1. )
CALL lbc_lnk( 'usrdef_zgr', pe3t , 'T', 1. )
CALL lbc_lnk( 'usrdef_zgr', pe3w , 'T', 1. )
CALL lbc_lnk( 'usrdef_zgr', pe3u , 'U', 1. )
CALL lbc_lnk( 'usrdef_zgr', pe3uw, 'U', 1. )
CALL lbc_lnk( 'usrdef_zgr', pe3f , 'F', 1. )
CALL lbc_lnk( 'usrdef_zgr', pe3v , 'V', 1. )
CALL lbc_lnk( 'usrdef_zgr', pe3vw, 'V', 1. )
WHERE( pe3t (:,:,:) == 0._wp ) pe3t (:,:,:) = 1._wp
WHERE( pe3u (:,:,:) == 0._wp ) pe3u (:,:,:) = 1._wp
WHERE( pe3v (:,:,:) == 0._wp ) pe3v (:,:,:) = 1._wp
WHERE( pe3f (:,:,:) == 0._wp ) pe3f (:,:,:) = 1._wp
WHERE( pe3w (:,:,:) == 0._wp ) pe3w (:,:,:) = 1._wp
WHERE( pe3uw(:,:,:) == 0._wp ) pe3uw(:,:,:) = 1._wp
WHERE( pe3vw(:,:,:) == 0._wp ) pe3vw(:,:,:) = 1._wp
ENDIF
!
!
!
END SUBROUTINE usr_def_zgr
SUBROUTINE zgr_z( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d ) ! 1D reference vertical coordinate
!!----------------------------------------------------------------------
!! *** ROUTINE zgr_z ***
!!
!! ** Purpose : set the depth of model levels and the resulting
!! vertical scale factors.
!!
!! ** Method : 1D z-coordinate system (use in all type of coordinate)
!! The depth of model levels is set from dep(k), an analytical function:
!! w-level: depw_1d = dep(k)
!! t-level: dept_1d = dep(k+0.5)
!! The scale factors are the discrete derivative of the depth:
!! e3w_1d(jk) = dk[ dept_1d ]
!! e3t_1d(jk) = dk[ depw_1d ]
!!
!! === Here constant vertical resolution ===
!!
!! ** Action : - pdept_1d, pdepw_1d : depth of T- and W-point (m)
!! - pe3t_1d , pe3w_1d : scale factors at T- and W-levels (m)
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(:), INTENT(out) :: pdept_1d, pdepw_1d ! 1D grid-point depth [m]
REAL(wp), DIMENSION(:), INTENT(out) :: pe3t_1d , pe3w_1d ! 1D vertical scale factors [m]
!
INTEGER :: jk ! dummy loop indices
REAL(wp) :: zt, zw ! local scalar
!!----------------------------------------------------------------------
!
IF(lwp) THEN ! Parameter print
WRITE(numout,*)
WRITE(numout,*) ' zgr_z : Reference vertical z-coordinates: uniform dz = ', rn_dz
WRITE(numout,*) ' ~~~~~~~'
ENDIF
!
! Reference z-coordinate (depth - scale factor at T- and W-points) ! Madec & Imbard 1996 function
! ----------------------
DO jk = 1, jpk
zw = REAL( jk , wp )
zt = REAL( jk , wp ) + 0.5_wp
pdepw_1d(jk) = rn_dz * REAL( jk-1 , wp )
pdept_1d(jk) = rn_dz * ( REAL( jk-1 , wp ) + 0.5_wp )
pe3w_1d (jk) = rn_dz
pe3t_1d (jk) = rn_dz
END DO
!
IF(lwp) THEN ! control print
WRITE(numout,*)
WRITE(numout,*) ' Reference 1D z-coordinate depth and scale factors:'
WRITE(numout, "(9x,' level gdept_1d gdepw_1d e3t_1d e3w_1d ')" )
WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
ENDIF
!
END SUBROUTINE zgr_z
!!======================================================================
END MODULE usrdef_zgr