Context Navigation

moistadj_generic.F90 @ 3205

Last change on this file since 3205 was 3104, checked in by emillour, 2 years ago
new routine moistadj_generic
File size: 14.0 KB

Line
1	subroutine moistadj_generic(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	use generic_cloud_common_h
6	use generic_tracer_index_mod, only: generic_tracer_index
7	use tracer_h
8	use ioipsl_getin_p_mod, only: getin_p !-> to get the metallicity
9	use comcstfi_mod, only: r, cpp, mugaz
10
11	implicit none
12
13
14	!=====================================================================
15	!
16	! Purpose
17	! -------
18	! Calculates moist convective adjustment by the method of Manabe.
19	!
20	! Authors
21	! -------
22	! Adapted from the moistadj.F90 routine
23	! for generic tracers (condensable species)
24	! by Noe CLEMENT (2023)
25	!
26	!=====================================================================
27
28	INTEGER,INTENT(IN) :: ngrid, nlayer, nq
29
30	REAL,INTENT(IN) :: pt(ngrid,nlayer) ! temperature (K)
31	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracer (kg/kg)
32	REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq)
33	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
34	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
35	REAL,INTENT(IN) :: ptimestep ! physics timestep (s)
36	REAL,INTENT(OUT) :: pdqmana(ngrid,nlayer,nq) ! tracer tendencies (kg/kg.s-1)
37	REAL,INTENT(OUT) :: pdtmana(ngrid,nlayer) ! temperature increment(K/s)
38	REAL,INTENT(OUT) :: rneb(ngrid,nlayer) ! cloud fraction
39
40	! Options :
41	!real, save :: metallicity ! metallicity of planet
42	REAL,SAVE :: metallicity = 0.0
43	!$OMP THREADPRIVATE(metallicity)
44
45	! local variables
46	REAL zt(ngrid,nlayer) ! temperature (K)
47	REAL zq(ngrid,nlayer) ! humidite specifique (kg/kg)
48
49	REAL d_t(ngrid,nlayer) ! temperature increment
50	REAL d_q(ngrid,nlayer) ! incrementation pour vapeur d'eau
51	REAL d_ql(ngrid,nlayer) ! incrementation pour l'eau liquide
52
53	! REAL t_coup
54	! PARAMETER (t_coup=234.0)
55	REAL seuil_vap
56	PARAMETER (seuil_vap=1.0E-10)
57
58	! Local variables
59
60	logical call_ice_vap_generic ! to call only one time the ice/vap pair of a tracer
61	integer igcm_generic_vap, igcm_generic_ice ! index of the vap and ice of generic_tracer
62
63	INTEGER i, k, iq, nn
64	INTEGER, PARAMETER :: niter = 1
65	INTEGER k1, k1p, k2, k2p
66	LOGICAL itest(ngrid)
67	REAL delta_q(ngrid, nlayer)
68	DOUBLE PRECISION :: cp_new_t(nlayer), v_cptt(ngrid,nlayer)
69	REAL cp_delta_t(nlayer)
70	DOUBLE PRECISION :: v_cptj(nlayer), v_cptjk1, v_ssig
71	REAL v_p, v_t, v_zqs,v_cptt2,v_pratio,v_dlnpsat
72	REAL zqs(ngrid,nlayer), zdqs(ngrid,nlayer),zpsat(ngrid,nlayer),zdlnpsat(ngrid,nlayer)
73	REAL zq1(ngrid), zq2(ngrid)
74	DOUBLE PRECISION :: gamcpdz(ngrid,2:nlayer)
75	DOUBLE PRECISION :: zdp, zdpm
76
77	real q_cri(ngrid,nlayer)
78
79	REAL zsat ! super-saturation
80	REAL zflo ! flotabilite
81
82	DOUBLE PRECISION :: local_q(ngrid,nlayer),local_t(ngrid,nlayer)
83
84	REAL zdelta, zcor, zcvm5
85
86	REAL dEtot, dqtot, masse ! conservation diagnostics
87	real dL1tot, dL2tot
88
89	! Indices of generic vapour and ice tracers
90	real,save :: RCPD=0.0
91	INTEGER,SAVE :: i_vap_generic=0 ! Generic Condensable Species vapour
92	INTEGER,SAVE :: i_ice_generic=0 ! Generic Condensable Species ice
93	!$OMP THREADPRIVATE(i_vap_generic,i_ice_generic,RCPD)
94
95	LOGICAL,SAVE :: firstcall=.TRUE.
96	!$OMP THREADPRIVATE(firstcall)
97
98	IF (firstcall) THEN
99
100	RCPD = cpp
101	! RCPV = 1.88e3 ! specific heat capacity of water vapor at 350K ! MUST BE CHANGED !!!
102	! 1.709e3 for methane - should be added in table_tracers_condensable (datagcm)
103
104	write(,) "value for metallicity? "
105	metallicity=0.0 ! default value
106	call getin_p("metallicity",metallicity)
107	write(,) " metallicity = ",metallicity
108
109	do iq=1, nq
110	call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
111	if(call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
112	write(,) "moistadjustment : taking back the values you have set in 'table_tracers_condensable' for thermodynamics. If you have forgotten one, it will crash."
113	m = constants_mass(iq)
114	delta_vapH = constants_delta_vapH(iq)
115	epsi_generic = constants_epsi_generic(iq)
116	RLVTT_generic = constants_RLVTT_generic(iq)
117	RCPV_generic = constants_RCPV_generic(iq)
118	endif
119	enddo
120
121	i_vap_generic=igcm_generic_vap
122	i_ice_generic=igcm_generic_ice
123
124	write(,) "rain: i_ice_generic=",i_ice_generic
125	write(,) " i_vap_generic=",i_vap_generic
126
127	firstcall = .FALSE.
128	ENDIF
129
130	! GCM -----> subroutine variables
131	zq(1:ngrid,1:nlayer) = pq(1:ngrid,1:nlayer,i_vap_generic)+ pdq(1:ngrid,1:nlayer,i_vap_generic)*ptimestep
132	zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer)
133	pdqmana(1:ngrid,1:nlayer,1:nq)=0.0
134
135	DO k = 1, nlayer
136	DO i = 1, ngrid
137	if(zq(i,k).lt.0.)then
138	zq(i,k)=0.0
139	endif
140	ENDDO
141	ENDDO
142
143	local_q(1:ngrid,1:nlayer) = zq(1:ngrid,1:nlayer)
144	local_t(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer)
145	rneb(1:ngrid,1:nlayer) = 0.0
146	d_ql(1:ngrid,1:nlayer) = 0.0
147	d_t(1:ngrid,1:nlayer) = 0.0
148	d_q(1:ngrid,1:nlayer) = 0.0
149
150	! Calculate v_cptt
151	DO k = 1, nlayer
152	DO i = 1, ngrid
153	v_cptt(i,k) = RCPD * local_t(i,k)
154	v_t = MAX(local_t(i,k),15.)
155	v_p = pplay(i,k)
156
157	call Psat_generic(v_t,v_p,metallicity,zpsat(i,k),zqs(i,k))
158	call Lcpdqsat_generic(v_t,v_p,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k))
159
160	!call Psat_water(v_t,v_p,zpsat(i,k),zqs(i,k))
161	!call Lcpdqsat_water(v_t,v_p,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k))
162	ENDDO
163	ENDDO
164
165	! Calculate Gamma * Cp * dz: (gamma is the critical gradient)
166	DO k = 2, nlayer
167	DO i = 1, ngrid
168	zdp = pplev(i,k)-pplev(i,k+1)
169	zdpm = pplev(i,k-1)-pplev(i,k)
170	! gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)zdpm + v_cptt(i,k)zdp) &
171	! + RLVTT /(zdpm+zdp) * (zqs(i,k-1)zdpm + zqs(i,k)zdp) ) &
172	!* (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) &
173	! / (1.0+ (zdqs(i,k-1)zdpm + zdqs(i,k)zdp)/(zdpm+zdp) )
174	! general case where water is not a trace gas (JL13)
175	v_zqs = (zqs(i,k-1)zdpm + zqs(i,k)zdp)/(zdpm+zdp)
176	v_cptt2 = (v_cptt(i,k-1)zdpm + v_cptt(i,k)zdp)/(zdpm+zdp)
177	v_pratio = ((1.-zpsat(i,k-1)/pplay(i,k-1))zdpm + (1.-zpsat(i,k)/pplay(i,k))zdp)/(zdpm+zdp)
178	v_dlnpsat = (zdlnpsat(i,k-1)zdpm + zdlnpsat(i,k)zdp)/(zdpm+zdp)
179	gamcpdz(i,k) = ( (R/RCPDv_cptt2(1.- v_zqs) + RLVTT_genericv_zqs) (pplay(i,k-1)-pplay(i,k))/pplev(i,k) ) &
180	/ (((1.- v_zqs) + v_zqs * RCPV_generic/RCPD)v_pratio + v_zqs v_dlnpsat)
181	ENDDO
182	ENDDO
183
184	! calculates moist convection inhibition criterion
185	DO k = 1, nlayer
186	DO i = 1, ngrid
187	q_cri(i,k) = ( 1 / (1 - 1/epsi_generic)) * r * mugaz/1000. / delta_vapH * zt(i,k)
188	ENDDO
189	ENDDO
190
191	write(,) "q_cri at 80K : ", ( 1 / (1 - 1/epsi_generic)) * (r * mugaz/1000.) / delta_vapH * 80
192
193	!------------------------------------ modification of unstable profile
194	DO 9999 i = 1, ngrid
195
196	itest(i) = .FALSE.
197
198	! print*,'we in the loop'
199	! stop
200
201	k1 = 0
202	k2 = 1
203
204	810 CONTINUE ! look for k1, the base of the column
205	k2 = k2 + 1
206	IF (k2 .GT. nlayer) GOTO 9999
207	zflo = v_cptt(i,k2-1) - v_cptt(i,k2) - gamcpdz(i,k2)
208	zsat=(local_q(i,k2-1)-zqs(i,k2-1))*(pplev(i,k2-1)-pplev(i,k2)) &
209	+(local_q(i,k2)-zqs(i,k2))*(pplev(i,k2)-pplev(i,k2+1))
210
211	IF ( zflo.LE.0.0 .OR. zsat.LE.0.0 .OR. local_q(i,k2-1).GT.q_cri(i,k2-1)) GOTO 810
212	k1 = k2 - 1
213	itest(i) = .TRUE.
214
215	820 CONTINUE !! look for k2, the top of the column
216	IF (k2 .EQ. nlayer) GOTO 821
217	k2p = k2 + 1
218	zsat=zsat+(pplev(i,k2p)-pplev(i,k2p+1))*(local_q(i,k2p)-zqs(i,k2p))
219	zflo = v_cptt(i,k2p-1) - v_cptt(i,k2p) - gamcpdz(i,k2p)
220
221	IF (zflo.LE.0.0 .OR. zsat.LE.0.0 .OR. local_q(i,k2p).GT.q_cri(i,k2p)) GOTO 821
222	k2 = k2p
223	GOTO 820
224	821 CONTINUE
225
226	!------------------------------------------------------ local adjustment
227	830 CONTINUE ! actual adjustment
228	Do nn=1,niter
229	v_cptj(k1) = 0.0
230	zdp = pplev(i,k1)-pplev(i,k1+1)
231	v_cptjk1 = ( (1.0+zdqs(i,k1))(v_cptt(i,k1)+v_cptj(k1)) + RLVTT_generic(local_q(i,k1)-zqs(i,k1)) ) * zdp
232	v_ssig = zdp * (1.0+zdqs(i,k1))
233
234	k1p = k1 + 1
235	DO k = k1p, k2
236	zdp = pplev(i,k)-pplev(i,k+1)
237	v_cptj(k) = v_cptj(k-1) + gamcpdz(i,k)
238	v_cptjk1 = v_cptjk1 + zdp * ( (1.0+zdqs(i, k))(v_cptt(i,k)+v_cptj(k)) + RLVTT_generic(local_q(i,k)-zqs(i,k)) )
239	v_ssig = v_ssig + zdp *(1.0+zdqs(i,k))
240	ENDDO
241
242
243	! this right here is where the adjustment is done???
244	DO k = k1, k2
245	cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k)
246	cp_delta_t(k) = cp_new_t(k) - v_cptt(i,k)
247	v_cptt(i,k)=cp_new_t(k)
248	local_q(i,k) = zqs(i,k) + zdqs(i,k)*cp_delta_t(k)/RLVTT_generic
249	local_t(i,k) = cp_new_t(k) / RCPD
250
251	v_t = MAX(local_t(i,k),15.)
252	v_p = pplay(i,k)
253
254	call Psat_generic(v_t,v_p,metallicity,zpsat(i,k),zqs(i,k))
255	call Lcpdqsat_generic(v_t,v_p,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k))
256
257	!call Psat_water(v_t,v_p,zpsat(i,k),zqs(i,k))
258	!call Lcpdqsat_water(v_t,v_p,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k))
259
260	ENDDO
261	Enddo ! nn=1,niter
262
263
264	!--------------------------------------------------- sounding downwards
265	! -- we refine the prognostic variables in
266	! -- the layer about to be adjusted
267
268	! DO k = k1, k2
269	! v_cptt(i,k) = RCPD * local_t(i,k)
270	! v_t = local_t(i,k)
271	! v_p = pplay(i,k)
272	!
273	! call Psat_water(v_t,v_p,zpsat,zqs(i,k))
274	! call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
275	! ENDDO
276
277	DO k = 2, nlayer
278	zdpm = pplev(i,k-1) - pplev(i,k)
279	zdp = pplev(i,k) - pplev(i,k+1)
280	! gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)zdpm + v_cptt(i,k)zdp) &
281	! + RLVTT /(zdpm+zdp) * (zqs(i,k-1)zdpm + zqs(i,k)zdp) ) &
282	! * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) &
283	! / (1.0+ (zdqs(i,k-1)zdpm + zdqs(i,k)zdp)/(zdpm+zdp) )
284	! general case where water is not a trace gas (JL13)
285	v_zqs = (zqs(i,k-1)zdpm + zqs(i,k)zdp)/(zdpm+zdp)
286	v_cptt2 = (v_cptt(i,k-1)zdpm + v_cptt(i,k)zdp)/(zdpm+zdp)
287	v_pratio = ((1.-zpsat(i,k-1)/pplay(i,k-1))zdpm + (1.-zpsat(i,k)/pplay(i,k))zdp)/(zdpm+zdp)
288	v_dlnpsat = (zdlnpsat(i,k-1)zdpm + zdlnpsat(i,k)zdp)/(zdpm+zdp)
289	gamcpdz(i,k) = ( (R/RCPDv_cptt2(1.- v_zqs) + RLVTT_genericv_zqs) (pplay(i,k-1)-pplay(i,k))/pplev(i,k) ) &
290	/ (((1.- v_zqs) + v_zqs * RCPV_generic/RCPD)v_pratio + v_zqs v_dlnpsat)
291	ENDDO
292
293	! Test to see if we've reached the bottom
294
295	IF (k1 .EQ. 1) GOTO 841 ! yes we have!
296	zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
297	zsat=(local_q(i,k1-1)-zqs(i,k1-1))*(pplev(i,k1-1)-pplev(i,k1)) &
298	+ (local_q(i,k1)-zqs(i,k1))*(pplev(i,k1)-pplev(i,k1+1))
299	IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 841 ! yes we have!
300
301	840 CONTINUE
302	k1 = k1 - 1
303	IF (k1 .EQ. 1) GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
304	zsat = zsat + (local_q(i,k1-1)-zqs(i,k1-1)) &
305	*(pplev(i,k1-1)-pplev(i,k1))
306	zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
307	IF (zflo.GT.0.0 .AND. zsat.GT.0.0) THEN
308	GOTO 840
309	ELSE
310	GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
311	ENDIF
312	841 CONTINUE
313
314	GOTO 810 ! look for other layers higher up
315
316	9999 CONTINUE ! loop over all the points
317
318	!-----------------------------------------------------------------------
319	! Determine the cloud fraction (hypothese: la nebulosite a lieu
320	! a l'endroit ou la vapeur d'eau est diminuee par l'ajustement):
321
322	DO k = 1, nlayer
323	DO i = 1, ngrid
324	IF (itest(i)) THEN
325	delta_q(i,k) = local_q(i,k) - zq(i,k)
326	IF (delta_q(i,k).LT.0.) rneb(i,k) = 1.0
327	ENDIF
328	ENDDO
329	ENDDO
330
331	! Distribuer l'eau condensee en eau liquide nuageuse (hypothese:
332	! l'eau liquide est distribuee aux endroits ou la vapeur d'eau
333	! diminue et d'une maniere proportionnelle a cet diminution):
334	DO i = 1, ngrid
335	IF (itest(i)) THEN
336	zq1(i) = 0.0
337	zq2(i) = 0.0
338	ENDIF
339	ENDDO
340	DO k = 1, nlayer
341	DO i = 1, ngrid
342	IF (itest(i)) THEN
343	zdp = pplev(i,k)-pplev(i,k+1)
344	zq1(i) = zq1(i) - delta_q(i,k) * zdp
345	zq2(i) = zq2(i) - MIN(0.0, delta_q(i,k)) * zdp
346	ENDIF
347	ENDDO
348	ENDDO
349	DO k = 1, nlayer
350	DO i = 1, ngrid
351	IF (itest(i)) THEN
352	IF (zq2(i).NE.0.0) d_ql(i,k) = - MIN(0.0,delta_q(i,k))*zq1(i)/zq2(i)
353	ENDIF
354	ENDDO
355	ENDDO
356
357	DO k = 1, nlayer
358	DO i = 1, ngrid
359	local_q(i, k) = MAX(local_q(i, k), seuil_vap)
360	ENDDO
361	ENDDO
362
363	DO k = 1, nlayer
364	DO i = 1, ngrid
365	d_t(i,k) = local_t(i,k) - zt(i,k)
366	d_q(i,k) = local_q(i,k) - zq(i,k)
367	ENDDO
368	ENDDO
369
370	! now subroutine -----> GCM variables
371	DO k = 1, nlayer
372	DO i = 1, ngrid
373
374	pdtmana(i,k) = d_t(i,k)/ptimestep
375	pdqmana(i,k,i_vap_generic) = d_q(i,k)/ptimestep
376	pdqmana(i,k,i_ice_generic) = d_ql(i,k)/ptimestep
377
378	ENDDO
379	ENDDO
380
381	end subroutine moistadj_generic

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

Download in other formats:

Original Format