Context Navigation

usrdef_zgr.F90 @ 8403

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

Add in ROMS WAD option ln_rwd+changes for implicit Bed Friction for ln_wd option
Note no ramp placed on ROMS bed friction yet
CS15mini case added as a Test CASE
at this revision AMM15 with Pure sigma coords barotorpic runs for 4 days without failure
in with ROMS option with 20cm min deoth and 50 vertical levels
Both run for CS15mini
In real domains nothing done on reference level yet so real domains
must have not negative depth points yet.
But a basic test has been done in WAD channel test cases (WAD7)

No changes in Main line source yet. See the MY_SRC sub dir of CS15 and TEST_CASES/WAD
for actual code changes.

File size: 18.6 KB

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

Note: See TracBrowser for help on using the repository browser.

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/ROMS_WAD_7832/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_SRC/usrdef_zgr.F90 @ 8403

Download in other formats: