Context Navigation

ch4cloud.F90 @ 3251

Last change on this file since 3251 was 3247, checked in by afalco, 2 years ago
Added missing files from pluto.old
File size: 11.8 KB

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

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

Download in other formats:

Original Format