Context Navigation

ch4cloud.F90 @ 3948

Last change on this file since 3948 was 3948, checked in by tbertrand, 6 weeks ago
PLUTO PCM : Simple scheme for CH4 clouds with number of nuclei determined by the haze number mixing ratio and not a fixed number in callphys.def (Nmix) TB
File size: 12.6 KB

Line
1	SUBROUTINE ch4cloud(ngrid,nlay,ptimestep, &
2	pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt, &
3	pq,pdq,pdqcloud,pdqscloud,pdtcloud, &
4	nq,rice_ch4)
5
6	use comgeomfi_h
7	use comcstfi_mod, only: pi, g, cpp
8	use tracer_h, only: igcm_ch4_gas, igcm_ch4_ice, &
9	rho_ch4_ice,lw_ch4,micro_indx
10	use callkeys_mod, only: Nmix_ch4, callmufi
11
12	IMPLICIT NONE
13
14	!=======================================================================
15	! Treatment of saturation of METHANE
16	!
17	!
18	! Modif de zq si saturation dans l'atmosphere
19	! si zq(ig,l)> zqsat(ig,l) -> zq(ig,l)=zqsat(ig,l)
20	! Le test est effectue de bas en haut. CH4 condensee
21	! (si saturation) est remise dans la couche en dessous.
22	! Le methane condensee dans la couche du bas est deposee a la surface
23	!
24	! Melanie Vangvichith (adapted from Mars water clouds)
25	! Tanguy Bertrand
26	! Completely rewritten by Francois Forget (to properly estimate the
27	! latent heat release.
28	!
29	!=======================================================================
30
31	!-----------------------------------------------------------------------
32	! declarations:
33	! -------------
34
35	! Inputs:
36	! ------
37
38	INTEGER ngrid,nlay
39	REAL ptimestep ! pas de temps physique (s)
40	REAL pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa)
41	REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
42	REAL pdpsrf(ngrid) ! tendance surf pressure
43	REAL pzlev(ngrid,nlay+1) ! altitude at layer boundaries
44	REAL pzlay(ngrid,nlay) ! altitude at the middle of the layers
45	REAL pt(ngrid,nlay) ! temperature at the middle of the layers (K)
46	REAL pdt(ngrid,nlay) ! tendance temperature des autres param.
47
48	real pq(ngrid,nlay,nq) ! traceur (kg/kg)
49	real pdq(ngrid,nlay,nq) ! tendance avant condensation (kg/kg.s-1)
50	integer nq ! nombre de traceurs
51
52	! Outputs:
53	! -------
54
55	real pdqcloud(ngrid,nlay,nq) ! tendance de la condensation CH4(kg/kg.s-1)
56	real pdqscloud(ngrid,nq) ! flux en surface (kg.m-2.s-1)
57	REAL pdtcloud(ngrid,nlay) ! tendance temperature due
58	! a la chaleur latente
59	REAL rice_ch4(ngrid,nlay) ! CH4 Ice mass mean radius (m)
60	! (r_c in montmessin_2004)
61
62	! local:
63	! ------
64	! REAL Nmix ! Cloud condensation nuclei
65	! parameter (Nmix=1.E5) ! /kg
66	real rnuclei ! Nuclei geometric mean radius (m)
67	parameter (rnuclei=2.E-7) ! m
68
69	REAL CBRT
70	EXTERNAL CBRT
71	INTEGER ig,l
72
73
74	REAL zq(ngrid,nlay,nq) ! local value of tracers
75	REAL zqsat(ngrid,nlay) ! saturation
76	REAL zt(ngrid,nlay) ! local value of temperature
77	REAL nmix(ngrid,nlay) ! Profil of Nmix
78
79	REAL vecnull(ngrid*nlay)
80
81	LOGICAL,SAVE :: firstcall=.true.
82
83	! indexes of methane gas, methane ice and dust tracers:
84	INTEGER,SAVE :: i_ch4=0 ! methane gas
85	INTEGER,SAVE :: i_ice=0 ! methane ice
86
87	REAL Tfin,qch4_fin
88	REAL temp_fin ! function used to compute Tfin at the end of the timestep
89
90	! TEMPORAIRE :
91
92	real pdtcloudmax,pdtcloudmin
93	integer igmin, igmax,lmin,lmax
94
95
96
97	! ** un petit test de coherence
98	! --------------------------
99
100	IF (firstcall) THEN
101	IF(ngrid.NE.ngrid) THEN
102	PRINT*,'STOP dans ch4cloud'
103	PRINT*,'probleme de dimensions :'
104	PRINT*,'ngrid =',ngrid
105	PRINT*,'ngrid =',ngrid
106	STOP
107	ENDIF
108
109	if (nq.gt.nq) then
110	write(,) 'stop in ch4cloud (nq.gt.nq)!'
111	write(,) 'nq=',nq,' nq=',nq
112	stop
113	endif
114
115	i_ch4=igcm_ch4_gas
116	i_ice=igcm_ch4_ice
117
118	if (i_ch4.eq.0) then
119	write(,) "stop in ch4cloud: i_ch4=0. "
120	write(,) "Maybe put ch4 in traceur.def?"
121	stop
122	endif
123
124	write(,) "methanecloud: i_ch4=",i_ch4
125	write(,) " i_ice=",i_ice
126
127	firstcall=.false.
128	ENDIF ! of IF (firstcall)
129
130
131	!-----------------------------------------------------------------------
132	! 1. initialisation
133	! -----------------
134
135	! On "update" la valeur de q(nq) (methane vapor) et temperature.
136	! On effectue qqes calculs preliminaires sur les couches :
137
138	zt(:,:)=pt(:,:)+ pdt(:,:)*ptimestep
139	zq(:,:,i_ch4)=pq(:,:,i_ch4)+pdq(:,:,i_ch4)*ptimestep
140	zq(:,:,i_ice)=pq(:,:,i_ice)+pdq(:,:,i_ice)*ptimestep
141
142	if (callmufi) then
143	zq(:,:,micro_indx(1))=pq(:,:,micro_indx(1))+ &
144	pdq(:,:,micro_indx(1))*ptimestep
145	zq(:,:,micro_indx(3))=pq(:,:,micro_indx(3))+ &
146	pdq(:,:,micro_indx(3))*ptimestep
147	endif
148
149	do l=1,nlay
150	do ig=1,ngrid
151	zq(ig,l,i_ch4)=max(zq(ig,l,i_ch4),1.E-30)
152	zq(ig,l,i_ice)=max(zq(ig,l,i_ice),0.)
153	enddo
154	enddo
155	if (callmufi) then
156	do l=1,nlay
157	do ig=1,ngrid
158	zq(ig,l,micro_indx(1))=max(zq(ig,l,micro_indx(1)),1.E-30)
159	zq(ig,l,micro_indx(3))=max(zq(ig,l,micro_indx(3)),1.E-30)
160	enddo
161	enddo
162	endif
163
164	pdqscloud(1:ngrid,1:nq)=0
165	pdqcloud(1:ngrid,1:nlay,1:nq)=0
166	pdtcloud(1:ngrid,1:nlay)=0
167
168	! Nmix profile
169	if (callmufi) then
170	nmix(:,:)=zq(ig,l,micro_indx(1))+zq(ig,l,micro_indx(3))
171	else
172	nmix(:,:)=Nmix_ch4
173	endif
174
175	! ----------------------------------------------
176	!
177	!
178	! q à saturation dans la couche l à temperature prédite zt
179	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
180	call methanesat(ngrid*nlay,zt,pplay,zqsat,vecnull)
181
182
183	! Loop to work where CH4 condensation/sublimation is occuring
184	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
185	do l=1,nlay
186	do ig=1,ngrid
187	if ((zq(ig,l,i_ch4).gt.zqsat(ig,l)).or. &
188	(zq(ig,l,i_ice).gt.1.E-10) ) then ! phase change
189
190	if(zq(ig,l,i_ice).lt.(zqsat(ig,l)-zq(ig,l,i_ch4))) then
191	!case when everything is sublimed:
192	qch4_fin = zq(ig,l,i_ch4) + zq(ig,l,i_ice)
193	Tfin = zt(ig,l) - (qch4_fin-zq(ig,l,i_ch4))*lw_ch4/cpp
194	else
195	!case when CH4 ice is present at the end
196	qch4_fin = zqsat(ig,l)
197	Tfin = zt(ig,l) - (qch4_fin-zq(ig,l,i_ch4))*lw_ch4/cpp
198	if(abs(Tfin-zt(ig,l)).gt.1.) then
199	! Improved calculation if the temperature change is significant (1K?)
200
201	! Estimating the temperature Tfin and condensation at the end of the
202	! timestep due to condensation-sublimation process, taking into acount
203	! latent heat, given by: Tfin - zt = (zq - qsat(Tfin))*L/cp
204	! Tfin, solution, of the equation, is estimated by iteration
205	! We assume that improved Tfin is between zt and the 1st estimation of Tfin:
206	Tfin = temp_fin(zt(ig,l),Tfin,0.01,pplay(ig,l), &
207	zt(ig,l),zq(ig,l,i_ch4),qch4_fin)
208	! Security check for small changes
209	if(abs(zq(ig,l,i_ch4)-qch4_fin).gt. &
210	abs(zq(ig,l,i_ch4)- zqsat(ig,l)) ) then
211	write(,) "warning: inaccuracy in CH4 clouds"
212	write(,) 'zq(ig,l,i_ch4),qch4_fin, zqsat(ig,l)', &
213	zq(ig,l,i_ch4),qch4_fin, zqsat(ig,l)
214
215	qch4_fin = zqsat(ig,l)
216
217	end if
218	end if
219	end if
220	! Final tendancies
221	pdqcloud(ig,l,i_ch4)=(qch4_fin-zq(ig,l,i_ch4))/ptimestep
222
223	!TB16
224	IF (zq(ig,l,i_ch4)+pdqcloud(ig,l,i_ch4)*ptimestep.lt.0.) &
225	then
226	pdqcloud(ig,l,i_ch4)=-zq(ig,l,i_ch4)/ptimestep
227	END IF
228
229	pdqcloud(ig,l,i_ice) = - pdqcloud(ig,l,i_ch4)
230	pdtcloud(ig,l)=(Tfin - zt(ig,l))/ptimestep
231
232	! Ice particle radius
233	zq(ig,l,i_ice)= &
234	zq(ig,l,i_ice)+pdqcloud(ig,l,i_ice)*ptimestep
235	rice_ch4(ig,l)=max( CBRT ( (zq(ig,l,i_ice)/rho_ch4_ice &
236	+nmix(ig,l)(4./3.)pirnuclei3.)/(nmix(ig,l)4./3.*pi)), &
237	rnuclei)
238	end if
239	enddo ! of do ig=1,ngrid
240	enddo ! of do l=1,nlay
241
242	! A correction if a lot of subliming N2 fills the 1st layer FF04/2005
243	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
244	! Then that should not affect the ice particle radius
245	do ig=1,ngrid
246	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,2)))then
247	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,3))) &
248	rice_ch4(ig,2)=rice_ch4(ig,3)
249	rice_ch4(ig,1)=rice_ch4(ig,2)
250	end if
251	end do
252
253	!********** ANALYSE pour Icarus **************
254	!
255	! pdtcloudmax =0.
256	! pdtcloudmin =0.
257	! igmin=999
258	! igmax=999
259	! lmin =999
260	! lmax =999
261	! do l=1, nlay
262	! do ig=1,ngrid
263	! if(pdtcloud(ig,l).gt.pdtcloudmax) then
264	! igmax=ig
265	! lmax = l
266	! pdtcloudmax=pdtcloud(ig,l)
267	! else if(pdtcloud(ig,l).lt.pdtcloudmin) then
268	! igmin=ig
269	! lmin = l
270	! pdtcloudmin=pdtcloud(ig,l)
271	! end if
272	! end do
273	! end do
274	!
275	! write(,)
276	! write(,) "MAX , ig, l, Tbefore, Tafter, Tfin_ini"
277	! write(,) igmax,lmax, zt(igmax,lmax),
278	! & zt(igmax,lmax)+ pdtcloud(igmax,lmax)*ptimestep,
279	! &zt(igmax,lmax)-(zqsat(igmax,lmax)-zq(igmax,lmax,i_ch4))*lw_ch4/cpp
280	! write(,) "MAX , Qch4gas_before, Qch4gas_after, Qsat ini"
281	! write(,) pq(igmax,lmax,i_ch4)+pdq(igmax,lmax,i_ch4)*ptimestep,
282	! & pq(igmax,lmax,i_ch4)+
283	! & (pdq(igmax,lmax,i_ch4)+ pdqcloud(igmax,lmax,i_ch4))*ptimestep
284	! & ,zqsat(igmax,lmax)
285	! write(,) "MAX , Qch4ice_before, Qch4ice_after"
286	! write(,) pq(igmax,lmax,i_ice)+pdq(igmax,lmax,i_ice)*ptimestep,
287	! & pq(igmax,lmax,i_ice)+
288	! & (pdq(igmax,lmax,i_ice)+ pdqcloud(igmax,lmax,i_ice))*ptimestep
289	! write(,) 'Rice=', rice_ch4(igmax,lmax)
290	!
291	! write(,)
292	! write(,) "MIN , ig, l, Tbefore, Tafter, Tfin_ini"
293	! write(,) igmin,lmin, zt(igmin,lmin),
294	! & zt(igmin,lmin)+ pdtcloud(igmin,lmin)*ptimestep,
295	! &zt(igmin,lmin)-(zqsat(igmin,lmin)-zq(igmin,lmin,i_ch4))*lw_ch4/cpp
296	! write(,) "MIN , Qch4gas_before, Qch4gas_after, Qsat ini"
297	! write(,) pq(igmin,lmin,i_ch4)+pdq(igmin,lmin,i_ch4)*ptimestep,
298	! & pq(igmin,lmin,i_ch4)+
299	! & (pdq(igmin,lmin,i_ch4)+ pdqcloud(igmin,lmin,i_ch4))*ptimestep
300	! & ,zqsat(igmin,lmin)
301	! write(,) "MIN , Qch4ice_before, Qch4ice_after"
302	! write(,) pq(igmin,lmin,i_ice)+pdq(igmin,lmin,i_ice)*ptimestep,
303	! & pq(igmin,lmin,i_ice)+
304	! & (pdq(igmin,lmin,i_ice)+ pdqcloud(igmin,lmin,i_ice))*ptimestep
305	! write(,) 'Rice=', rice_ch4(igmin,lmin)
306	! write(,) "pdtcloud(igmin,lmin) " , pdtcloud(igmin,lmin)
307	! write(,) "pdqcloud(i_ch4)", pdqcloud(igmin,lmin,i_ch4)
308	! write(,) "pdqcloud(i_ice)", pdqcloud(igmin,lmin,i_ice)
309	!
310
311	!********** FIN ANALYSE pour Icarus **************
312
313	!**************************************************
314	! Output --- removed
315	!**************************************************
316	! NB: for diagnostics use zq(), the updated value of tracers
317	! Computing ext visible optical depth in each layer
318
319	RETURN
320	END
321
322	!================================================================================
323
324	FUNCTION temp_fin(X1,X2,XACC,pres,temp,zq,qsat)
325	use comcstfi_mod, only: pi, g, cpp
326	use tracer_h, only: igcm_ch4_gas, igcm_ch4_ice, rho_ch4_ice, lw_ch4
327	implicit none
328
329	real pres, temp, zq,qsat
330	real X1,X2,XACC, F
331	real FMID,D,temp_fin,DX, XMID
332	integer JMAX,J
333	real func
334
335	PARAMETER (JMAX=40)
336
337	call methanesat(1,X2,pres,qsat,0.)
338	FMID = X2 -temp -(zq-qsat)*lw_ch4/cpp
339
340	call methanesat(1,X1,pres,qsat,0.)
341	F = X1 -temp -(zq-qsat)*lw_ch4/cpp
342
343	IF(FFMID.GE.0.) write(,*) 'Fix Tfin firt guesses in ch4cloud'
344	IF(F.LT.0.)THEN
345	temp_fin=X1
346	DX=X2-X1
347	ELSE
348	temp_fin=X2
349	DX=X1-X2
350	ENDIF
351	DO 11 J=1,JMAX
352	DX=DX*.5
353	XMID=temp_fin+DX
354
355	call methanesat(1,XMID,pres,qsat,0.)
356	FMID = XMID -temp -(zq-qsat)*lw_ch4/cpp
357
358	IF(FMID.LE.0.)temp_fin=XMID
359	IF(ABS(DX).LT.XACC .OR. FMID.EQ.0.) RETURN
360	11 CONTINUE
361	! PAUSE 'too many bisections'
362	END
363
364

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

Download in other formats:

Original Format