Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

watercloud.F @ 3094

Last change on this file since 3094 was 11, checked in by aslmd, 14 years ago
spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 17.1 KB

Line
1	SUBROUTINE watercloud(ngrid,nlay, ptimestep,
2	& pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt,
3	& pq,pdq,pdqcloud,qsurf,pdqscloud,pdtcloud,
4	& nq,nsize,tau,icount,zls)
5	IMPLICIT NONE
6
7	c=======================================================================
8	c Treatment of saturation of water vapor
9	c
10	c
11	c Modif de zq si saturation dans l'atmopshere
12	c si zq(ig,l)> zqsat(ig,l) -> zq(ig,l)=zqsat(ig,l)
13	c Le test est effectue de bas en haut. L'eau condensee
14	c (si saturation) est remise dans la couche en dessous.
15	c L'eau condensee dans la couche du bas est deposee a la surface
16	c
17	c WARNING : water vapor mixing ratio is assumed to be
18	c q(nqmx)
19	c
20	c Modification: Franck Montmessin water ice scheme
21	c
22	c=======================================================================
23
24	c-----------------------------------------------------------------------
25	c declarations:
26	c -------------
27
28	#include "dimensions.h"
29	#include "dimphys.h"
30	#include "comcstfi.h"
31	#include "callkeys.h"
32	#include "tracer.h"
33
34	#include "fisice.h"
35	#include "comgeomfi.h"
36	c
37	c arguments:
38	c ----------
39
40	INTEGER ngrid,nlay,nsize
41	INTEGER icount
42	REAL ptimestep ! pas de temps physique (s)
43	REAL pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa)
44	REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
45	REAL pzlev(ngrid,nlay+1) ! altitude at layer boundaries
46	REAL pzlay(ngrid,nlay) ! altitude at the middle of the layers
47	REAL pdpsrf(ngrid) ! tendance surf pressure
48	REAL pt(ngrid,nlay) ! temperature au centre des couches (K)
49	REAL pdt(ngrid,nlay) ! tendance temperature des autres param.
50	REAL pdtcloud(ngrid,nlay) ! tendance temperature due a la chaleur latente
51	REAL tau(ngridmx,nsize)
52	REAL zls
53
54	c Traceurs :
55	integer nq ! nombre de traceurs
56	real pq(ngrid,nlay,nq) ! traceur (kg/kg)
57	real pdq(ngrid,nlay,nq) ! tendance avant condensation (kg/kg.s-1)
58	real pdqcloud(ngrid,nlay,nq) ! tendance de la condesation H2O(kg/kg.s-1)
59	real qsurf(ngrid,nq)
60	real pdqscloud(ngrid,nq) ! flux en surface (kg.m-2.s-1)
61
62	c local:
63	c ------
64
65	REAL CBRT
66	EXTERNAL CBRT
67	INTEGER ig,l
68
69	LOGICAL firstcall,improved
70	SAVE firstcall
71
72	REAL zq(ngridmx,nlayermx,nqmx)
73	REAL zq0(ngridmx,nlayermx,nqmx)
74	REAL zqsat(ngridmx,nlayermx)
75	REAL zt(ngridmx,nlayermx)
76
77
78	REAL masse (ngridmx,nlayermx)
79	REAL epaisseur (ngridmx,nlayermx)
80	REAL dustcores(ngridmx,nlayermx) !Dust number density (#/kg)
81	REAL rfinal ! Ice crystal radius after condensation(m)
82	REAL seq ! Equilibrium saturation ration (accounting for curvature effect)
83	REAL dzq ! masse de glace echangee (kg/kg)
84	REAL lw !Latent heat of sublimation (J.kg-1)
85	REAL To
86
87	REAL Ctot
88	REAL*8 ph2o,satu
89	REAL gr,Cste,up,dwn,newvap
90
91	DATA firstcall/.true./
92	cc DATA improved/.false./ ! Actionne une microphysique plus raffinee
93	DATA improved/.true./
94
95	c Reference temperature, T=273,15 K
96	data To/273.15/
97
98	SAVE improved
99
100	c POur diagnostique :
101	c ~~~~~~~~~~~~~~~~~
102	REAL taucond(ngridmx,nlayermx) ! taux de condensation (kg/kg/s-1)
103	ccc****WRF
104	c REAL mtot(ngridmx) ! Total mass of water vapor (kg/m2)
105	c REAL icetot(ngridmx) ! Total mass of water ice (kg/m2)
106	c REAL rave(ngridmx) ! Mean crystal radius in a column (m)
107	ccc****WRF
108	REAL cloudmass(ngridmx,nlayermx) ! Integrated mass of water ice in each layer (kg/m2)
109	REAL waterhem(2,3)
110
111	REAL icetot2pm(ngridmx)
112	REAL mtot2pm(ngridmx)
113	REAL rave2pm(ngridmx)
114	REAL qsurf2pm(ngridmx)
115	INTEGER i,nit
116	REAL ecart,loctime
117	integer univtime,lon2pm(ngridmx)
118
119	integer ig_vl2
120	real latvl2,lonvl2
121
122	logical tesfield ! output of TES like data
123	data tesfield /.false./
124
125	c ** un petit test de coherence
126	c --------------------------
127
128	IF (firstcall) THEN
129	IF(ngrid.NE.ngridmx) THEN
130	PRINT*,'STOP dans watercloud'
131	PRINT*,'probleme de dimensions :'
132	PRINT*,'ngrid =',ngrid
133	PRINT*,'ngridmx =',ngridmx
134	STOP
135	ENDIF
136	firstcall=.false.
137	if (tesfield) then
138	latvl2= 47.57
139	lonvl2= 134.26
140	ig_vl2= 1+ int( (1.5-(latvl2-90.)jjm/180.) -2 )iim +
141	& int(1.5+(lonvl2+180)*iim/360.)
142	open(201,file='VL2_data')
143	open(199,file='water_budgets')
144	print*,'File water_budgets created'
145	c Find the universal time corresponding to 2PM (TES obs)
146	print*,'Universal time corresponding to 2PM (TES obs)'
147	print*,'of each longitude - Used when iceparty is on.'
148	print*,'2PM outputs are written in a specific file.'
149	print*,'\| Longitude \| U.T = 2PM LT \|'
150	nit= 48
151	do ig=1,ngridmx
152	ecart=24.
153	do i=1,nit-1
154	loctime = i/float(nit)24. + long(ig)12./pi
155	if (loctime.gt.24.) loctime = loctime - 24.
156	if (loctime.lt.0) loctime = loctime + 24.
157	if (abs(14.-loctime).lt.ecart) then
158	ecart = abs(14.-loctime)
159	lon2pm(ig) = i
160	endif
161	enddo
162	if(int(lati(ig)*180./pi).eq.0)
163	* print,long(ig)/pi180.,'\|',lon2pm(ig)/2.
164	enddo
165	open(200,file='2pm_outputs',form='unformatted')
166	print*,'File 2pm_outputs created'
167	write(200) long,lati
168	endif ! (tesfield)
169	ENDIF
170
171
172	c-----------------------------------------------------------------------
173	c 1. initialisation
174	c -----------------
175
176	c On "update" la valeur de q(nqmx) (water vapor) et temperature.
177	c On effectue qqes calculs preliminaires sur les couches :
178	c masse (kg.m-2), epaisseur(m).
179
180	do l=1,nlay
181	do ig=1,ngrid
182	zq(ig,l,nq)=pq(ig,l,nq)+pdq(ig,l,nq)*ptimestep
183	zq(ig,l,nq)=max(zq(ig,l,nq),1.E-30) ! FF 12/2004
184	zq0(ig,l,nq)=zq(ig,l,nq)
185	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
186	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
187	epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l)
188	if(iceparty) then
189	zq(ig,l,nq-1)=pq(ig,l,nq-1)+pdq(ig,l,nq-1)*ptimestep
190	zq(ig,l,nq-1)=max(zq(ig,l,nq-1),0.) ! FF 12/2004
191	zq0(ig,l,nq-1)=zq(ig,l,nq-1)
192	rdust(ig,l)= max(.8e-6*exp(-pzlay(ig,l)/18000.),1.e-9)
193	c Calcul du rayon moyen des particules de glace.
194	c Hypothese : Dans une couche, la glace presente se
195	c repartit uniformement autour du nbre de poussieres
196	c dont le rayon moyen est prescrit par rdust.
197	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
198	if ((dustbin.eq.1).and..not.(doubleq)) then
199	zq(ig,l,1)=pq(ig,l,1)+pdq(ig,l,1)*ptimestep
200	dustcores(ig,l)=max(zq(ig,l,1)/
201	& (rho_dust4./3.piradius(1)*3.),1.e-9)
202	else
203	dustcores(ig,l)=( epaisseur(ig,l)/masse(ig,l) ) *
204	& 2.e+6/0.1max(tau(ig,1),0.001)exp(-pzlay(ig,l)/10000.)
205
206	!!!PATCH WATER
207	c TEMPORAIRE : réduction du nombre de nuclei FF 04/2008 :
208	dustcores(ig,l) = dustcores(ig,l) / 8.
209	! reduction facteur 2
210
211	endif
212	rice(ig,l)=CBRT( ( zq(ig,l,nq-1)/rho_ice+
213	& dustcores(ig,l)(4./3.)pirdust(ig,l)*3. )
214	& / (dustcores(ig,l)4./3.pi) )
215	rice(ig,l)=max(rice(ig,l),rdust(ig,l))
216	endif
217	enddo
218	enddo
219
220	call zerophys(ngrid*nq,pdqscloud)
221	call zerophys(ngridnlaynq,pdqcloud)
222	call zerophys(ngrid*nlay,pdtcloud)
223	call zerophys(ngrid,mtot)
224
225	call zerophys(ngrid,icetot)
226	call zerophys(ngrid,rave)
227	call zerophys(ngrid,cloudmass)
228
229	c ----------------------------------------------
230	c
231	c
232	c Rapport de melange a saturation dans la couche l : -------
233	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
234
235	call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
236
237	c taux de condensation (kg/kg/s-1) dans les differentes couches
238	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
239	c (Pour diagnostic seulement !)
240	if(.not.iceparty)then
241	do l=1, nlay
242	do ig=1,ngridmx
243	taucond(ig,l)= max((zq(ig,l,nq)-zqsat(ig,l))/ptimestep,0.)
244	end do
245	end do
246	endif
247
248
249	if(iceparty) then
250	do l=1,nlay
251	do ig=1,ngrid
252
253	If (improved) then
254	c Improved microphysics scheme
255	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
256
257	Ctot = zq(ig,l,nq) + zq(ig,l,nq-1)
258	ph2o = zq(ig,l,nq) * 44. / 18. * pplay(ig,l)
259	satu = zq(ig,l,nq) / zqsat(ig,l)
260
261	call growthrate(ptimestep,zt(ig,l),pplay(ig,l),
262	& ph2o,ph2o/satu,seq,rice(ig,l),gr)
263	Cste = ptimestep * 4. * pi * rice(ig,l)
264	* * rho_ice * dustcores(ig,l)
265	up = zq(ig,l,nq) + Cste * gr * seq
266	dwn = 1. + Cste * gr / zqsat(ig,l)
267	newvap = min(up/dwn,Ctot)
268
269	gr = gr * ( newvap/zqsat(ig,l) - seq )
270	dzq = min( max( Cste * gr,-zq(ig,l,nq-1) )
271	* , zq(ig,l,nq) )
272
273	c Nucleation (sat ratio must be larger than a critical value)
274	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
275	if (satu.gt.1.) then
276	if (satu.le.1.4.and.zq(ig,l,nq-1).lt.1.e-8)
277	* dzq = 0.
278	endif
279
280	Else
281	c Old version
282	c ~~~~~~~~~~~
283	if (zq(ig,l,nq).ge.zqsat(ig,l))then ! Condensation
284	dzq=zq(ig,l,nq)-zqsat(ig,l)
285	elseif(zq(ig,l,nq).lt.zqsat(ig,l))then ! Sublimation
286	dzq=-min(zqsat(ig,l)-zq(ig,l,nq),zq(ig,l,nq-1))
287	endif
288
289	Endif
290
291	c Water Mass change
292	c ~~~~~~~~~~~~~~~~~
293	zq(ig,l,nq-1)=zq(ig,l,nq-1)+dzq
294	zq(ig,l,nq)=zq(ig,l,nq)-dzq
295
296	rice(ig,l)=max( CBRT ( (zq(ig,l,nq-1)/rho_ice
297	& +dustcores(ig,l)(4./3.)pirdust(ig,l)*3.)
298	& /(dustcores(ig,l)4./3.pi) ), rdust(ig,l))
299
300	if(activice.and.pclc(ig).gt.0)
301	& rice(ig,l)=rice(ig,l)/CBRT(pclc(ig))
302
303	enddo
304	enddo
305
306	else ! if not iceparty
307
308	c Saturation couche nlay a 2 :
309	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
310	do l=nlay,2, -1
311	do ig=1,ngrid
312	if (zq(ig,l,nq).gt.zqsat(ig,l))then
313	zq(ig,l-1,nq) = zq(ig,l-1,nq)+(zq(ig,l,nq)-zqsat(ig,l))
314	& *(pplev(ig,l)-pplev(ig,l+1))/(pplev(ig,l-1)-pplev(ig,l))
315	zq(ig,l,nq)=zqsat(ig,l)
316	endif
317	enddo
318	enddo
319
320	c Saturation couche l=1 si pas iceparty
321	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
322	do ig=1,ngridmx
323	if (zq(ig,1,nq).gt.zqsat(ig,1))then
324	pdqscloud(ig,nq)= (zq(ig,1,nq)-zqsat(ig,1))
325	& (pplev(ig,1)-pplev(ig,2))/(gptimestep)
326	zq(ig,1,nq)=zqsat(ig,1)
327	endif
328	end do
329
330	endif ! (iceparty)
331
332	c Tendance finale
333	c ~~~~~~~~~~~~~~~
334	do l=1, nlay
335	do ig=1,ngridmx
336	pdqcloud(ig,l,nq) = (zq(ig,l,nq) - zq0(ig,l,nq))/ptimestep
337	if(iceparty) then
338	pdqcloud(ig,l,nq-1) =
339	& (zq(ig,l,nq-1) - zq0(ig,l,nq-1))/ptimestep
340	endif
341	lw=(2834.3-0.28(zt(ig,l)-To)-0.004(zt(ig,l)-To)*2)1.e+3
342	pdtcloud(ig,l)=-pdqcloud(ig,l,nq)*lw/cpp
343	end do
344	end do
345
346	c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
347	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
348	c Then that should not affect the ice particle radius
349	do ig=1,ngridmx
350	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,2)))then
351	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,3)))
352	& rice(ig,2)=rice(ig,3)
353	rice(ig,1)=rice(ig,2)
354	end if
355	end do
356
357	c**************************************************
358	c Output
359	c**************************************************
360	univtime = icount - int(float(icount)/48.) * 48
361
362	c if "tesfield"
363	c ~~~~~~~~~~~~~
364	if (tesfield.and.mod(icount,12).eq.0)
365	. call zerophys(6,waterhem)
366	c****WRF: ligne trop longue
367	do ig=1,ngridmx
368	do l=1 ,nlay
369	c masse de vapeur d'eau dans la couche l
370	mtot(ig)=mtot(ig)+pq(ig,l,nq)*masse(ig,l)
371	if(iceparty) then
372	c masse de glace d'eau dans la couche l
373	icetot(ig)=icetot(ig)+masse(ig,l)*pq(ig,l,nq-1)
374	c rayon moyen des cristaux dans la colonne ig
375	rave(ig)=rave(ig)+masse(ig,l)pq(ig,l,nq-1)rice(ig,l)
376	cloudmass(ig,l)=masse(ig,l)*pq(ig,l,nq-1)
377	endif
378	enddo
379	if(iceparty) then
380	rave(ig)=rave(ig)/max(icetot(ig),1.e-30)
381	if (icetot(ig)*1000.lt.0.01) rave(ig)=0.
382	endif
383	if (tesfield) then
384	if (univtime.eq.lon2pm(ig)) mtot2pm(ig) = mtot(ig)
385	if (univtime.eq.lon2pm(ig).and.iceparty)
386	* icetot2pm(ig) = icetot(ig)
387	if (univtime.eq.lon2pm(ig).and.iceparty)
388	* rave2pm(ig) = rave(ig)
389	if (univtime.eq.lon2pm(ig)) qsurf2pm(ig) = qsurf(ig,nqmx)
390	endif
391
392	c Outputs total mass of water and clouds *
393	c ******************************************
394	IF (tesfield.and.mod(icount,12).EQ.0) THEN
395	if (int(lati(ig)*180/pi).gt.0) then
396	waterhem(1,1)=waterhem(1,1)+area(ig)*mtot(ig)
397	waterhem(1,2)=waterhem(1,2)+area(ig)*qsurf(ig,nqmx)
398	if (iceparty)waterhem(1,3)=waterhem(1,3)+area(ig)*icetot(ig)
399	endif
400	if (int(lati(ig)*180/pi).lt.0) then
401	waterhem(2,1)=waterhem(2,1)+area(ig)*mtot(ig)
402	waterhem(2,2)=waterhem(2,2)+area(ig)*qsurf(ig,nqmx)
403	if (iceparty)waterhem(2,3)=waterhem(2,3)+area(ig)*icetot(ig)
404	endif
405	if (int(lati(ig)*180/pi).eq.0) then
406	waterhem(1,1)=waterhem(1,1)+area(ig)*mtot(ig)/2.
407	waterhem(2,1)=waterhem(2,1)+area(ig)*mtot(ig)/2.
408	waterhem(1,2)=waterhem(1,2)+area(ig)*qsurf(ig,nqmx)/2.
409	waterhem(2,2)=waterhem(2,2)+area(ig)*qsurf(ig,nqmx)/2.
410	if (iceparty) then
411	waterhem(1,3)=waterhem(1,3)+area(ig)*icetot(ig)/2.
412	waterhem(2,3)=waterhem(2,3)+area(ig)*icetot(ig)/2.
413	endif
414	endif
415	ENDIF ! (total mass)
416	enddo ! (ngridmx)
417
418	if (tesfield.and.mod(icount,12).eq.0) then
419	write(199,6000)zls*180./pi,waterhem(1,1),waterhem(1,2),
420	& waterhem(1,3),waterhem(2,1),waterhem(2,2),waterhem(2,3)
421	!6000 format(f8.3,x,6(e12.5,x))
422	!!!!!!!!!!!!!!!!g95
423	6000 format(f8.3,1x,6(e12.5,1x))
424	write(201,6001) zls*180./pi,pplev(ig_vl2,1),icetot(ig_vl2)
425	& ,zt(ig_vl2,1)
426	!6001 format(f8.3,x,3(e12.5,x))
427	!!!!!!!!!!!!!!!!g95
428	6001 format(f8.3,1x,3(e12.5,1x))
429	endif
430
431	***************************************************
432	* Write outputs at 2PM LT in a dedicated file
433	***************************************************
434	if (tesfield) then
435	if (mod(icount,48).eq.0.and.icount.ne.1.and.iceparty) then
436	print*,'ecriture dans 2pm_outputs'
437	write(200) zls
438	write(200) mtot2pm
439	if (iceparty) write(200) icetot2pm
440	write(200) qsurf2pm
441	if (iceparty) write(200) rave2pm
442	endif
443	endif
444
445	c ******************************************
446	c Output in diagfi
447	c ******************************************
448
449	if (ngrid.gt.1) then
450
451	c CALL WRITEDIAGFI(ngridmx,'mtot',
452	c & 'total mass of water vapor','kg/m2',2,mtot)
453	c
454	c if (callstats) then
455	c call wstats(ngrid,"mtot",
456	c . "total mass of water vapor","kg/m2",2,mtot)
457	c endif
458	c
459	cc if(.not.iceparty)
460	cc & call WRITEDIAGFI(ngridmx,'taucond',
461	cc & 'taux cond H2O ice','kg/kg/s',3,taucond)
462	cc
463	c if(iceparty) then
464	c CALL WRITEDIAGFI(ngridmx,'icetot',
465	c & 'total mass of water ice','kg/m2',2,icetot)
466	c if (callstats) then
467	c call wstats(ngrid,"icetot",
468	c . "total mass of water ice","kg/m2",2,icetot)
469	c endif
470	cc CALL WRITEDIAGFI(ngridmx,'cloudmass',
471	cc & 'mass of water ice/layer','kg/m2',3,cloudmass)
472	cc CALL WRITEDIAGFI(ngridmx,'rice','ice radius',
473	cc & 'meter',3,rice)
474	cc CALL WRITEDIAGFI(ngridmx,'rave','Mean ice radius',
475	cc & 'meter',2,rave)
476	c endif
477
478	else
479
480	c CALL WRITEG1D(ngridmx,1,mtot,'gas','kg/m2')
481	c CALL WRITEG1D(ngridmx,1,icetot,'ice','kg/m2')
482	c call WRITEG1D(ngridmx,nlay,rice,'rice','um')
483	c call WRITEG1D(ngridmx,nlay,dustcores,'pouss','#/layer')
484
485	endif
486	99 continue
487
488	RETURN
489	END

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

Download in other formats:

Original Format