Context Navigation

ch4cloud.F90 @ 3504

Last change on this file since 3504 was 3252, checked in by afalco, 9 months ago
Pluto PCM: Fixed 1d ch4 and co AF
File size: 12.0 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	if (i_ch4.eq.0) then
118	write(,) "stop in ch4cloud: i_ch4=0. "
119	write(,) "Maybe put ch4 in traceur.def?"
120	stop
121	endif
122
123	write(,) "methanecloud: i_ch4=",i_ch4
124	write(,) " i_ice=",i_ice
125
126	firstcall=.false.
127	ENDIF ! of IF (firstcall)
128
129
130	!-----------------------------------------------------------------------
131	! 1. initialisation
132	! -----------------
133
134	! On "update" la valeur de q(nq) (methane vapor) et temperature.
135	! On effectue qqes calculs preliminaires sur les couches :
136
137	do l=1,nlay
138	do ig=1,ngrid
139	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
140	zq(ig,l,i_ch4)=pq(ig,l,i_ch4)+pdq(ig,l,i_ch4)*ptimestep
141	zq(ig,l,i_ch4)=max(zq(ig,l,i_ch4),1.E-30)
142	zq(ig,l,i_ice)=pq(ig,l,i_ice)+pdq(ig,l,i_ice)*ptimestep
143	zq(ig,l,i_ice)=max(zq(ig,l,i_ice),0.)
144	enddo
145	enddo
146
147	pdqscloud(1:ngrid,1:nq)=0
148	pdqcloud(1:ngrid,1:nlay,1:nq)=0
149	pdtcloud(1:ngrid,1:nlay)=0
150
151	! ----------------------------------------------
152	!
153	!
154	! q à saturation dans la couche l à temperature prédite zt
155	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
156	call methanesat(ngrid*nlay,zt,pplay,zqsat,vecnull)
157
158
159	! Loop to work where CH4 condensation/sublimation is occuring
160	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
161	do l=1,nlay
162	do ig=1,ngrid
163	if ((zq(ig,l,i_ch4).gt.zqsat(ig,l)).or. &
164	(zq(ig,l,i_ice).gt.1.E-10) ) then ! phase change
165
166	if(zq(ig,l,i_ice).lt.(zqsat(ig,l)-zq(ig,l,i_ch4))) then
167	!case when everything is sublimed:
168	qch4_fin = zq(ig,l,i_ch4) + zq(ig,l,i_ice)
169	Tfin = zt(ig,l) - (qch4_fin-zq(ig,l,i_ch4))*lw_ch4/cpp
170	else
171	!case when CH4 ice is present at the end
172	qch4_fin = zqsat(ig,l)
173	Tfin = zt(ig,l) - (qch4_fin-zq(ig,l,i_ch4))*lw_ch4/cpp
174	if(abs(Tfin-zt(ig,l)).gt.1.) then
175	! Improved calculation if the temperature change is significant (1K?)
176
177	! Estimating the temperature Tfin and condensation at the end of the
178	! timestep due to condensation-sublimation process, taking into acount
179	! latent heat, given by: Tfin - zt = (zq - qsat(Tfin))*L/cp
180	! Tfin, solution, of the equation, is estimated by iteration
181	! We assume that improved Tfin is between zt and the 1st estimation of Tfin:
182	Tfin = temp_fin(zt(ig,l),Tfin,0.01,pplay(ig,l), &
183	zt(ig,l),zq(ig,l,i_ch4),qch4_fin)
184	! Security check for small changes
185	if(abs(zq(ig,l,i_ch4)-qch4_fin).gt. &
186	abs(zq(ig,l,i_ch4)- zqsat(ig,l)) ) then
187	write(,) "warning: inaccuracy in CH4 clouds"
188	write(,) 'zq(ig,l,i_ch4),qch4_fin, zqsat(ig,l)', &
189	zq(ig,l,i_ch4),qch4_fin, zqsat(ig,l)
190
191	qch4_fin = zqsat(ig,l)
192
193	end if
194	end if
195	end if
196	! Final tendancies
197	pdqcloud(ig,l,i_ch4)=(qch4_fin-zq(ig,l,i_ch4))/ptimestep
198
199	!TB16
200	IF (zq(ig,l,i_ch4)+pdqcloud(ig,l,i_ch4)*ptimestep.lt.0.) &
201	then
202	pdqcloud(ig,l,i_ch4)=-zq(ig,l,i_ch4)/ptimestep
203	END IF
204
205	pdqcloud(ig,l,i_ice) = - pdqcloud(ig,l,i_ch4)
206	pdtcloud(ig,l)=(Tfin - zt(ig,l))/ptimestep
207
208	! Ice particle radius
209	zq(ig,l,i_ice)= &
210	zq(ig,l,i_ice)+pdqcloud(ig,l,i_ice)*ptimestep
211	rice_ch4(ig,l)=max( CBRT ( (zq(ig,l,i_ice)/rho_ch4_ice &
212	+Nmix_ch4(4./3.)pirnuclei3.)/(Nmix_ch44./3.*pi)), &
213	rnuclei)
214	end if
215	enddo ! of do ig=1,ngrid
216	enddo ! of do l=1,nlay
217
218	! A correction if a lot of subliming N2 fills the 1st layer FF04/2005
219	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
220	! Then that should not affect the ice particle radius
221	do ig=1,ngrid
222	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,2)))then
223	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,3))) &
224	rice_ch4(ig,2)=rice_ch4(ig,3)
225	rice_ch4(ig,1)=rice_ch4(ig,2)
226	end if
227	end do
228
229	!********** ANALYSE pour Icarus **************
230	!
231	! pdtcloudmax =0.
232	! pdtcloudmin =0.
233	! igmin=999
234	! igmax=999
235	! lmin =999
236	! lmax =999
237	! do l=1, nlay
238	! do ig=1,ngrid
239	! if(pdtcloud(ig,l).gt.pdtcloudmax) then
240	! igmax=ig
241	! lmax = l
242	! pdtcloudmax=pdtcloud(ig,l)
243	! else if(pdtcloud(ig,l).lt.pdtcloudmin) then
244	! igmin=ig
245	! lmin = l
246	! pdtcloudmin=pdtcloud(ig,l)
247	! end if
248	! end do
249	! end do
250	!
251	! write(,)
252	! write(,) "MAX , ig, l, Tbefore, Tafter, Tfin_ini"
253	! write(,) igmax,lmax, zt(igmax,lmax),
254	! & zt(igmax,lmax)+ pdtcloud(igmax,lmax)*ptimestep,
255	! &zt(igmax,lmax)-(zqsat(igmax,lmax)-zq(igmax,lmax,i_ch4))*lw_ch4/cpp
256	! write(,) "MAX , Qch4gas_before, Qch4gas_after, Qsat ini"
257	! write(,) pq(igmax,lmax,i_ch4)+pdq(igmax,lmax,i_ch4)*ptimestep,
258	! & pq(igmax,lmax,i_ch4)+
259	! & (pdq(igmax,lmax,i_ch4)+ pdqcloud(igmax,lmax,i_ch4))*ptimestep
260	! & ,zqsat(igmax,lmax)
261	! write(,) "MAX , Qch4ice_before, Qch4ice_after"
262	! write(,) pq(igmax,lmax,i_ice)+pdq(igmax,lmax,i_ice)*ptimestep,
263	! & pq(igmax,lmax,i_ice)+
264	! & (pdq(igmax,lmax,i_ice)+ pdqcloud(igmax,lmax,i_ice))*ptimestep
265	! write(,) 'Rice=', rice_ch4(igmax,lmax)
266	!
267	! write(,)
268	! write(,) "MIN , ig, l, Tbefore, Tafter, Tfin_ini"
269	! write(,) igmin,lmin, zt(igmin,lmin),
270	! & zt(igmin,lmin)+ pdtcloud(igmin,lmin)*ptimestep,
271	! &zt(igmin,lmin)-(zqsat(igmin,lmin)-zq(igmin,lmin,i_ch4))*lw_ch4/cpp
272	! write(,) "MIN , Qch4gas_before, Qch4gas_after, Qsat ini"
273	! write(,) pq(igmin,lmin,i_ch4)+pdq(igmin,lmin,i_ch4)*ptimestep,
274	! & pq(igmin,lmin,i_ch4)+
275	! & (pdq(igmin,lmin,i_ch4)+ pdqcloud(igmin,lmin,i_ch4))*ptimestep
276	! & ,zqsat(igmin,lmin)
277	! write(,) "MIN , Qch4ice_before, Qch4ice_after"
278	! write(,) pq(igmin,lmin,i_ice)+pdq(igmin,lmin,i_ice)*ptimestep,
279	! & pq(igmin,lmin,i_ice)+
280	! & (pdq(igmin,lmin,i_ice)+ pdqcloud(igmin,lmin,i_ice))*ptimestep
281	! write(,) 'Rice=', rice_ch4(igmin,lmin)
282	! write(,) "pdtcloud(igmin,lmin) " , pdtcloud(igmin,lmin)
283	! write(,) "pdqcloud(i_ch4)", pdqcloud(igmin,lmin,i_ch4)
284	! write(,) "pdqcloud(i_ice)", pdqcloud(igmin,lmin,i_ice)
285	!
286
287	!********** FIN ANALYSE pour Icarus **************
288
289	!**************************************************
290	! Output --- removed
291	!**************************************************
292	! NB: for diagnostics use zq(), the updated value of tracers
293	! Computing ext visible optical depth in each layer
294
295	RETURN
296	END
297
298	!================================================================================
299
300	FUNCTION temp_fin(X1,X2,XACC,pres,temp,zq,qsat)
301	use comcstfi_mod, only: pi, g, cpp
302	use tracer_h, only: igcm_ch4_gas, igcm_ch4_ice, rho_ch4_ice, lw_ch4
303	implicit none
304
305	real pres, temp, zq,qsat
306	real X1,X2,XACC, F
307	real FMID,D,temp_fin,DX, XMID
308	integer JMAX,J
309	real func
310
311	PARAMETER (JMAX=40)
312
313	call methanesat(1,X2,pres,qsat,0.)
314	FMID = X2 -temp -(zq-qsat)*lw_ch4/cpp
315
316	call methanesat(1,X1,pres,qsat,0.)
317	F = X1 -temp -(zq-qsat)*lw_ch4/cpp
318
319	IF(FFMID.GE.0.) write(,*) 'Fix Tfin firt guesses in ch4cloud'
320	IF(F.LT.0.)THEN
321	temp_fin=X1
322	DX=X2-X1
323	ELSE
324	temp_fin=X2
325	DX=X1-X2
326	ENDIF
327	DO 11 J=1,JMAX
328	DX=DX*.5
329	XMID=temp_fin+DX
330
331	call methanesat(1,XMID,pres,qsat,0.)
332	FMID = XMID -temp -(zq-qsat)*lw_ch4/cpp
333
334	IF(FMID.LE.0.)temp_fin=XMID
335	IF(ABS(DX).LT.XACC .OR. FMID.EQ.0.) RETURN
336	11 CONTINUE
337	! PAUSE 'too many bisections'
338	END
339
340

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

Download in other formats:

Original Format