Context Navigation

thermcell_dv2.F90 @ 3436

Last change on this file since 3436 was 2231, checked in by aboissinot, 5 years ago
Fix a bug in thermcell_dv2: now the plume height (zmin-zmax) is used for computations instead of maximal altitude (zmax).
File size: 6.6 KB

Line
1	!
2	!
3	!
4	SUBROUTINE thermcell_dv2(ngrid,nlay,ptimestep,fm,entr,detr,masse,fraca, &
5	zmax,zmin,u,v,du,dv,ua,va)
6
7
8	!===============================================================================
9	!
10	! Calcul du transport vertical dans la couche limite en presence
11	! de "thermiques" explicitement representes
12	! Calcul du dq/dt une fois qu'on connait les ascendances
13	!
14	! Vectorisation, FH : 2010/03/08
15	!
16	!===============================================================================
17
18	IMPLICIT none
19
20
21	!===============================================================================
22	! Declaration
23	!===============================================================================
24
25	! Inputs:
26	! -------
27
28	INTEGER, INTENT(in) :: ngrid
29	INTEGER, INTENT(in) :: nlay
30
31	REAL, INTENT(in) :: ptimestep
32	REAL, INTENT(in) :: masse(ngrid,nlay)
33	REAL, INTENT(in) :: fm(ngrid,nlay+1)
34	REAL, INTENT(in) :: entr(ngrid,nlay)
35	REAL, INTENT(in) :: detr(ngrid,nlay)
36	REAL, INTENT(in) :: fraca(ngrid,nlay+1)
37	REAL, INTENT(in) :: zmax(ngrid)
38	REAL, INTENT(in) :: zmin(ngrid)
39	REAL, INTENT(in) :: u(ngrid,nlay)
40	REAL, INTENT(in) :: v(ngrid,nlay)
41
42	! Outputs:
43	! --------
44
45	REAL, INTENT(out) :: ua(ngrid,nlay) ! u in the plume
46	REAL, INTENT(out) :: va(ngrid,nlay) ! v in the plume
47	REAL, INTENT(out) :: du(ngrid,nlay) ! large scale u variation
48	REAL, INTENT(out) :: dv(ngrid,nlay) ! large scale v variation
49
50	! Local:
51	! ------
52
53	INTEGER ig, l
54	INTEGER iter
55
56	REAL qa(ngrid,nlay)
57	REAL zf
58	REAL zf2
59	REAL wvd(ngrid,nlay+1)
60	REAL wud(ngrid,nlay+1)
61	REAL gamma0(ngrid,nlay+1)
62	REAL gamma(ngrid,nlay+1)
63	REAL ue(ngrid,nlay)
64	REAL ve(ngrid,nlay)
65	REAL dua(ngrid,nlay)
66	REAL dva(ngrid,nlay)
67	REAL plume_height(ngrid)
68
69	LOGICAL ltherm(ngrid,nlay)
70
71	!===============================================================================
72	! Initialization
73	!===============================================================================
74
75	DO ig=1,ngrid
76	plume_height(ig) = zmax(ig) - zmin(ig)
77	ENDDO
78
79	DO ig=1,ngrid
80	ua(ig,1)=u(ig,1)
81	va(ig,1)=v(ig,1)
82	ue(ig,1)=u(ig,1)
83	ve(ig,1)=v(ig,1)
84	ENDDO
85
86	gamma0(1:ngrid,1)=0.
87
88	DO l=2,nlay
89	DO ig=1,ngrid
90	ltherm(ig,l) = (fm(ig,l+1) + detr(ig,l)) * ptimestep > 1.e-5 * masse(ig,l)
91	! IF (ltherm(ig,l).and.(zmax(ig) > 0.)) THEN
92	IF (ltherm(ig,l).and.(plume_height(ig) > 0.)) THEN
93	! gamma0(ig,l) = masse(ig,l) * 0.5 / zmax(ig) &
94	gamma0(ig,l) = masse(ig,l) * 0.5 / plume_height(ig) &
95	& * sqrt(0.5 * (fraca(ig,l+1) + fraca(ig,l)))
96	ELSE
97	gamma0(ig,l) = 0.
98	ENDIF
99	ENDDO
100	ENDDO
101
102	gamma(:,:) = 0.
103
104	!===============================================================================
105	!
106	!===============================================================================
107
108	DO l=2,nlay
109
110	DO ig=1,ngrid
111	IF (ltherm(ig,l)) THEN
112	dua(ig,l) = ua(ig,l-1) - u(ig,l-1)
113	dva(ig,l) = va(ig,l-1) - v(ig,l-1)
114	ELSE
115	ua(ig,l) = u(ig,l)
116	va(ig,l) = v(ig,l)
117	ue(ig,l) = u(ig,l)
118	ve(ig,l) = v(ig,l)
119	ENDIF
120	ENDDO
121
122	DO iter=1,5
123	DO ig=1,ngrid
124	! FH: Calcul prenant en compte la fraction reelle
125	zf = 0.5 * (fraca(ig,l) + fraca(ig,l+1))
126	zf2 = 1./(1.-zf)
127	! FH: Calcul avec fraction infiniement petite
128	! zf = 0.
129	! zf2 = 1.
130
131	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
132	! FH: La première fois on multiplie le coefficient de freinage par le module du
133	! vent dans la couche en dessous. Mais pourquoi donc ?
134	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
135	IF (ltherm(ig,l)) THEN
136	! FH: On choisit une relaxation lineaire.
137	! gamma(ig,l) = gamma0(ig,l)
138	! FH: On choisit une relaxation quadratique.
139	gamma(ig,l) = gamma0(ig,l) * sqrt(dua(ig,l)2 + dva(ig,l)2)
140
141	ua(ig,l) = (fm(ig,l) * ua(ig,l-1) &
142	& + (zf2 * entr(ig,l) + gamma(ig,l)) * u(ig,l)) &
143	& / (fm(ig,l+1) + detr(ig,l) + entr(ig,l) * zf * zf2 &
144	& + gamma(ig,l))
145	va(ig,l) = (fm(ig,l) * va(ig,l-1) &
146	& + (zf2 * entr(ig,l) + gamma(ig,l)) * v(ig,l)) &
147	& / (fm(ig,l+1) + detr(ig,l) + entr(ig,l) * zf * zf2 &
148	& + gamma(ig,l))
149	dua(ig,l) = ua(ig,l) - u(ig,l)
150	dva(ig,l) = va(ig,l) - v(ig,l)
151	ue(ig,l) = (u(ig,l) - zf * ua(ig,l)) * zf2
152	ve(ig,l) = (v(ig,l) - zf * va(ig,l)) * zf2
153	ENDIF
154	ENDDO
155	ENDDO
156	ENDDO
157
158	!-------------------------------------------------------------------------------
159	! Calcul du flux vertical de moment dans l'environnement
160	!-------------------------------------------------------------------------------
161
162	DO l=2,nlay
163	DO ig=1,ngrid
164	wud(ig,l) = fm(ig,l) * ue(ig,l)
165	wvd(ig,l) = fm(ig,l) * ve(ig,l)
166	ENDDO
167	ENDDO
168
169	DO ig=1,ngrid
170	wud(ig,1) = 0.
171	wvd(ig,1) = 0.
172	wud(ig,nlay+1) = 0.
173	wvd(ig,nlay+1) = 0.
174	ENDDO
175
176	!===============================================================================
177	! Tendencies
178	!===============================================================================
179
180	DO l=1,nlay
181	DO ig=1,ngrid
182	du(ig,l) = ((detr(ig,l) + gamma(ig,l)) * ua(ig,l) &
183	& - (entr(ig,l) + gamma(ig,l)) * ue(ig,l) &
184	& - wud(ig,l) + wud(ig,l+1)) &
185	& / masse(ig,l)
186	dv(ig,l) = ((detr(ig,l) + gamma(ig,l)) * va(ig,l) &
187	& - (entr(ig,l) + gamma(ig,l)) * ve(ig,l) &
188	& - wvd(ig,l) + wvd(ig,l+1)) &
189	& / masse(ig,l)
190	ENDDO
191	ENDDO
192
193
194	RETURN
195	END

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

Download in other formats:

Original Format