Context Navigation

moistadj.F90 @ 1150

Last change on this file since 1150 was 1016, checked in by jleconte, 11 years ago

07/08/2013 == JL

Water cycle in double precision (largescale+moistadj)
Improved wate rayleigh.
First step for rayleigh with variable species. Now, just need to change optcv.
changed some interpolation indices in callcorrk to limit dependency of OLR on the number of layers

File size: 11.2 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	DOUBLE PRECISION :: cp_new_t(nlayermx), v_cptt(ngrid,nlayermx)
59	REAL cp_delta_t(nlayermx)
60	DOUBLE PRECISION :: v_cptj(nlayermx), v_cptjk1, v_ssig
61	REAL 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	DOUBLE PRECISION :: gamcpdz(ngrid,2:nlayermx)
65	DOUBLE PRECISION :: zdp, zdpm
66
67	REAL zsat ! super-saturation
68	REAL zflo ! flotabilite
69
70	DOUBLE PRECISION :: 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.-zpsat(i,k-1)/pplay(i,k-1))zdpm + (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) = ( (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_pratio + v_zqs 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	ENDDO
213	Enddo ! nn=1,niter
214
215
216	!--------------------------------------------------- sounding downwards
217	! -- we refine the prognostic variables in
218	! -- the layer about to be adjusted
219
220	! DO k = k1, k2
221	! v_cptt(i,k) = RCPD * local_t(i,k)
222	! v_t = local_t(i,k)
223	! v_p = pplay(i,k)
224	!
225	! call Psat_water(v_t,v_p,zpsat,zqs(i,k))
226	! call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
227	! ENDDO
228
229	DO k = 2, nlayermx
230	zdpm = pplev(i,k-1) - pplev(i,k)
231	zdp = pplev(i,k) - pplev(i,k+1)
232	! gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)zdpm + v_cptt(i,k)zdp) &
233	! + RLVTT /(zdpm+zdp) * (zqs(i,k-1)zdpm + zqs(i,k)zdp) ) &
234	! * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) &
235	! / (1.0+ (zdqs(i,k-1)zdpm + zdqs(i,k)zdp)/(zdpm+zdp) )
236	! general case where water is not a trace gas (JL13)
237	v_zqs = (zqs(i,k-1)zdpm + zqs(i,k)zdp)/(zdpm+zdp)
238	v_cptt2 = (v_cptt(i,k-1)zdpm + v_cptt(i,k)zdp)/(zdpm+zdp)
239	v_pratio = ((1.-zpsat(i,k-1)/pplay(i,k-1))zdpm + (1.-zpsat(i,k)/pplay(i,k))zdp)/(zdpm+zdp)
240	v_dlnpsat = (zdlnpsat(i,k-1)zdpm + zdlnpsat(i,k)zdp)/(zdpm+zdp)
241	gamcpdz(i,k) = ( (R/RCPDv_cptt2(1.- v_zqs) + RLVTTv_zqs) (pplay(i,k-1)-pplay(i,k))/pplev(i,k) ) &
242	/ (((1.- v_zqs) + v_zqs * RCPV/RCPD)v_pratio + v_zqs v_dlnpsat)
243	ENDDO
244
245	! Test to see if we've reached the bottom
246
247	IF (k1 .EQ. 1) GOTO 841 ! yes we have!
248	zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
249	zsat=(local_q(i,k1-1)-zqs(i,k1-1))*(pplev(i,k1-1)-pplev(i,k1)) &
250	+ (local_q(i,k1)-zqs(i,k1))*(pplev(i,k1)-pplev(i,k1+1))
251	IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 841 ! yes we have!
252
253	840 CONTINUE
254	k1 = k1 - 1
255	IF (k1 .EQ. 1) GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
256	zsat = zsat + (local_q(i,k1-1)-zqs(i,k1-1)) &
257	*(pplev(i,k1-1)-pplev(i,k1))
258	zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
259	IF (zflo.GT.0.0 .AND. zsat.GT.0.0) THEN
260	GOTO 840
261	ELSE
262	GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
263	ENDIF
264	841 CONTINUE
265
266	GOTO 810 ! look for other layers higher up
267
268	9999 CONTINUE ! loop over all the points
269
270	!-----------------------------------------------------------------------
271	! Determine the cloud fraction (hypothese: la nebulosite a lieu
272	! a l'endroit ou la vapeur d'eau est diminuee par l'ajustement):
273
274	DO k = 1, nlayermx
275	DO i = 1, ngrid
276	IF (itest(i)) THEN
277	delta_q(i,k) = local_q(i,k) - zq(i,k)
278	IF (delta_q(i,k).LT.0.) rneb(i,k) = 1.0
279	ENDIF
280	ENDDO
281	ENDDO
282
283	! Distribuer l'eau condensee en eau liquide nuageuse (hypothese:
284	! l'eau liquide est distribuee aux endroits ou la vapeur d'eau
285	! diminue et d'une maniere proportionnelle a cet diminution):
286
287	DO i = 1, ngrid
288	IF (itest(i)) THEN
289	zq1(i) = 0.0
290	zq2(i) = 0.0
291	ENDIF
292	ENDDO
293	DO k = 1, nlayermx
294	DO i = 1, ngrid
295	IF (itest(i)) THEN
296	zdp = pplev(i,k)-pplev(i,k+1)
297	zq1(i) = zq1(i) - delta_q(i,k) * zdp
298	zq2(i) = zq2(i) - MIN(0.0, delta_q(i,k)) * zdp
299	ENDIF
300	ENDDO
301	ENDDO
302	DO k = 1, nlayermx
303	DO i = 1, ngrid
304	IF (itest(i)) THEN
305	IF (zq2(i).NE.0.0) d_ql(i,k) = - MIN(0.0,delta_q(i,k))*zq1(i)/zq2(i)
306	ENDIF
307	ENDDO
308	ENDDO
309
310	DO k = 1, nlayermx
311	DO i = 1, ngrid
312	local_q(i, k) = MAX(local_q(i, k), seuil_vap)
313	ENDDO
314	ENDDO
315
316	DO k = 1, nlayermx
317	DO i = 1, ngrid
318	d_t(i,k) = local_t(i,k) - zt(i,k)
319	d_q(i,k) = local_q(i,k) - zq(i,k)
320	ENDDO
321	ENDDO
322
323	! now subroutine -----> GCM variables
324	DO k = 1, nlayermx
325	DO i = 1, ngrid
326
327	pdtmana(i,k) = d_t(i,k)/ptimestep
328	pdqmana(i,k,i_h2o) = d_q(i,k)/ptimestep
329	pdqmana(i,k,i_ice) = d_ql(i,k)/ptimestep
330
331	ENDDO
332	ENDDO
333
334
335	RETURN
336	END

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

Download in other formats:

Original Format