Context Navigation

moistadj.F90 @ 723

Last change on this file since 723 was 650, checked in by jleconte, 13 years ago

Corrected the temperature used to differentiate sublimation and evaporation in watersat_grad
Minor name changes in watercommon
Better physical parametrization of the effective radius of liquid and icy water cloud particles in callcorrk

(for radfixed=true)

Added consistency check in inifis
Moved 1d water initialization from physiqu to rcm1d

File size: 12.2 KB

Line
1	subroutine moistadj(t, pq, pplev, pplay, dtmana, dqmana, ptimestep, rneb)
2
3	use watercommon_h, only: T_h2O_ice_liq, RLVTT, RCPD
4
5	implicit none
6
7
8	!=====================================================================
9	!
10	! Purpose
11	! -------
12	! Calculates moist convective adjustment by the method of Manabe.
13	!
14	! Authors
15	! -------
16	! Adapted from the LMDTERRE code by R. Wordsworth (2010)
17	! Original author Z. X. Li (1993)
18	!
19	!=====================================================================
20
21	#include "dimensions.h"
22	#include "dimphys.h"
23	#include "tracer.h"
24	#include "comcstfi.h"
25
26	! Pre-arguments (for universal model)
27	real pq(ngridmx,nlayermx,nqmx) ! tracer (kg/kg)
28	REAL pdq(ngridmx,nlayermx,nqmx)
29
30	real dqmana(ngridmx,nlayermx,nqmx) ! tendency of tracers (kg/kg.s-1)
31	REAL dtmana(ngridmx,nlayermx) ! temperature increment
32
33	! Arguments
34	REAL t(ngridmx,nlayermx) ! temperature (K)
35	REAL q(ngridmx,nlayermx) ! humidite specifique (kg/kg)
36	REAL pplev(ngridmx,nlayermx+1) ! pression a inter-couche (Pa)
37	REAL pplay(ngridmx,nlayermx) ! pression au milieu de couche (Pa)
38
39	REAL d_t(ngridmx,nlayermx) ! temperature increment
40	REAL d_q(ngridmx,nlayermx) ! incrementation pour vapeur d'eau
41	REAL d_ql(ngridmx,nlayermx) ! incrementation pour l'eau liquide
42	REAL rneb(ngridmx,nlayermx) ! cloud fraction
43	REAL ptimestep
44
45	! REAL t_coup
46	! PARAMETER (t_coup=234.0)
47	REAL seuil_vap
48	PARAMETER (seuil_vap=1.0E-10)
49
50	! Local variables
51	INTEGER i, k, iq
52	INTEGER k1, k1p, k2, k2p
53	LOGICAL itest(ngridmx)
54	REAL delta_q(ngridmx, nlayermx)
55	REAL cp_new_t(nlayermx)
56	REAL cp_delta_t(nlayermx)
57	REAL new_qb(nlayermx)
58	REAL v_cptj(nlayermx), v_cptjk1, v_ssig
59	REAL v_cptt(ngridmx,nlayermx), v_p, v_t
60	REAL v_qs(ngridmx,nlayermx), v_qsd(ngridmx,nlayermx)
61	REAL zq1(ngridmx), zq2(ngridmx)
62	REAL gamcpdz(ngridmx,2:nlayermx)
63	REAL zdp, zdpm
64
65	REAL zsat ! super-saturation
66	REAL zflo ! flotabilite
67
68	REAL local_q(ngridmx,nlayermx),local_t(ngridmx,nlayermx)
69
70	REAL zdelta, zcor, zcvm5
71
72	REAL dEtot, dqtot, masse ! conservation diagnostics
73	real dL1tot, dL2tot
74
75	! Indices of water vapour and water ice tracers
76	INTEGER,SAVE :: i_h2o=0 ! water vapour
77	INTEGER,SAVE :: i_ice=0 ! water ice
78
79	LOGICAL firstcall
80	SAVE firstcall
81
82	DATA firstcall /.TRUE./
83
84	IF (firstcall) THEN
85
86	i_h2o=igcm_h2o_vap
87	i_ice=igcm_h2o_ice
88
89	write(,) "rain: i_ice=",i_ice
90	write(,) " i_h2o=",i_h2o
91
92	firstcall = .FALSE.
93	ENDIF
94
95	! GCM -----> subroutine variables
96	DO k = 1, nlayermx
97	DO i = 1, ngridmx
98
99	q(i,k) = pq(i,k,i_h2o)
100
101	if(q(i,k).lt.0.)then
102	q(i,k)=0.0
103	endif
104	DO iq = 1, nqmx
105	dqmana(i,k,iq)=0.0
106	ENDDO
107	ENDDO
108	ENDDO
109
110	DO k = 1, nlayermx
111	DO i = 1, ngridmx
112	local_q(i,k) = q(i,k)
113	local_t(i,k) = t(i,k)
114	rneb(i,k) = 0.0
115	d_ql(i,k) = 0.0
116	d_t(i,k) = 0.0
117	d_q(i,k) = 0.0
118	ENDDO
119	new_qb(k)=0.0
120	ENDDO
121
122	! Calculate v_cptt
123	DO k = 1, nlayermx
124	DO i = 1, ngridmx
125	v_cptt(i,k) = RCPD * local_t(i,k)
126	v_t = local_t(i,k)
127	v_p = pplay(i,k)
128
129	call watersat(v_t,v_p,v_qs(i,k))
130	call watersat_grad(v_t,v_qs(i,k),v_qsd(i,k))
131	ENDDO
132	ENDDO
133
134	! TEST: RH DIAGNOSTIC
135	! DO k = 1, nlayermx
136	! DO i = 1, ngridmx
137	! v_t = local_t(i,k)
138	! IF (v_t.LT.T_h2O_ice_liq) THEN
139	! print*,'RHs=',q(i,k) / v_qs(i,k)
140	! ELSE
141	! print*,'RHl=',q(i,k) / v_qs(i,k)
142	! ENDIF
143	! ENDDO
144	! ENDDO
145
146	! Calculate Gamma * Cp * dz: (gamma is the critical gradient)
147	DO k = 2, nlayermx
148	DO i = 1, ngridmx
149	zdp = pplev(i,k)-pplev(i,k+1)
150	zdpm = pplev(i,k-1)-pplev(i,k)
151	! gamcpdz(i,k) = ( ( RD/RCPD /(zdpm+zdp) *
152	gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * &
153	(v_cptt(i,k-1)zdpm + v_cptt(i,k)zdp) &
154	+RLVTT /(zdpm+zdp) * &
155	(v_qs(i,k-1)zdpm + v_qs(i,k)zdp) &
156	)* (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) &
157	/ (1.0+(v_qsd(i,k-1)*zdpm+ &
158	v_qsd(i,k)*zdp)/(zdpm+zdp) )
159	ENDDO
160	ENDDO
161
162	!------------------------------------ modification of unstable profile
163	DO 9999 i = 1, ngridmx
164	itest(i) = .FALSE.
165
166	! print*,'we in the loop'
167	! stop
168
169	k1 = 0
170	k2 = 1
171
172	810 CONTINUE ! look for k1, the base of the column
173	k2 = k2 + 1
174	IF (k2 .GT. nlayermx) GOTO 9999
175	zflo = v_cptt(i,k2-1) - v_cptt(i,k2) - gamcpdz(i,k2)
176	zsat=(local_q(i,k2-1)-v_qs(i,k2-1))*(pplev(i,k2-1)-pplev(i,k2)) &
177	+(local_q(i,k2)-v_qs(i,k2))*(pplev(i,k2)-pplev(i,k2+1))
178
179	IF ( zflo.LE.0.0 .OR. zsat.LE.0.0 ) GOTO 810
180	k1 = k2 - 1
181	itest(i) = .TRUE.
182
183	820 CONTINUE !! look for k2, the top of the column
184	IF (k2 .EQ. nlayermx) GOTO 821
185	k2p = k2 + 1
186	zsat=zsat+(pplev(i,k2p)-pplev(i,k2p+1))*(local_q(i,k2p)-v_qs(i,k2p))
187	zflo = v_cptt(i,k2p-1) - v_cptt(i,k2p) - gamcpdz(i,k2p)
188
189	IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 821
190	k2 = k2p
191	GOTO 820
192	821 CONTINUE
193
194	!------------------------------------------------------ local adjustment
195	830 CONTINUE ! actual adjustment
196	v_cptj(k1) = 0.0
197	zdp = pplev(i,k1)-pplev(i,k1+1)
198	v_cptjk1 = ( (1.0+v_qsd(i,k1))*(v_cptt(i,k1)+v_cptj(k1)) &
199	+ RLVTT(local_q(i,k1)-v_qs(i,k1)) ) zdp
200	v_ssig = zdp * (1.0+v_qsd(i,k1))
201
202	k1p = k1 + 1
203	DO k = k1p, k2
204	zdp = pplev(i,k)-pplev(i,k+1)
205	v_cptj(k) = v_cptj(k-1) + gamcpdz(i,k)
206	v_cptjk1 = v_cptjk1 + zdp &
207	* ( (1.0+v_qsd(i, k))*(v_cptt(i,k)+v_cptj(k)) &
208	+ RLVTT*(local_q(i,k)-v_qs(i,k)) )
209	v_ssig = v_ssig + zdp *(1.0+v_qsd(i,k))
210	ENDDO
211
212
213	! this right here is where the adjustment is done???
214	DO k = k1, k2
215	cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k)
216	cp_delta_t(k) = cp_new_t(k) - v_cptt(i,k)
217	new_qb(k) = v_qs(i,k) + v_qsd(i,k)*cp_delta_t(k)/RLVTT
218	local_q(i,k) = new_qb(k)
219	local_t(i,k) = cp_new_t(k) / RCPD
220
221	! print*,'v_qs in loop=',v_qs
222	! print*,'v_qsd in loop=',v_qsd
223	! print*,'new_qb in loop=',new_qb
224	! print*,'cp_delta_t in loop=',cp_delta_t
225	ENDDO
226
227
228	!--------------------------------------------------- sounding downwards
229	! -- we refine the prognostic variables in
230	! -- the layer about to be adjusted
231
232	DO k = k1, k2
233	v_cptt(i,k) = RCPD * local_t(i,k)
234	v_t = local_t(i,k)
235	v_p = pplay(i,k)
236
237	! IF (v_t.LT.t_coup) THEN
238	! v_qs(i,k) = qsats(v_t) / v_p
239	! v_qsd(i,k) = dqsats(v_t,v_qs(i,k))
240	! ELSE
241	! v_qs(i,k) = qsatl(v_t) / v_p
242	! v_qsd(i,k) = dqsatl(v_t,v_qs(i,k))
243	! ENDIF
244
245	call watersat(v_t,v_p,v_qs(i,k))
246	call watersat_grad(v_t,v_qs(i,k),v_qsd(i,k))
247
248	ENDDO
249	DO k = 2, nlayermx
250	zdpm = pplev(i,k-1) - pplev(i,k)
251	zdp = pplev(i,k) - pplev(i,k+1)
252	! gamcpdz(i,k) = ( ( RD/RCPD /(zdpm+zdp) *
253	gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * &
254	(v_cptt(i,k-1)zdpm+v_cptt(i,k)zdp) &
255	+RLVTT /(zdpm+zdp) * &
256	(v_qs(i,k-1)zdpm+v_qs(i,k)zdp) &
257	)* (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) &
258	/ (1.0+(v_qsd(i,k-1)zdpm+v_qsd(i,k)zdp) &
259	/(zdpm+zdp) )
260	ENDDO
261
262	! Test to see if we've reached the bottom
263
264	IF (k1 .EQ. 1) GOTO 841 ! yes we have!
265	zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
266	zsat=(local_q(i,k1-1)-v_qs(i,k1-1))*(pplev(i,k1-1)-pplev(i,k1)) &
267	+ (local_q(i,k1)-v_qs(i,k1))*(pplev(i,k1)-pplev(i,k1+1))
268	IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 841 ! yes we have!
269
270	840 CONTINUE
271	k1 = k1 - 1
272	IF (k1 .EQ. 1) GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
273	zsat = zsat + (local_q(i,k1-1)-v_qs(i,k1-1)) &
274	*(pplev(i,k1-1)-pplev(i,k1))
275	zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
276	IF (zflo.GT.0.0 .AND. zsat.GT.0.0) THEN
277	GOTO 840
278	ELSE
279	GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
280	ENDIF
281	841 CONTINUE
282
283	GOTO 810 ! look for other layers higher up
284
285	9999 CONTINUE ! loop over all the points
286
287	! print*,'k1=',k1
288	! print*,'k2=',k2
289
290	! print*,'local_t=',local_t
291	! print*,'v_cptt=',v_cptt
292	! print*,'gamcpdz=',gamcpdz
293
294	!-----------------------------------------------------------------------
295	! Determine the cloud fraction (hypothese: la nebulosite a lieu
296	! a l'endroit ou la vapeur d'eau est diminuee par l'ajustement):
297
298	DO k = 1, nlayermx
299	DO i = 1, ngridmx
300	IF (itest(i)) THEN
301	delta_q(i,k) = local_q(i,k) - q(i,k)
302	IF (delta_q(i,k).LT.0.) rneb(i,k) = 1.0
303	ENDIF
304	ENDDO
305	ENDDO
306
307	! Distribuer l'eau condensee en eau liquide nuageuse (hypothese:
308	! l'eau liquide est distribuee aux endroits ou la vapeur d'eau
309	! diminue et d'une maniere proportionnelle a cet diminution):
310
311	DO i = 1, ngridmx
312	IF (itest(i)) THEN
313	zq1(i) = 0.0
314	zq2(i) = 0.0
315	ENDIF
316	ENDDO
317	DO k = 1, nlayermx
318	DO i = 1, ngridmx
319	IF (itest(i)) THEN
320	zdp = pplev(i,k)-pplev(i,k+1)
321	zq1(i) = zq1(i) - delta_q(i,k) * zdp
322	zq2(i) = zq2(i) - MIN(0.0, delta_q(i,k)) * zdp
323	ENDIF
324	ENDDO
325	ENDDO
326	DO k = 1, nlayermx
327	DO i = 1, ngridmx
328	IF (itest(i)) THEN
329	IF (zq2(i).NE.0.0) &
330	d_ql(i,k) = - MIN(0.0,delta_q(i,k))*zq1(i)/zq2(i)
331	ENDIF
332	ENDDO
333	ENDDO
334
335	! print*,'local_q BEFORE=',local_q
336
337	DO k = 1, nlayermx
338	DO i = 1, ngridmx
339	local_q(i, k) = MAX(local_q(i, k), seuil_vap)
340	ENDDO
341	ENDDO
342
343	DO k = 1, nlayermx
344	DO i = 1, ngridmx
345	d_t(i,k) = local_t(i,k) - t(i,k)
346	d_q(i,k) = local_q(i,k) - q(i,k)
347	ENDDO
348	ENDDO
349
350	! now subroutine -----> GCM variables
351	DO k = 1, nlayermx
352	DO i = 1, ngridmx
353
354	dtmana(i,k) = d_t(i,k)/ptimestep
355	dqmana(i,k,i_h2o) = d_q(i,k)/ptimestep
356	dqmana(i,k,i_ice) = d_ql(i,k)/ptimestep
357
358	ENDDO
359	ENDDO
360
361	! print*,'IN MANABE:'
362	! print*,'pplev=',pplev
363	! print*,'t=',t
364	! print*,'d_t=',d_t
365	! print*,'d_q=',d_q
366	! print*,'local_q=',local_q
367	! print*,'q=',q
368	! print*,'v_qs(i,k)=',v_qs
369	! print*,'v_qsd(i,k)=',v_qsd
370	! print*,'cp_delta_t(k)=',cp_delta_t
371
372	! print*,'d_ql=',d_ql
373	! print*,'delta_q=',delta_q
374	! print*,'zq1=',zq1
375	! print*,'zq2=',zq2
376	!! print*,'i_h2o=',i_h2o
377	! print*,'i_ice=',i_ice
378	!
379	! print*,'IN MANABE:'
380	! print*,'d_q=',d_q
381	! print*,'d_ql=',d_ql
382	! print*,'dtmana=',d_t
383	! stop
384	! print*,'gamcpdz at end=',gamcpdz
385	! stop
386
387	! Some conservation diagnostics...
388	! dEtot=0.0
389	! dL1tot=0.0
390	! dL2tot=0.0
391	! dqtot=0.0
392	! masse=0.0
393	! DO k = 1, nlayermx
394	! DO i = 1, ngridmx
395	!
396	! masse = (pplev(i,k) - pplev(i,k+1))/g
397	!
398	! dEtot = dEtot + cppd_t(i,k)masse
399	! dL1tot = dL1tot + RLVTTd_ql(i,k)masse
400	! dL2tot = dL2tot + RLVTTd_q(i,k)masse ! is this line necessary?
401	!
402	! dqtot = dqtot + (d_q(i,k) + d_ql(i,k))*masse
403	!
404	! ENDDO
405	! ENDDO
406
407	! print*,'In manabe energy change=',dEtot
408	! print*,'In manabe condense energy change 1 =',dL1tot
409	! print*,'In manabe condense energy change 2 =',dL2tot
410	! print*,'In manabe water change=',dqtot
411
412	RETURN
413	END

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

Download in other formats:

Original Format