Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

PHY_vdfTKE_RUN.f90 @ 4031

Last change on this file since 4031 was 3990, checked in by millour, 8 years ago
An updated LMDZ5 (based on LMDZ rev 2786) to play with. EM
File size: 10.2 KB

Line
1	subroutine PHY_vdfTKE_RUN
2
3	!------------------------------------------------------------------------------+
4	! Mon 17-Jun-2013 MAR \|
5	! SubRoutine PHY_vdfTKE_RUN includes Vertical Diffusion of TKE \|
6	! and epsilon \|
7	! \|
8	! \|
9	! version 3.p.4.1 created by H. Gallee, Tue 19-Mar-2013 \|
10	! Last Modification by H. Gallee, Mon 17-Jun-2013 \|
11	! \|
12	!------------------------------------------------------------------------------+
13	! \|
14	! INPUT: Kzm_AT Vertical Turbulent Coeffic.(momentum) [m2/s2] \|
15	! ^^^^^^ \|
16	! \|
17	! INPUT / OUTPUT: The Vertical Turbulent Fluxes are included for: \|
18	! ^^^^^^^^^^^^^^ \|
19	! a) Turbulent Kinetic Energy TKE_AT(_xyz) [m2/s2] \|
20	! b) Turbulent Kinetic Energy Dissipation eps_AT(_xyz) [m2/s3] \|
21	! \|
22	! #OPTIONS: #De: Dirichlet Type Top Boundary Condit. for TKE_AT (TKE ) \|
23	! #^^^^^^^^ & eps_AT (epsilon) \|
24	! \|
25	!------------------------------------------------------------------------------+
26
27	use Mod_Real
28	use Mod_PHY____dat
29	use Mod_PHY____grd
30	use Mod_PHY_AT_grd
31	use Mod_PHY_AT_kkl
32	use Mod_PHY_DY_kkl
33
34
35
36	! Local Variables
37	! ================
38
39	use Mod_vdfTKE_RUN
40
41
42	IMPLICIT NONE
43
44
45	real(kind=real8) :: S3DSBC
46	real(kind=real8) :: sige2k
47	real(kind=real8) :: psa_sq
48	! #De real(kind=real8) :: TKEtop = 0.
49
50	integer :: i ,j ,k
51	integer :: k1 ,ikl
52
53
54
55
56	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
57	! !
58	! ALLOCATION
59	! ==========
60
61	IF (it_RUN.EQ.1 .OR. FlagDALLOC) THEN !
62	allocate ( S3D__1(mzp-1) )
63	allocate ( S3D__2(mzp-1) )
64	allocate ( S3D__3(mzp-1) )
65	allocate ( S3D__4(mzp-1) )
66	allocate ( S3D__5(mzp-1) )
67	allocate ( S3D__6(mzp-1) )
68	allocate ( S3D__7(mzp-1) )
69	END IF !
70	! !
71	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
72
73
74
75
76	! INITIALIZATION
77	! ==============
78
79	sige2k = sige / sigk
80
81
82
83	! =====================
84	DO ikl = 1,kcolp
85	! =====================
86
87	i = ii__AP(ikl)
88	j = jj__AP(ikl)
89
90	psa_sq = psa_DY(ikl) *psa_DY(ikl)
91
92
93
94
95	! Vertical Diffusion of Turbulent Kinetic Energy
96	! ==============================================
97
98
99	! Tridiagonal Matrix Coefficients - TKE_AT
100	! ----------------------------------------
101
102	k=mzp-1
103	S3DSBC = -GravF20.5(Kzm_AT(ikl,k)+Kzm_AT(ikl,k+1)) &! SBC: TKE and epsilon, Atm SBL
104	& * sigk betaAT roamDY(ikl,k)*roa_DY(ikl,k ) &!
105	& /(psa_sq * dsigmi(k) *dsigma( k+1))
106	S3D__1(mzp-1) = S3DSBC ! SBC: TKE and epsilon, Atm SBL
107
108	DO k=mzp-2,1,-1
109	S3D__1(k) =-GravF20.5(Kzm_AT(ikl,k)+Kzm_AT(ikl,k+1)) &
110	& * sigk betaAT roamDY(ikl,k)*roa_DY(ikl,k ) &
111	& /(psa_sq * dsigmi(k) *dsigma( k+1))
112
113	S3D__3(k+1) = S3D__1(k) * dsigmi( k)/dsigmi( k+1) &
114	& / roa_DY(ikl,k)*roa_DY(ikl,k+1)
115	END DO
116
117	S3D__3(1) = 0.0 ! UBC: Von Neuman , Atm Top
118	DO k= 1,mzp-1
119	S3D__1(k) = S3D__1(k) * dt__AT
120	S3D__3(k) = S3D__3(k) * dt__AT
121	S3D__2(k) = 1.0 -S3D__3(k) -S3D__1(k)
122	END DO
123
124
125	! Second Member of the Tridiagonal System - TKE_AT
126	! ------------------------------------------------
127
128	S3D__4(1) = &
129	& S3D__1(1) a_b_AT(TKE_AT(ikl,1)-TKE_AT(ikl,k1p(1)))
130	! #De S3D__1(1) = 0.0
131	! #De S3D__2(1) = 1.0
132	! #De S3D__4(1) = TKEtop
133
134	DO k=k1p(1),mzp-2
135	S3D__4(k) = &
136	& S3D__1(k) a_b_AT(TKE_AT(ikl, k )-TKE_AT(ikl,k1p(1))) &
137	& -S3D__3(k) a_b_AT(TKE_AT(ikl,k1m(k))-TKE_AT(ikl, k ))
138	END DO
139
140	k= mzp-1
141	S3D__4(k) = &
142	& S3D__1(k)(a_b_AT TKE_AT(ikl, k )-TKE_AT(ikl,k1p(1)) &
143	& /betaAT ) &
144	& -S3D__3(k) a_b_AT(TKE_AT(ikl,k1m(k))-TKE_AT(ikl, k ))
145
146	S3D__1(k) = 0.000
147
148	! S3D__4(mzp-1)=-(alphAT* TKE_AT(ikl,mzp-1)-TKE_AT(ikl,mzp )) &
149	! & * GravF2 0.5000(Kzm_AT(ikl,mzp-1)+Kzm_AT(ikl,mzp-2)) &
150	! & * roamDY(ikl,mzp-1)*roamDY(ikl,mzp-1) &
151	! & / (psa_sq* dsigmi( mzp-1)*dsigma( mzp-1)) &
152	! & -S3D__3(mzp-1)* a_b_AT*(TKE_AT(ikl,mzp-2)-TKE_AT(ikl,mzp-1))
153
154
155	! Tridiagonal Matrix Inversion - TKE_AT
156	! -------------------------------------
157
158	k1= 1
159	! #De k1= 2
160	DO k=k1,mzp-1
161	S3D__4(k) = S3D__4(k) + TKE_AT(ikl,k)
162	END DO
163
164	! Forward Sweep
165	! ~~~~~~~~~~~~~~
166	S3D__5(1) = S3D__2(1)
167	S3D__6(1) =-S3D__1(1) /S3D__5(1)
168	DO k=k1p(1),mzp-1
169	S3D__5(k) = S3D__3(k) *S3D__6(k-1)+S3D__2(k)
170	S3D__6(k) =-S3D__1(k) /S3D__5(k)
171	END DO
172	S3D__7(1) = S3D__4(1) /S3D__5(1)
173	DO k=k1p(1),mzp-1
174	S3D__7(k) =(S3D__4(k) -S3D__3(k) &
175	& *S3D__7(k-1))/S3D__5(k)
176	END DO
177
178	! Backward Sweep
179	! ~~~~~~~~~~~~~~
180	DO k=k1m(mzp-1),1,-1
181	S3D__7(k) = S3D__6(k) *S3D__7(k+1)+S3D__7(k)
182	END DO
183
184
185	DO k=1,mzp-1
186	TrT_AT(ikl,k) = TrT_AT(ikl,k) &
187	& +(S3D__7(k) -TKE_AT(ikl,k)) /dt__AT
188	TKE_AT(ikl,k) = S3D__7(k)
189	END DO
190
191
192	! Vertical Diffusion of Dissipation
193	! =================================
194
195
196	! Update Tridiagonal Matrix Coefficients - eps_AT
197	! -----------------------------------------------
198
199	S3D__1(mzp-1) = S3DSBC * dt__AT ! SBC: TKE and epsilon, Atm SBL
200	DO k=1,mzp-1
201	S3D__1(k) = S3D__1(k) * sige2k
202	S3D__3(k) = S3D__3(k) * sige2k
203	S3D__2(k) = 1.0 - S3D__3(k) - S3D__1(k)
204	END DO
205
206
207	! Second Member of the Tridiagonal System - eps_AT
208	! ------------------------------------------------
209
210	S3D__4(1) = &
211	& S3D__1(1) a_b_AT(eps_AT(ikl,1)-eps_AT(ikl,k1p(1)))
212	! #De S3D__1(1) = 0.0
213	! #De S3D__2(1) = 1.0
214	! #De S3D__4(1) = eps_DI(i,j)
215
216	DO k=k1p(1),mzp-2
217	S3D__4(k) = &
218	& S3D__1(k) a_b_AT(eps_AT(ikl,k)-eps_AT(ikl,k1p(k))) &
219	& -S3D__3(k) a_b_AT(eps_AT(ikl,k1m(k))-eps_AT(ikl,k))
220	END DO
221
222	k= mzp-1
223	S3D__4(k) = &
224	& S3D__1(k)(a_b_AT eps_AT(ikl, k )-eps_AT(ikl,k1p(1)) &
225	& /betaAT ) &
226	& -S3D__3(k) a_b_AT(eps_AT(ikl,k1m(k))-eps_AT(ikl, k ))
227
228	S3D__1(k) = 0.000
229
230	! S3D__4(mzp-1)=-(alphAT* eps_AT(ikl,mzp-1)-eps_AT(ikl,mzp)) &
231	! & * GravF2* 0.5000*(Kzm_AT(ikl,mzp-1)+Kzm_AT(ikl,mzp-2)) &
232	! & * roamDY(ikl,mzp-1)*roamDY(ikl,mzp-1) &
233	! & / (psa_sq* dsigmi( mzp-1)*dsigma( mzp-1)) &
234	! & -S3D__3(mzp-1)* a_b_AT*(eps_AT(ikl,mzp-2)-eps_AT(ikl,mzp-1))
235
236
237	! Tridiagonal Matrix Inversion - eps_AT
238	! -------------------------------------
239
240	k1= 1
241	! #De k1= 2
242	DO k=k1,mzp-1
243	S3D__4(k) = S3D__4(k) + eps_AT(ikl,k)
244	END DO
245
246	! Forward Sweep
247	! ~~~~~~~~~~~~~~
248	S3D__5(1) = S3D__2(1)
249	S3D__6(1) =-S3D__1(1) /S3D__5(1)
250	DO k=k1p(1),mzp-1
251	S3D__5(k) = S3D__3(k) *S3D__6(k-1)+S3D__2(k)
252	S3D__6(k) =-S3D__1(k) /S3D__5(k)
253	END DO
254	S3D__7(1) = S3D__4(1) /S3D__5(1)
255	DO k=k1p(1),mzp-1
256	S3D__7(k) =(S3D__4(k) -S3D__3(k) &
257	& *S3D__7(k-1))/S3D__5(k)
258	END DO
259
260	! Backward Sweep
261	! ~~~~~~~~~~~~~~
262	DO k=k1m(mzp-1),1,-1
263	S3D__7(k) = S3D__6(k) *S3D__7(k+1)+S3D__7(k)
264	END DO
265
266	DO k=1,mzp-1
267	eps_AT(ikl,k) = S3D__7(k)
268	END DO
269
270
271
272
273	! =====================
274	ENDDO ! ikl = 1,kcolp
275	! =====================
276
277
278
279
280	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
281	! !
282	! DE-ALLOCATION
283	! =============
284
285	IF (FlagDALLOC) THEN !
286	deallocate ( S3D__1 )
287	deallocate ( S3D__2 )
288	deallocate ( S3D__3 )
289	deallocate ( S3D__4 )
290	deallocate ( S3D__5 )
291	deallocate ( S3D__6 )
292	deallocate ( S3D__7 )
293	ENDIF
294	! !
295	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
296
297
298
299
300	return
301	end subroutine PHY_vdfTKE_RUN

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

Download in other formats:

Original Format