Context Navigation

watercloud.F @ 47

Last change on this file since 47 was 38, checked in by emillour, 14 years ago
Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs. EM
File size: 11.9 KB

Line
1	SUBROUTINE watercloud(ngrid,nlay, ptimestep,
2	& pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt,
3	& pq,pdq,pdqcloud,pdqscloud,pdtcloud,
4	& nq,naersize,tau,
5	& ccn,rdust,rice,nuice)
6	IMPLICIT NONE
7
8	c=======================================================================
9	c Treatment of saturation of water vapor
10	c
11	c
12	c Modif de zq si saturation dans l'atmosphere
13	c si zq(ig,l)> zqsat(ig,l) -> zq(ig,l)=zqsat(ig,l)
14	c Le test est effectue de bas en haut. L'eau condensee
15	c (si saturation) est remise dans la couche en dessous.
16	c L'eau condensee dans la couche du bas est deposee a la surface
17	c
18	c Modification: Franck Montmessin water ice scheme
19	c Francois Forget : change nuclei density & outputs
20	c Ehouarn Millour: sept.2008, tracers are now handled
21	c by name (and not fixed index)
22	c
23	c=======================================================================
24
25	c-----------------------------------------------------------------------
26	c declarations:
27	c -------------
28
29	#include "dimensions.h"
30	#include "dimphys.h"
31	#include "comcstfi.h"
32	#include "callkeys.h"
33	#include "tracer.h"
34	#include "comgeomfi.h"
35
36	c Inputs:
37	c ------
38
39	INTEGER ngrid,nlay
40	REAL ptimestep ! pas de temps physique (s)
41	REAL pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa)
42	REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
43	REAL pdpsrf(ngrid) ! tendance surf pressure
44	REAL pzlev(ngrid,nlay+1) ! altitude at layer boundaries
45	REAL pzlay(ngrid,nlay) ! altitude at the middle of the layers
46	REAL pt(ngrid,nlay) ! temperature at the middle of the layers (K)
47	REAL pdt(ngrid,nlay) ! tendance temperature des autres param.
48
49	real pq(ngrid,nlay,nq) ! traceur (kg/kg)
50	real pdq(ngrid,nlay,nq) ! tendance avant condensation (kg/kg.s-1)
51
52	integer nq ! nombre de traceurs
53	integer naersize ! nombre de traceurs radiativement actifs (=naerkind)
54	REAL tau(ngridmx,naersize)
55	REAL ccn(ngridmx,nlayermx) ! Cloud condensation nuclei
56	! (particules kg-1)
57	c CCN reduction factor
58	REAL, PARAMETER :: ccn_factor = 1.
59	real rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)
60
61	c Outputs:
62	c -------
63
64	real pdqcloud(ngrid,nlay,nq) ! tendance de la condensation H2O(kg/kg.s-1)
65	real pdqscloud(ngrid,nq) ! flux en surface (kg.m-2.s-1)
66	REAL pdtcloud(ngrid,nlay) ! tendance temperature due
67	! a la chaleur latente
68
69	REAL rice(ngrid,nlay) ! Ice mass mean radius (m)
70	! (r_c in montmessin_2004)
71	REAL nuice(ngrid,nlay) ! Estimated effective variance
72	! of the size distribution
73
74	c local:
75	c ------
76
77	REAL CBRT
78	EXTERNAL CBRT
79	INTEGER ig,iq,l
80
81
82	REAL zq(ngridmx,nlayermx,nqmx) ! local value of tracers
83	REAL zq0(ngridmx,nlayermx,nqmx) ! local initial value of tracers
84	REAL zqsat(ngridmx,nlayermx) ! saturation
85	REAL zt(ngridmx,nlayermx) ! local value of temperature
86
87	REAL masse (ngridmx,nlayermx)
88	REAL epaisseur (ngridmx,nlayermx)
89	REAL rfinal ! Ice crystal radius after condensation(m)
90	REAL seq ! Equilibrium saturation ration (accounting for curvature effect)
91	REAL dzq ! masse de glace echangee (kg/kg)
92	REAL lw !Latent heat of sublimation (J.kg-1)
93	REAL,PARAMETER :: To=273.15 ! reference temperature, T=273.15 K
94
95	REAL Ctot
96	REAL*8 ph2o,satu
97	REAL gr,Cste,up,dwn,newvap
98
99	LOGICAL,SAVE :: firstcall=.true.
100	! To use more refined microphysics, set improved to .true.
101	LOGICAL,PARAMETER :: improved=.true.
102
103	c Pour diagnostique :
104	c ~~~~~~~~~~~~~~~~~
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
108	INTEGER i
109
110	! indexes of water vapour, water ice and dust tracers:
111	INTEGER,SAVE :: i_h2o=0 ! water vapour
112	INTEGER,SAVE :: i_ice=0 ! water ice
113	CHARACTER(LEN=20) :: tracername ! to temporarly store text
114
115	c ** un petit test de coherence
116	c --------------------------
117
118	IF (firstcall) THEN
119	IF(ngrid.NE.ngridmx) THEN
120	PRINT*,'STOP dans watercloud'
121	PRINT*,'probleme de dimensions :'
122	PRINT*,'ngrid =',ngrid
123	PRINT*,'ngridmx =',ngridmx
124	STOP
125	ENDIF
126
127	if (nq.gt.nqmx) then
128	write(,) 'stop in watercloud (nq.gt.nqmx)!'
129	write(,) 'nq=',nq,' nqmx=',nqmx
130	stop
131	endif
132
133	i_h2o=igcm_h2o_vap
134	i_ice=igcm_h2o_ice
135
136	write(,) "watercloud: i_h2o=",i_h2o
137	write(,) " i_ice=",i_ice
138
139	write(,) "water_param CCN reduction factor:", ccn_factor
140
141	firstcall=.false.
142	ENDIF ! of IF (firstcall)
143
144
145	c-----------------------------------------------------------------------
146	c 1. initialisation
147	c -----------------
148
149	c On "update" la valeur de q(nqmx) (water vapor) et temperature.
150	c On effectue qqes calculs preliminaires sur les couches :
151	c masse (kg.m-2), epaisseur(m).
152
153	do l=1,nlay
154	do ig=1,ngrid
155	zq(ig,l,i_h2o)=pq(ig,l,i_h2o)+pdq(ig,l,i_h2o)*ptimestep
156	zq(ig,l,i_h2o)=max(zq(ig,l,i_h2o),1.E-30) ! FF 12/2004
157	zq0(ig,l,i_h2o)=zq(ig,l,i_h2o)
158	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
159	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
160	epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l)
161
162	zq(ig,l,i_ice)=pq(ig,l,i_ice)+pdq(ig,l,i_ice)*ptimestep
163	zq(ig,l,i_ice)=max(zq(ig,l,i_ice),0.) ! FF 12/2004
164	zq0(ig,l,i_ice)=zq(ig,l,i_ice)
165
166	c This typical profile is not used anymore; rdust is
167	c set up in updatereffrad.F.
168	c rdust(ig,l)= max(.8e-6*exp(-pzlay(ig,l)/18000.),1.e-9)
169	enddo
170	enddo
171
172	do l=1,nlay
173	do ig=1,ngrid
174	ccn(ig,l) = ccn(ig,l) / ccn_factor
175	c TEMPORAIRE : réduction du nombre de nuclei FF 04/2008 :
176	c reduction facteur 3
177	c ccn(ig,l) = ccn(ig,l) / 27.
178	c reduction facteur 2
179	c ccn(ig,l) = ccn(ig,l) / 8.
180	c Calcul du rayon moyen des particules de glace.
181	c Hypothese : Dans une couche, la glace presente se
182	c repartit uniformement autour du nbre de poussieres
183	c dont le rayon moyen est prescrit par rdust.
184	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
185	rice(ig,l)=CBRT( ( zq(ig,l,i_ice)/rho_ice+
186	& ccn(ig,l)(4./3.)pirdust(ig,l)*3. )
187	& / (ccn(ig,l)4./3.pi) )
188	rice(ig,l)=max(rice(ig,l),rdust(ig,l))
189	c Effective variance of the size distribution
190	nuice(ig,l)=nuice_ref
191	enddo ! of do ig=1,ngrid
192	enddo ! of dol=1,nlay
193
194	pdqscloud(1:ngrid,1:nq)=0
195	pdqcloud(1:ngrid,1:nlay,1:nq)=0
196	pdtcloud(1:ngrid,1:nlay)=0
197
198	c icetot(1:ngrid)=0
199	c rave(1:ngrid)=0
200
201	c ----------------------------------------------
202	c
203	c
204	c Rapport de melange a saturation dans la couche l : -------
205	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
206
207	call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
208
209	c taux de condensation (kg/kg/s-1) dans les differentes couches
210	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
211
212	c Iceparty is not used anymore: water=>iceparty (JBM).
213	c if(iceparty) then
214
215	do l=1,nlay
216	do ig=1,ngrid
217
218	IF (improved) then
219	c Improved microphysics scheme
220	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
221
222	Ctot = zq(ig,l,i_h2o) + zq(ig,l,i_ice)
223	ph2o = zq(ig,l,i_h2o) * 44. / 18. * pplay(ig,l)
224	satu = zq(ig,l,i_h2o) / zqsat(ig,l)
225
226	call growthrate(ptimestep,zt(ig,l),pplay(ig,l),
227	& ph2o,ph2o/satu,seq,rice(ig,l),gr)
228	Cste = ptimestep * 4. * pi * rice(ig,l)
229	* * rho_ice * ccn(ig,l)
230	up = zq(ig,l,i_h2o) + Cste * gr * seq
231	dwn = 1. + Cste * gr / zqsat(ig,l)
232	newvap = min(up/dwn,Ctot)
233
234	gr = gr * ( newvap/zqsat(ig,l) - seq )
235	dzq = min( max( Cste * gr,-zq(ig,l,i_ice) )
236	* , zq(ig,l,i_h2o) )
237
238	c Nucleation (sat ratio must be larger than a critical value)
239	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
240	if (satu.gt.1.) then
241	if (satu.le.1.4.and.zq(ig,l,i_ice).lt.1.e-8)
242	* dzq = 0.
243	endif
244
245	ELSE
246	c Old version
247	c ~~~~~~~~~~~
248	if (zq(ig,l,i_h2o).ge.zqsat(ig,l))then ! Condensation
249	dzq=zq(ig,l,i_h2o)-zqsat(ig,l)
250	elseif(zq(ig,l,i_h2o).lt.zqsat(ig,l))then ! Sublimation
251	dzq=-min(zqsat(ig,l)-zq(ig,l,i_h2o),zq(ig,l,i_ice))
252	endif
253
254	ENDIF ! of IF (improved)
255
256	c Water Mass change
257	c ~~~~~~~~~~~~~~~~~
258	zq(ig,l,i_ice)=zq(ig,l,i_ice)+dzq
259	zq(ig,l,i_h2o)=zq(ig,l,i_h2o)-dzq
260
261	rice(ig,l)=max( CBRT ( (zq(ig,l,i_ice)/rho_ice
262	& +ccn(ig,l)(4./3.)pirdust(ig,l)*3.)
263	& /(ccn(ig,l)4./3.pi) ), rdust(ig,l))
264
265	enddo ! of do ig=1,ngrid
266	enddo ! of do l=1,nlay
267
268	c The following part have been commented because iceparty
269	c is not used anymore: water=>iceparty (JBM).
270
271	c else ! if not iceparty
272
273	c Saturation couche nlay a 2 :
274	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
275	c do l=nlay,2, -1
276	c do ig=1,ngrid
277	c if (zq(ig,l,i_h2o).gt.zqsat(ig,l))then
278	c zq(ig,l-1,i_h2o)= zq(ig,l-1,i_h2o)+
279	c & (zq(ig,l,i_h2o)-zqsat(ig,l))
280	c & *(pplev(ig,l)-pplev(ig,l+1))/(pplev(ig,l-1)-pplev(ig,l))
281	c zq(ig,l,i_h2o)=zqsat(ig,l)
282	c endif
283	c enddo
284	c enddo
285
286	c Saturation couche l=1 si pas iceparty
287	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
288	c do ig=1,ngridmx
289	c if (zq(ig,1,i_h2o).gt.zqsat(ig,1))then
290	c pdqscloud(ig,i_ice)=(zq(ig,1,i_h2o)-zqsat(ig,1))
291	c & (pplev(ig,1)-pplev(ig,2))/(gptimestep)
292	c zq(ig,1,i_h2o)=zqsat(ig,1)
293	c endif
294	c enddo
295
296	c endif ! of if (iceparty)
297
298	c Tendance finale
299	c ~~~~~~~~~~~~~~~
300	do l=1, nlay
301	do ig=1,ngridmx
302	pdqcloud(ig,l,i_h2o)=(zq(ig,l,i_h2o)
303	& -zq0(ig,l,i_h2o))/ptimestep
304	pdqcloud(ig,l,i_ice) =
305	& (zq(ig,l,i_ice) - zq0(ig,l,i_ice))/ptimestep
306	lw=(2834.3-0.28(zt(ig,l)-To)-0.004(zt(ig,l)-To)*2)1.e+3
307	pdtcloud(ig,l)=-pdqcloud(ig,l,i_h2o)*lw/cpp
308	end do
309	end do
310
311	c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
312	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
313	c Then that should not affect the ice particle radius
314	do ig=1,ngridmx
315	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,2)))then
316	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,3)))
317	& rice(ig,2)=rice(ig,3)
318	rice(ig,1)=rice(ig,2)
319	end if
320	end do
321
322	c**************************************************
323	c Output
324	c**************************************************
325	! NB: for diagnostics use zq(), the updated value of tracers
326
327	c do ig=1,ngridmx
328	c do l=1 ,nlay
329	c masse de glace d'eau dans la couche l
330	c icetot(ig)=icetot(ig)+masse(ig,l)*zq(ig,l,i_ice)
331	c rayon moyen des cristaux dans la colonne ig
332	c rave(ig)=rave(ig)+masse(ig,l)zq(ig,l,i_ice)rice(ig,l)
333	c enddo
334	c rave(ig)=rave(ig)/max(icetot(ig),1.e-30)
335	c if (icetot(ig)*1000.lt.0.01) rave(ig)=0.
336	c enddo ! (ngridmx)
337	c**************************************************
338
339	RETURN
340	END
341

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

Download in other formats:

Original Format