Context Navigation

thermcell_height.F90 @ 2130

Last change on this file since 2130 was 2127, checked in by aboissinot, 7 years ago
Fix in convadj.F New version of the thermal plume model (new entrainment and detrainment formulae, alimentation is removed)
File size: 8.4 KB

Line
1	!
2	!
3	!
4	SUBROUTINE thermcell_height(ngrid,nlay,lmin,linter,lmix,zw2, &
5	zlev,lmax,zmax,zmix,wmax,f_star)
6
7
8	!===============================================================================
9	! Purpose: calcul des caracteristiques du thermique: zmax,wmax,zmix
10	!===============================================================================
11
12	USE thermcell_mod
13	USE print_control_mod, ONLY: prt_level
14
15	IMPLICIT NONE
16
17
18	!===============================================================================
19	! Declaration
20	!===============================================================================
21
22	! Inputs:
23	! -------
24
25	INTEGER ngrid, nlay
26	INTEGER lmin(ngrid)
27
28	REAL zlev(ngrid,nlay+1)
29	REAL f_star(ngrid,nlay+1)
30
31	! Outputs:
32	! --------
33
34	INTEGER lmix(ngrid)
35	INTEGER lmax(ngrid)
36
37	REAL linter(ngrid)
38	REAL zmax(ngrid) ! maximal vertical speed
39	REAL zmix(ngrid) ! altitude of maximal vertical speed
40	REAL wmax(ngrid) ! maximal vertical speed
41	REAL zw2(ngrid,nlay+1) ! square vertical speed (becomes its squere root)
42
43	! Local:
44	! ------
45
46	INTEGER ig, l
47
48	REAL num(ngrid)
49	REAL denom(ngrid)
50	REAL zlevinter(ngrid)
51	REAL zlev2(ngrid,nlay+1)
52
53	!===============================================================================
54	! Initialization
55	!===============================================================================
56
57	DO ig=1,ngrid
58	lmax(ig) = lmin(ig)
59	lmix(ig) = lmin(ig)
60	ENDDO
61
62	DO ig=1,ngrid
63	wmax(ig) = 0.
64	ENDDO
65
66	DO ig=1,ngrid
67	zmax(ig) = 0.
68	zlevinter(ig) = zlev(ig,1)
69	ENDDO
70
71	DO ig=1,ngrid
72	zlev2(ig,:) = zlev(ig,:) - zlev(ig,lmin(ig))
73	ENDDO
74
75	!===============================================================================
76	! Calcul de la hauteur max du thermique
77	!===============================================================================
78
79	!-------------------------------------------------------------------------------
80	! Calcul de lmax
81	!-------------------------------------------------------------------------------
82
83	DO ig=1,ngrid
84	DO l=nlay,lmin(ig)+1,-1
85	IF (zw2(ig,l).le.0..or.f_star(ig,l).le.0.) THEN
86	lmax(ig) = l - 1
87	ENDIF
88	ENDDO
89	ENDDO
90
91	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
92	! On traite le cas particulier qu'il faudrait eviter ou le thermique
93	! atteind le haut du modele ...
94	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
95	DO ig=1,ngrid
96	IF (zw2(ig,nlay).gt.0.) THEN
97	print *, 'WARNING: a thermal plume reaches the highest layer!'
98	print *, 'ig,l', ig, nlay
99	print *, 'zw2', zw2(ig,nlay)
100	lmax(ig) = nlay
101	ENDIF
102	ENDDO
103
104	!-------------------------------------------------------------------------------
105	! Calcul de zmax continu via le linter
106	!-------------------------------------------------------------------------------
107
108	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
109	! AB : lmin=lmax means the plume is not active and then zw2=0 at each level.
110	! Otherwise we have lmax>lmin.
111	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
112
113	DO ig=1,ngrid
114	l = lmax(ig)
115
116	IF (l==lmin(ig)) THEN
117	linter(ig) = l
118	ELSE
119	linter(ig) = l - zw2(ig,l) / (zw2(ig,l+1) - zw2(ig,l))
120	ENDIF
121	ENDDO
122
123	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
124	! AB : zmax is now equal to zlevinter-zlev(lmin) because lmin can be different
125	! from 1.
126	! We have to substract this level because zmax must be the plume height (to
127	! derive the good closure), not the plume highest altitude.
128	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
129
130	DO ig=1,ngrid
131	zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig) &
132	& + zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1) &
133	& - zlev(ig,lmax(ig)))
134	zmax(ig) = max(zmax(ig), zlevinter(ig) - zlev(ig,lmin(ig)))
135	ENDDO
136
137	!===============================================================================
138	! Calcul de wmax et du niveau d'inversion.
139	!===============================================================================
140
141	!-------------------------------------------------------------------------------
142	! Calcul de lmix et wmax
143	!-------------------------------------------------------------------------------
144
145	DO l=1,nlay
146	DO ig=1,ngrid
147	IF(l.le.lmax(ig) .and. l.gt.lmin(ig)) THEN
148	IF (zw2(ig,l).lt.0.) THEN
149	print *, 'ERROR: zw2 has negative value(s)!'
150	print *, 'ig,l', ig, l
151	print *, 'zw2', zw2(ig,l)
152	CALL abort
153	ENDIF
154
155	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
156	! AB : WARNING zw2 becomes its square root!
157	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
158	zw2(ig,l) = sqrt(zw2(ig,l))
159
160	IF (zw2(ig,l).gt.wmax(ig)) THEN
161	wmax(ig) = zw2(ig,l)
162	lmix(ig) = l
163	ENDIF
164	ELSE
165	zw2(ig,l) = 0.
166	ENDIF
167	ENDDO
168	ENDDO
169
170	!-------------------------------------------------------------------------------
171	! Calcul de zmix continu (profil parabolique des vitesses)
172	!-------------------------------------------------------------------------------
173
174	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
175	! AB : We substract zlev(lmin) in zmix formula because we have to
176	! compare it with zmax which is zlev(linter)-zlev(lmin).
177	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
178
179	DO ig=1,ngrid
180	IF (lmix(ig).gt.lmin(ig)) THEN
181	IF (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig))) &
182	& *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1))) &
183	& -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1)) &
184	& *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10) &
185	& THEN
186	zmix(ig) = ((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig))) &
187	& ((zlev(ig,lmix(ig)))2-(zlev(ig,lmix(ig)+1))*2) &
188	& -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1)) &
189	& ((zlev(ig,lmix(ig)-1))2-(zlev(ig,lmix(ig)))*2)) &
190	& /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig))) &
191	& *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1))) &
192	& -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1)) &
193	& *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig)))))) &
194	& - zlev(ig,lmin(ig))
195	ELSE
196	zmix(ig) = zlev(ig,lmix(ig)) - zlev(ig,lmin(ig))
197	print *, 'WARNING: problematic zmix value!'
198	ENDIF
199	ELSE
200	zmix(ig) = 0.
201	ENDIF
202	ENDDO
203
204	!===============================================================================
205	! Verification zmax > zmix
206	!===============================================================================
207
208	DO ig=1,ngrid
209	IF ((zmax(ig)-zmix(ig)).lt.0.) THEN
210	zmix(ig) = 0.9 * zmax(ig)
211	print *, 'WARNING: we have zmix > zmax!'
212	print *, 'zmax,zmix', zmax(ig), zmix(ig)
213	ENDIF
214	ENDDO
215
216	!-------------------------------------------------------------------------------
217	! Calcul du nouveau lmix correspondant
218	!-------------------------------------------------------------------------------
219
220	DO ig=1,ngrid
221	DO l=1,nlay
222	IF (zmix(ig).ge.zlev2(ig,l).and.zmix(ig).lt.zlev2(ig,l+1)) THEN
223	lmix(ig) = l
224	ENDIF
225	ENDDO
226	ENDDO
227
228
229	RETURN
230	END

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

Download in other formats:

Original Format