Context Navigation

moistadj.F90 @ 1351

Last change on this file since 1351 was 1315, checked in by milmd, 10 years ago
LMDZ.GENERIC. OpenMP directives added in generic physic. When running in pure OpenMP or hybrid OpenMP/MPI, may have some bugs with condense_cloud and wstats routines.
File size: 11.3 KB

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

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

Download in other formats:

Original Format