Context Navigation

moistadj.F90 @ 937

Last change on this file since 937 was 875, checked in by jleconte, 12 years ago

11/02/2013 == JL

Updated moist convection scheme to handle situations with a large water vapor content
Added a keyword to enable ocean runoff in callphys.def (activerunoff)

File size: 11.4 KB

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

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

Download in other formats:

Original Format