Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

PHY_Atm_AT_RUN.f90 @ 3956

Last change on this file since 3956 was 2089, checked in by Laurent Fairhead, 10 years ago

Inclusion de la physique de MAR

Integration of MAR physics

File size: 13.8 KB

Line
1	subroutine PHY_Atm_AT_RUN(FlagSV_KzT,FlagCM)
2
3	!------------------------------------------------------------------------------+
4	! Sat 22-Jun-2013 MAR \|
5	! MAR PHY_Atm_AT_RUN \|
6	! subroutine PHY_Atm_AT_RUN performs Turbulent Vertical Diffusion Scheme \|
7	! \|
8	! version 3.p.4.1 created by H. Gallee, Wed 13-Mar-2013 \|
9	! Last Modification by H. Gallee, Sat 22-Jun-2013 \|
10	! \|
11	!------------------------------------------------------------------------------+
12
13	use Mod_Real
14	use Mod_PHY____dat
15	use Mod_PHY____grd
16	use Mod_PHY_AT_ctr
17	use Mod_PHY_AT_grd
18	use Mod_PHY_AT_kkl
19	use Mod_SISVAT_gpt
20
21
22	IMPLICIT NONE
23
24
25	logical :: FlagSV_KzT
26	logical :: FlagCM
27
28	integer :: i, j, k, ikl
29
30
31
32
33	! Assignation of Mod_PHY_AT_kkl
34	! ================================
35
36	DO ikl = 1,kcolp
37
38	i = ii__AP(ikl)
39	j = jj__AP(ikl)
40
41	LMO_AT(ikl) = LMO_SV_gpt(ikl)
42
43	ENDDO
44
45
46
47
48	! Turbulent Kinetic Energy
49	! ========================
50
51	IF (AT_TKE) THEN
52
53	! **************
54	CALL PHY_vdfTKE_RUN
55	CALL PHY_genTKE_RUN
56	! **************
57
58	ELSE
59	DO ikl = 1,kcolp
60	DO k = 1,mzp
61	TKE_AT(ikl,k)= eps6
62	eps_AT(ikl,k)= epsn
63	ENDDO
64	ENDDO
65	END IF
66
67
68
69
70	! Turbulent Vertical Diffusion
71	! ============================
72
73	NewAAT = .TRUE.
74	SchmAT = 'u'
75	! ******
76	CALL Atm_AT
77	! ******
78
79
80	NewAAT = .FALSE.
81	SchmAT = 'v'
82	! ******
83	CALL Atm_AT
84	! ******
85
86
87	IF (.NOT.FlagSV_KzT) THEN
88
89	SchmAT = 'T'
90	! ******
91	CALL Atm_AT
92	! ******
93	ELSE
94
95	DO k =1,mzp
96	DO ikl=1,kcolp
97	dpktAT(ikl,k) = 0.0000
98	ENDDO
99	ENDDO
100
101	END IF
102
103
104	SchmAT = 'q'
105	! ******
106	CALL Atm_AT
107	! ******
108
109
110	IF (FlagCM) THEN
111
112	SchmAT = 'w'
113	! ******
114	CALL Atm_AT
115	! ******
116
117
118	SchmAT = 'i'
119	! ******
120	CALL Atm_AT
121	! ******
122
123
124	SchmAT = 'C'
125	! ******
126	CALL Atm_AT
127	! ******
128
129
130	SchmAT = 's'
131	! ******
132	CALL Atm_AT
133	! ******
134
135
136	SchmAT = 'r'
137	! ******
138	CALL Atm_AT
139	! ******
140
141
142	END IF
143
144
145	end subroutine PHY_Atm_AT_RUN
146
147
148
149
150	subroutine Atm_AT
151
152	!------------------------------------------------------------------------------+
153	! Sat 22-Jun-2013 MAR \|
154	! MAR Atm_AT \|
155	! subroutine Atm_AT solves Turbulent Vertical Diffusion Scheme \|
156	! \|
157	! version 3.p.4.1 created by H. Gallee, Wed 13-Mar-2013 \|
158	! Last Modification by H. Gallee, Sat 22-Jun-2013 \|
159	! \|
160	!------------------------------------------------------------------------------+
161
162	use Mod_Real
163	use Mod_PHY____dat
164	use Mod_PHY____grd
165	use Mod_PHY_AT_grd
166	use Mod_PHY_AT_kkl
167	use Mod_PHY_CM_kkl
168	use Mod_PHY_DY_kkl
169	use Mod_SISVAT_gpt
170
171
172
173	IMPLICIT NONE
174
175
176
177	integer :: i, j, k, ikl
178	real(kind=real8) :: psa_sq
179
180
181
182
183	! ++++++++++++++++
184	DO ikl = 1,kcolp
185	! ++++++++++++++++
186
187	i = ii__AP(ikl)
188	j = jj__AP(ikl)
189
190
191
192	! Tri-Diagonal Matrix: Variable
193	! =============================
194
195	! Wind Speed
196	! ----------
197
198	IF (schmAT.EQ.'a') THEN
199	DO k=1,mzp
200	Kz__AT(ikl,k) = 1.0
201	var_AT(ikl,k) = ua__DY(ikl,k)
202	ENDDO
203	k= mzpp
204	var_AT(ikl,k) = 0.
205
206	ELSE IF (schmAT.EQ.'b') THEN
207	DO k=1,mzp
208	Kz__AT(ikl,k) = 1.0
209	var_AT(ikl,k) = va__DY(ikl,k)
210	ENDDO
211	k= mzpp
212	var_AT(ikl,k) = 0.
213
214	ELSE IF (schmAT.EQ.'u') THEN
215	DO k=1,mzp
216	Kz__AT(ikl,k) = Kzm_AT(ikl,k)
217	var_AT(ikl,k) = ua__DY(ikl,k)
218	ENDDO
219	k= mzpp
220	var_AT(ikl,k) = 0.
221
222	ELSE IF (schmAT.EQ.'v') THEN
223	DO k=1,mzp
224	Kz__AT(ikl,k) = Kzm_AT(ikl,k)
225	var_AT(ikl,k) = va__DY(ikl,k)
226	ENDDO
227	k= mzpp
228	var_AT(ikl,k) = 0.
229
230
231	! Temperature
232	! -----------
233
234	ELSE IF (schmAT.EQ.'T') THEN
235	DO k=1,mzp
236	Kz__AT(ikl,k) = Kzh_AT(ikl,k)
237	var_AT(ikl,k) = pkt_DY(ikl,k)
238	ENDDO
239	k= mzpp
240	var_AT(ikl,k) = pkt_DY(ikl,k)
241
242
243	! Water Species
244	! -------------
245
246	ELSE IF (schmAT.EQ.'q') THEN
247	DO k=1,mzp
248	Kz__AT(ikl,k) = Kzh_AT(ikl,k)
249	var_AT(ikl,k) = qv__DY(ikl,k)
250	ENDDO
251	k= mzpp
252	var_AT(ikl,k) = qv__DY(ikl,mzp) &
253	& -(Z___DY(ikl,mzp)-Z___DY(ikl,mzpp)) &
254	& * uqs_SV_gpt(ikl) /Kzh_AT(ikl,mzp )
255	var_AT(ikl,k) = max(epsq,var_AT(ikl,k ))
256	var_AT(ikl,k) = min( var_AT(ikl,k ),qvsiCM(ikl,k))
257	qv__DY(ikl,mzpp) = var_AT(ikl,k )
258
259	ELSE IF (schmAT.EQ.'w') THEN
260	DO k=1,mzp
261	Kz__AT(ikl,k) = Kzh_AT(ikl,k)
262	var_AT(ikl,k) = qw__CM(ikl,k)
263	ENDDO
264	k= mzpp
265	var_AT(ikl,k) = qw__CM(ikl,mzp)
266
267	ELSE IF (schmAT.EQ.'i') THEN
268	DO k=1,mzp
269	Kz__AT(ikl,k) = Kzh_AT(ikl,k)
270	var_AT(ikl,k) = qi__CM(ikl,k)
271	ENDDO
272	k= mzpp
273	var_AT(ikl,k) = qi__CM(ikl,mzp)
274
275	ELSE IF (schmAT.EQ.'C') THEN
276	DO k=1,mzp
277	Kz__AT(ikl,k) = Kzh_AT(ikl,k)
278	var_AT(ikl,k) = CCNiCM(ikl,k)
279	ENDDO
280	k= mzpp
281	var_AT(ikl,k) = CCNiCM(ikl,mzp)
282
283	ELSE IF (schmAT.EQ.'s') THEN
284	DO k=1,mzp
285	Kz__AT(ikl,k) = Kzh_AT(ikl,k)
286	var_AT(ikl,k) = qs__CM(ikl,k)
287	ENDDO
288	k= mzpp
289	var_AT(ikl,k) = qs__CM(ikl,k)
290
291	ELSE IF (schmAT.EQ.'r') THEN
292	DO k=1,mzp
293	Kz__AT(ikl,k) = Kzh_AT(ikl,k)
294	var_AT(ikl,k) = qr__CM(ikl,k)
295	ENDDO
296	k= mzpp
297	var_AT(ikl,k) = qr__CM(ikl,mzp)
298
299	ELSE
300	write(6,*)
301	write(6,*) ' schmAT = ',schmAT
302	write(6,*) ' Inadequate..., NOR for Ekman Spiral, NOR for momentum, NOR for scalars'
303	write(6,*) ' '
304
305	! **************
306	call MAR________BUG
307	! **************
308	STOP
309	END IF
310
311
312
313	! Tri-Diagonal Matrix Coefficients
314	! ================================
315
316	IF (NewAAT) THEN
317	psa_sq = psa_DY(ikl) *psa_DY(ikl)
318	Ac__AT(ikl,mzp) = Ac0_AT(mzp) &
319	& * roa_DY(ikl,mzp) *roamDY(ikl,mzp) &
320	& / psa_sq
321	DO k=mzp-1,1,-1
322	Ac__AT(ikl,k ) = Ac0_AT(k ) &
323	& * roa_DY(ikl,k ) *roamDY(ikl,k ) &
324	& / psa_sq
325	Cc__AT(ikl,k+1) = Ac__AT(ikl,k) &
326	& * roa_DY(ikl,k+1) /roa_DY(ikl,k )
327	ENDDO
328	ENDIF
329
330	A___AT(ikl,mzp) = Ac__AT(ikl,mzp) * Kz__AT(ikl,mzp)
331	DO k=mzp-1,1,-1
332	A___AT(ikl,k ) = Ac__AT(ikl,k ) * Kz__AT(ikl,k )
333	C___AT(ikl,k+1) = Cc__AT(ikl,k+1) * Kz__AT(ikl,k )
334	ENDDO
335	DO k= 2,mzp
336	B___AT(ikl,k ) = 1.0000 - A___AT(ikl,k ) - C___AT(ikl,k )
337	ENDDO
338	B___AT(ikl,1 ) = 1.0000 - A___AT(ikl,1 )
339	C___AT(ikl,1 ) = 0.0000
340
341	D___AT(ikl,1 ) = &
342	& A___AT(ikl,1 ) * a_b_AT *(var_AT(ikl,1 ) - var_AT(ikl,2 ))
343
344	k= mzp
345	D___AT(ikl,k ) = &
346	& A___AT(ikl,k ) * (a_b_AT * var_AT(ikl,k ) - var_AT(ikl,k+1) &
347	& / betaAT) &
348	& - C___AT(ikl,k ) * a_b_AT *(var_AT(ikl,k-1) - var_AT(ikl,k )) &
349	& + var_AT(ikl,k )
350	A___AT(ikl,k ) = 0.0000
351	DO k= 2,mzp-1
352	D___AT(ikl,k ) = &
353	& A___AT(ikl,k ) * a_b_AT *(var_AT(ikl,k ) - var_AT(ikl,k+1)) &
354	& - C___AT(ikl,k ) * a_b_AT *(var_AT(ikl,k-1) - var_AT(ikl,k )) &
355	& + var_AT(ikl,k )
356	ENDDO
357	D___AT(ikl,1 ) = &
358	& A___AT(ikl,1 ) * a_b_AT *(var_AT(ikl,1 ) - var_AT(ikl,2 )) &
359	& + var_AT(ikl,1 )
360
361
362
363	! Tri-Diagonal Matrix Resolution (Gaussian Elimination / Thomas Algorithm)
364	! ==============================
365
366	! Forward Sweep
367	! --------------
368
369	P___AT(1) = B___AT(ikl,1)
370	Q___AT(1) = -A___AT(ikl,1) / P___AT(1)
371	DO k=2,mzp
372	P___AT(k) = C___AT(ikl,k) * Q___AT(k-1) + B___AT(ikl,k)
373	Q___AT(k) = -A___AT(ikl,k) / P___AT(k)
374	ENDDO
375	X___AT(1) = D___AT(ikl,1) / P___AT(1)
376	DO k=2,mzp
377	X___AT(k) = (D___AT(ikl,k) - C___AT(ikl,k) * X___AT(k-1)) / P___AT(k)
378	END DO
379
380
381	! Backward Sweep
382	! --------------
383
384	DO k= mzp-1,1,-1
385	X___AT(k) = Q___AT(k) * X___AT(k+1) + X___AT(k)
386	END DO
387
388
389
390	! Elimination of too small values
391	! -------------------------------
392
393	IF (SchmAT.EQ.'w'.OR.SchmAT.EQ.'i'.OR.SchmAT.EQ.'C'.OR. &
394	& SchmAT.EQ.'s'.OR.SchmAT.EQ.'r') THEN
395	DO k= mzp ,1,-1
396	IF (X___AT(k).NE.0.0000 .AND. X___AT(k).LT.1.e-18) &
397	& X___AT(k) = 0.0000
398	END DO
399	END IF
400
401
402
403
404	! Tendency
405	! ========
406
407	! Wind Speed
408	! ----------
409
410	IF (schmAT.EQ.'a') THEN
411	DO k=1,mzp
412	dua_AT(ikl,k) =(X___AT(k) - ua__DY(ikl,k)) / dt__AT
413	ENDDO
414	ELSE IF (schmAT.EQ.'b') THEN
415	DO k=1,mzp
416	dva_AT(ikl,k) =(X___AT(k) - va__DY(ikl,k)) / dt__AT
417	ENDDO
418	ELSE IF (schmAT.EQ.'u') THEN
419	DO k=1,mzp
420	dua_AT(ikl,k) =(X___AT(k) - ua__DY(ikl,k)) / dt__AT
421	ENDDO
422	ELSE IF (schmAT.EQ.'v') THEN
423	DO k=1,mzp
424	dva_AT(ikl,k) =(X___AT(k) - va__DY(ikl,k)) / dt__AT
425	ENDDO
426
427
428	! Temperature
429	! -----------
430
431	ELSE IF (schmAT.EQ.'T') THEN
432	DO k=1,mzp
433	dpktAT(ikl,k) =(X___AT(k) - pkt_DY(ikl,k)) / dt__AT
434	ENDDO
435
436
437	! Water Species
438	! -------------
439
440	ELSE IF (schmAT.EQ.'q') THEN
441	DO k=1,mzp
442	dqv_AT(ikl,k) =(X___AT(k) - qv__DY(ikl,k)) / dt__AT
443	ENDDO
444
445	ELSE IF (schmAT.EQ.'w') THEN
446	DO k=1,mzp
447	dqw_AT(ikl,k) =(X___AT(k) - qw__CM(ikl,k)) / dt__AT
448	ENDDO
449
450	ELSE IF (schmAT.EQ.'i') THEN
451	DO k=1,mzp
452	dqi_AT(ikl,k) =(X___AT(k) - qi__CM(ikl,k)) / dt__AT
453	ENDDO
454
455	ELSE IF (schmAT.EQ.'C') THEN
456	DO k=1,mzp
457	dCi_AT(ikl,k) =(X___AT(k) - CCNiCM(ikl,k)) / dt__AT
458	ENDDO
459
460	ELSE IF (schmAT.EQ.'s') THEN
461	DO k=1,mzp
462	dqs_AT(ikl,k) =(X___AT(k) - qs__CM(ikl,k)) / dt__AT
463	ENDDO
464
465	ELSE IF (schmAT.EQ.'r') THEN
466	DO k=1,mzp
467	dqr_AT(ikl,k) =(X___AT(k) - qr__CM(ikl,k)) / dt__AT
468	ENDDO
469
470	ELSE
471	write(6,*)
472	write(6,*) ' schmAT = ',schmAT
473	write(6,*) ' Inadequate..., NOR for Ekman Spiral, NOR for momentum, NOR for scalars'
474	write(6,*) ' '
475
476	! **************
477	call MAR________BUG
478	! **************
479	STOP
480	END IF
481
482
483	! ++++++
484	END DO
485	! ++++++
486
487
488
489
490	end subroutine Atm_AT

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

Download in other formats:

Original Format