Context Navigation

moistadj.F90 @ 778

Last change on this file since 778 was 773, checked in by jleconte, 13 years ago

05/09/2012 == JL

Correction of the calculation of the solar longitude in tlocked case.

-Can now handle any prograde resonance with nres=omega_rot/omega_orb.
-Sun now goes westward for the standard 2:1 case, as expected.

In the gray case, the separation between kappa_IR and VI is now set by

wave number, independently of the usual IR/VISIBLE calculation separation.
i.e. kappa_IR can be used in the calculation of the downward stellar flux

if the wavenumber in the band is low enough and vice versa.

In ave_stelspec, stellar flux averaging has been generalized to incorporate

very red/blue stellar spectra (great care must however be taken of the band
limit used for the corralated k distributions).

-Brown dwarf spectra from Allard et al. have been added.
-Any Black body temperature can now be used.

File size: 10.2 KB

Line
1	subroutine moistadj(pt, pq, pdq, pplev, pplay, pdtmana, pdqmana, ptimestep, rneb)
2
3	use watercommon_h, only: T_h2O_ice_liq, RLVTT, RCPD, Psat_water, Lcpdqsat_water
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	REAL pt(ngridmx,nlayermx) ! temperature (K)
27	REAL pq(ngridmx,nlayermx,nqmx) ! tracer (kg/kg)
28	REAL pdq(ngridmx,nlayermx,nqmx)
29
30	REAL pdqmana(ngridmx,nlayermx,nqmx) ! tendency of tracers (kg/kg.s-1)
31	REAL pdtmana(ngridmx,nlayermx) ! temperature increment
32
33	! local variables
34	REAL zt(ngridmx,nlayermx) ! temperature (K)
35	REAL zq(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, nn
52	INTEGER, PARAMETER :: niter = 1
53	INTEGER k1, k1p, k2, k2p
54	LOGICAL itest(ngridmx)
55	REAL delta_q(ngridmx, nlayermx)
56	REAL cp_new_t(nlayermx)
57	REAL cp_delta_t(nlayermx)
58	REAL v_cptj(nlayermx), v_cptjk1, v_ssig
59	REAL v_cptt(ngridmx,nlayermx), v_p, v_t, zpsat
60	REAL zqs(ngridmx,nlayermx), zdqs(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	zq(1:ngridmx,1:nlayermx) = pq(1:ngridmx,1:nlayermx,i_h2o)+ pdq(1:ngridmx,1:nlayermx,i_h2o)*ptimestep
97	zt(1:ngridmx,1:nlayermx) = pt(1:ngridmx,1:nlayermx)
98	pdqmana(1:ngridmx,1:nlayermx,1:nqmx)=0.0
99
100	DO k = 1, nlayermx
101	DO i = 1, ngridmx
102	if(zq(i,k).lt.0.)then
103	zq(i,k)=0.0
104	endif
105	ENDDO
106	ENDDO
107
108	local_q(1:ngridmx,1:nlayermx) = zq(1:ngridmx,1:nlayermx)
109	local_t(1:ngridmx,1:nlayermx) = zt(1:ngridmx,1:nlayermx)
110	rneb(1:ngridmx,1:nlayermx) = 0.0
111	d_ql(1:ngridmx,1:nlayermx) = 0.0
112	d_t(1:ngridmx,1:nlayermx) = 0.0
113	d_q(1:ngridmx,1:nlayermx) = 0.0
114
115	! Calculate v_cptt
116	DO k = 1, nlayermx
117	DO i = 1, ngridmx
118	v_cptt(i,k) = RCPD * local_t(i,k)
119	v_t = MAX(local_t(i,k),15.)
120	v_p = pplay(i,k)
121
122	call Psat_water(v_t,v_p,zpsat,zqs(i,k))
123	call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
124	ENDDO
125	ENDDO
126
127	! Calculate Gamma * Cp * dz: (gamma is the critical gradient)
128	DO k = 2, nlayermx
129	DO i = 1, ngridmx
130	zdp = pplev(i,k)-pplev(i,k+1)
131	zdpm = pplev(i,k-1)-pplev(i,k)
132	gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)zdpm + v_cptt(i,k)zdp) &
133	+ RLVTT /(zdpm+zdp) * (zqs(i,k-1)zdpm + zqs(i,k)zdp) ) &
134	* (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) &
135	/ (1.0+ (zdqs(i,k-1)zdpm + zdqs(i,k)zdp)/(zdpm+zdp) )
136	ENDDO
137	ENDDO
138
139	!------------------------------------ modification of unstable profile
140	DO 9999 i = 1, ngridmx
141	itest(i) = .FALSE.
142
143	! print*,'we in the loop'
144	! stop
145
146	k1 = 0
147	k2 = 1
148
149	810 CONTINUE ! look for k1, the base of the column
150	k2 = k2 + 1
151	IF (k2 .GT. nlayermx) GOTO 9999
152	zflo = v_cptt(i,k2-1) - v_cptt(i,k2) - gamcpdz(i,k2)
153	zsat=(local_q(i,k2-1)-zqs(i,k2-1))*(pplev(i,k2-1)-pplev(i,k2)) &
154	+(local_q(i,k2)-zqs(i,k2))*(pplev(i,k2)-pplev(i,k2+1))
155
156	IF ( zflo.LE.0.0 .OR. zsat.LE.0.0 ) GOTO 810
157	k1 = k2 - 1
158	itest(i) = .TRUE.
159
160	820 CONTINUE !! look for k2, the top of the column
161	IF (k2 .EQ. nlayermx) GOTO 821
162	k2p = k2 + 1
163	zsat=zsat+(pplev(i,k2p)-pplev(i,k2p+1))*(local_q(i,k2p)-zqs(i,k2p))
164	zflo = v_cptt(i,k2p-1) - v_cptt(i,k2p) - gamcpdz(i,k2p)
165
166	IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 821
167	k2 = k2p
168	GOTO 820
169	821 CONTINUE
170
171	!------------------------------------------------------ local adjustment
172	830 CONTINUE ! actual adjustment
173	Do nn=1,niter
174	v_cptj(k1) = 0.0
175	zdp = pplev(i,k1)-pplev(i,k1+1)
176	v_cptjk1 = ( (1.0+zdqs(i,k1))(v_cptt(i,k1)+v_cptj(k1)) + RLVTT(local_q(i,k1)-zqs(i,k1)) ) * zdp
177	v_ssig = zdp * (1.0+zdqs(i,k1))
178
179	k1p = k1 + 1
180	DO k = k1p, k2
181	zdp = pplev(i,k)-pplev(i,k+1)
182	v_cptj(k) = v_cptj(k-1) + gamcpdz(i,k)
183	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)) )
184	v_ssig = v_ssig + zdp *(1.0+zdqs(i,k))
185	ENDDO
186
187
188	! this right here is where the adjustment is done???
189	DO k = k1, k2
190	cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k)
191	cp_delta_t(k) = cp_new_t(k) - v_cptt(i,k)
192	v_cptt(i,k)=cp_new_t(k)
193	local_q(i,k) = zqs(i,k) + zdqs(i,k)*cp_delta_t(k)/RLVTT
194	local_t(i,k) = cp_new_t(k) / RCPD
195
196	v_t = MAX(local_t(i,k),15.)
197	v_p = pplay(i,k)
198
199	call Psat_water(v_t,v_p,zpsat,zqs(i,k))
200	call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
201
202
203
204	! print*,'i,k,zqs,cpdT=',i,k,zqs(i,k),cp_delta_t(k)
205	ENDDO
206	Enddo ! nn=1,niter
207
208
209	!--------------------------------------------------- sounding downwards
210	! -- we refine the prognostic variables in
211	! -- the layer about to be adjusted
212
213	! DO k = k1, k2
214	! v_cptt(i,k) = RCPD * local_t(i,k)
215	! v_t = local_t(i,k)
216	! v_p = pplay(i,k)
217	!
218	! call Psat_water(v_t,v_p,zpsat,zqs(i,k))
219	! call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
220	! ENDDO
221
222	DO k = 2, nlayermx
223	zdpm = pplev(i,k-1) - pplev(i,k)
224	zdp = pplev(i,k) - pplev(i,k+1)
225	gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)zdpm + v_cptt(i,k)zdp) &
226	+ RLVTT /(zdpm+zdp) * (zqs(i,k-1)zdpm + zqs(i,k)zdp) ) &
227	* (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) &
228	/ (1.0+ (zdqs(i,k-1)zdpm + zdqs(i,k)zdp)/(zdpm+zdp) )
229	ENDDO
230
231	! Test to see if we've reached the bottom
232
233	IF (k1 .EQ. 1) GOTO 841 ! yes we have!
234	zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
235	zsat=(local_q(i,k1-1)-zqs(i,k1-1))*(pplev(i,k1-1)-pplev(i,k1)) &
236	+ (local_q(i,k1)-zqs(i,k1))*(pplev(i,k1)-pplev(i,k1+1))
237	IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 841 ! yes we have!
238
239	840 CONTINUE
240	k1 = k1 - 1
241	IF (k1 .EQ. 1) GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
242	zsat = zsat + (local_q(i,k1-1)-zqs(i,k1-1)) &
243	*(pplev(i,k1-1)-pplev(i,k1))
244	zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
245	IF (zflo.GT.0.0 .AND. zsat.GT.0.0) THEN
246	GOTO 840
247	ELSE
248	GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
249	ENDIF
250	841 CONTINUE
251
252	GOTO 810 ! look for other layers higher up
253
254	9999 CONTINUE ! loop over all the points
255
256	! print*,'k1=',k1
257	! print*,'k2=',k2
258
259	! print*,'local_t=',local_t
260	! print*,'v_cptt=',v_cptt
261	! print*,'gamcpdz=',gamcpdz
262
263	!-----------------------------------------------------------------------
264	! Determine the cloud fraction (hypothese: la nebulosite a lieu
265	! a l'endroit ou la vapeur d'eau est diminuee par l'ajustement):
266
267	DO k = 1, nlayermx
268	DO i = 1, ngridmx
269	IF (itest(i)) THEN
270	delta_q(i,k) = local_q(i,k) - zq(i,k)
271	IF (delta_q(i,k).LT.0.) rneb(i,k) = 1.0
272	ENDIF
273	ENDDO
274	ENDDO
275
276	! Distribuer l'eau condensee en eau liquide nuageuse (hypothese:
277	! l'eau liquide est distribuee aux endroits ou la vapeur d'eau
278	! diminue et d'une maniere proportionnelle a cet diminution):
279
280	DO i = 1, ngridmx
281	IF (itest(i)) THEN
282	zq1(i) = 0.0
283	zq2(i) = 0.0
284	ENDIF
285	ENDDO
286	DO k = 1, nlayermx
287	DO i = 1, ngridmx
288	IF (itest(i)) THEN
289	zdp = pplev(i,k)-pplev(i,k+1)
290	zq1(i) = zq1(i) - delta_q(i,k) * zdp
291	zq2(i) = zq2(i) - MIN(0.0, delta_q(i,k)) * zdp
292	ENDIF
293	ENDDO
294	ENDDO
295	DO k = 1, nlayermx
296	DO i = 1, ngridmx
297	IF (itest(i)) THEN
298	IF (zq2(i).NE.0.0) d_ql(i,k) = - MIN(0.0,delta_q(i,k))*zq1(i)/zq2(i)
299	ENDIF
300	ENDDO
301	ENDDO
302
303	! print*,'local_q BEFORE=',local_q
304
305	DO k = 1, nlayermx
306	DO i = 1, ngridmx
307	local_q(i, k) = MAX(local_q(i, k), seuil_vap)
308	ENDDO
309	ENDDO
310
311	DO k = 1, nlayermx
312	DO i = 1, ngridmx
313	d_t(i,k) = local_t(i,k) - zt(i,k)
314	d_q(i,k) = local_q(i,k) - zq(i,k)
315	ENDDO
316	ENDDO
317
318	! now subroutine -----> GCM variables
319	DO k = 1, nlayermx
320	DO i = 1, ngridmx
321
322	pdtmana(i,k) = d_t(i,k)/ptimestep
323	pdqmana(i,k,i_h2o) = d_q(i,k)/ptimestep
324	pdqmana(i,k,i_ice) = d_ql(i,k)/ptimestep
325
326	ENDDO
327	ENDDO
328
329
330	RETURN
331	END

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

Download in other formats:

Original Format