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 Francois Forget : change nuclei density & outputs |
---|
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 | |
---|
35 | #include "fisice.h" |
---|
36 | #include "comgeomfi.h" |
---|
37 | c |
---|
38 | c arguments: |
---|
39 | c ---------- |
---|
40 | |
---|
41 | INTEGER ngrid,nlay,nsize |
---|
42 | INTEGER icount |
---|
43 | REAL ptimestep ! pas de temps physique (s) |
---|
44 | REAL pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa) |
---|
45 | REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa) |
---|
46 | REAL pzlev(ngrid,nlay+1) ! altitude at layer boundaries |
---|
47 | REAL pzlay(ngrid,nlay) ! altitude at the middle of the layers |
---|
48 | REAL pdpsrf(ngrid) ! tendance surf pressure |
---|
49 | REAL pt(ngrid,nlay) ! temperature au centre des couches (K) |
---|
50 | REAL pdt(ngrid,nlay) ! tendance temperature des autres param. |
---|
51 | REAL pdtcloud(ngrid,nlay) ! tendance temperature due a la chaleur latente |
---|
52 | REAL tau(ngridmx,nsize) |
---|
53 | REAL zls |
---|
54 | |
---|
55 | c Traceurs : |
---|
56 | integer nq ! nombre de traceurs |
---|
57 | real pq(ngrid,nlay,nq) ! traceur (kg/kg) |
---|
58 | real pdq(ngrid,nlay,nq) ! tendance avant condensation (kg/kg.s-1) |
---|
59 | real pdqcloud(ngrid,nlay,nq) ! tendance de la condesation H2O(kg/kg.s-1) |
---|
60 | real qsurf(ngrid,nq) |
---|
61 | real pdqscloud(ngrid,nq) ! flux en surface (kg.m-2.s-1) |
---|
62 | |
---|
63 | c local: |
---|
64 | c ------ |
---|
65 | |
---|
66 | REAL CBRT |
---|
67 | EXTERNAL CBRT |
---|
68 | INTEGER ig,l |
---|
69 | |
---|
70 | LOGICAL firstcall,improved |
---|
71 | SAVE firstcall |
---|
72 | |
---|
73 | REAL zq(ngridmx,nlayermx,nqmx) |
---|
74 | REAL zq0(ngridmx,nlayermx,nqmx) |
---|
75 | REAL zqsat(ngridmx,nlayermx) |
---|
76 | REAL zt(ngridmx,nlayermx) |
---|
77 | |
---|
78 | |
---|
79 | REAL masse (ngridmx,nlayermx) |
---|
80 | REAL epaisseur (ngridmx,nlayermx) |
---|
81 | REAL dustcores(ngridmx,nlayermx) !Dust number density (#/kg) |
---|
82 | REAL rfinal ! Ice crystal radius after condensation(m) |
---|
83 | REAL seq ! Equilibrium saturation ration (accounting for curvature effect) |
---|
84 | REAL dzq ! masse de glace echangee (kg/kg) |
---|
85 | REAL lw !Latent heat of sublimation (J.kg-1) |
---|
86 | REAL To |
---|
87 | |
---|
88 | REAL Ctot |
---|
89 | REAL*8 ph2o,satu |
---|
90 | REAL gr,Cste,up,dwn,newvap |
---|
91 | |
---|
92 | DATA firstcall/.true./ |
---|
93 | DATA improved/.false./ ! Actionne une microphysique plus raffinee |
---|
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 | REAL mtot(ngridmx) ! Total mass of water vapor (kg/m2) |
---|
104 | REAL icetot(ngridmx) ! Total mass of water ice (kg/m2) |
---|
105 | REAL rave(ngridmx) ! Mean crystal radius in a column (m) |
---|
106 | REAL cloudmass(ngridmx,nlayermx) ! mass of ice in each layer (kg/m2) |
---|
107 | REAL op825(ngridmx,nlayermx) ! abs optical depth at 825 cm-1 |
---|
108 | REAL tau825(ngridmx) ! abs optical depth at 825 cm-1 |
---|
109 | REAL waterhem(2,3) |
---|
110 | real a,b, Qabs |
---|
111 | |
---|
112 | REAL icetot2pm(ngridmx) |
---|
113 | REAL mtot2pm(ngridmx) |
---|
114 | REAL rave2pm(ngridmx) |
---|
115 | REAL qsurf2pm(ngridmx) |
---|
116 | INTEGER i,nit |
---|
117 | REAL ecart,loctime |
---|
118 | integer univtime,lon2pm(ngridmx) |
---|
119 | |
---|
120 | integer ig_vl2 |
---|
121 | real latvl2,lonvl2 |
---|
122 | |
---|
123 | logical tesfield ! output of TES like data |
---|
124 | data tesfield /.false./ |
---|
125 | |
---|
126 | c ** un petit test de coherence |
---|
127 | c -------------------------- |
---|
128 | |
---|
129 | IF (firstcall) THEN |
---|
130 | IF(ngrid.NE.ngridmx) THEN |
---|
131 | PRINT*,'STOP dans watercloud' |
---|
132 | PRINT*,'probleme de dimensions :' |
---|
133 | PRINT*,'ngrid =',ngrid |
---|
134 | PRINT*,'ngridmx =',ngridmx |
---|
135 | STOP |
---|
136 | ENDIF |
---|
137 | if(nq.gt.nqmx) then |
---|
138 | write(*,*) 'stop in watercloud (nq.gt.nqmx)!' |
---|
139 | write(*,*) 'nq=',nq,' nqmx=',nqmx |
---|
140 | stop |
---|
141 | endif |
---|
142 | firstcall=.false. |
---|
143 | ENDIF |
---|
144 | |
---|
145 | |
---|
146 | c----------------------------------------------------------------------- |
---|
147 | c 1. initialisation |
---|
148 | c ----------------- |
---|
149 | |
---|
150 | c On "update" la valeur de q(nqmx) (water vapor) et temperature. |
---|
151 | c On effectue qqes calculs preliminaires sur les couches : |
---|
152 | c masse (kg.m-2), epaisseur(m). |
---|
153 | |
---|
154 | do l=1,nlay |
---|
155 | do ig=1,ngrid |
---|
156 | zq(ig,l,nq)=pq(ig,l,nq)+pdq(ig,l,nq)*ptimestep |
---|
157 | zq(ig,l,nq)=max(zq(ig,l,nq),1.E-30) ! FF 12/2004 |
---|
158 | zq0(ig,l,nq)=zq(ig,l,nq) |
---|
159 | zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep |
---|
160 | masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g |
---|
161 | epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l) |
---|
162 | if(iceparty) then |
---|
163 | zq(ig,l,nq-1)=pq(ig,l,nq-1)+pdq(ig,l,nq-1)*ptimestep |
---|
164 | zq(ig,l,nq-1)=max(zq(ig,l,nq-1),0.) ! FF 12/2004 |
---|
165 | zq0(ig,l,nq-1)=zq(ig,l,nq-1) |
---|
166 | rdust(ig,l)= max(.8e-6*exp(-pzlay(ig,l)/18000.),1.e-9) |
---|
167 | c Calcul du rayon moyen des particules de glace. |
---|
168 | c Hypothese : Dans une couche, la glace presente se |
---|
169 | c repartit uniformement autour du nbre de poussieres |
---|
170 | c dont le rayon moyen est prescrit par rdust. |
---|
171 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
172 | if ((dustbin.eq.1).and..not.(doubleq)) then |
---|
173 | zq(ig,l,1)=pq(ig,l,1)+pdq(ig,l,1)*ptimestep |
---|
174 | dustcores(ig,l)=max(zq(ig,l,1)/ |
---|
175 | & (rho_dust*4./3.*pi*radius(1)**3.),1.e-9) |
---|
176 | else |
---|
177 | dustcores(ig,l)=( epaisseur(ig,l)/masse(ig,l) ) * |
---|
178 | & 2.e+6/0.1*max(tau(ig,1),0.001)*exp(-pzlay(ig,l)/10000.) |
---|
179 | |
---|
180 | c TEMPORAIRE : réduction du nombre de nuclei FF 04/2008 : |
---|
181 | c dustcores(ig,l) = dustcores(ig,l) / 27. ! reduction facteur 3 |
---|
182 | dustcores(ig,l) = dustcores(ig,l) / 8. ! reduction facteur 2 |
---|
183 | endif |
---|
184 | rice(ig,l)=CBRT( ( zq(ig,l,nq-1)/rho_ice+ |
---|
185 | & dustcores(ig,l)*(4./3.)*pi*rdust(ig,l)**3. ) |
---|
186 | & / (dustcores(ig,l)*4./3.*pi) ) |
---|
187 | rice(ig,l)=max(rice(ig,l),rdust(ig,l)) |
---|
188 | endif |
---|
189 | enddo |
---|
190 | enddo |
---|
191 | |
---|
192 | call zerophys(ngrid*nq,pdqscloud) |
---|
193 | call zerophys(ngrid*nlay*nq,pdqcloud) |
---|
194 | call zerophys(ngrid*nlay,pdtcloud) |
---|
195 | call zerophys(ngrid,mtot) |
---|
196 | |
---|
197 | call zerophys(ngrid,icetot) |
---|
198 | call zerophys(ngrid,rave) |
---|
199 | call zerophys(ngrid,cloudmass) |
---|
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 | c (Pour diagnostic seulement !) |
---|
212 | if(.not.iceparty)then |
---|
213 | do l=1, nlay |
---|
214 | do ig=1,ngridmx |
---|
215 | taucond(ig,l)= max((zq(ig,l,nq)-zqsat(ig,l))/ptimestep,0.) |
---|
216 | end do |
---|
217 | end do |
---|
218 | endif |
---|
219 | |
---|
220 | |
---|
221 | if(iceparty) then |
---|
222 | do l=1,nlay |
---|
223 | do ig=1,ngrid |
---|
224 | |
---|
225 | If (improved) then |
---|
226 | c Improved microphysics scheme |
---|
227 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
228 | |
---|
229 | Ctot = zq(ig,l,nq) + zq(ig,l,nq-1) |
---|
230 | ph2o = zq(ig,l,nq) * 44. / 18. * pplay(ig,l) |
---|
231 | satu = zq(ig,l,nq) / zqsat(ig,l) |
---|
232 | |
---|
233 | call growthrate(ptimestep,zt(ig,l),pplay(ig,l), |
---|
234 | & ph2o,ph2o/satu,seq,rice(ig,l),gr) |
---|
235 | Cste = ptimestep * 4. * pi * rice(ig,l) |
---|
236 | * * rho_ice * dustcores(ig,l) |
---|
237 | up = zq(ig,l,nq) + Cste * gr * seq |
---|
238 | dwn = 1. + Cste * gr / zqsat(ig,l) |
---|
239 | newvap = min(up/dwn,Ctot) |
---|
240 | |
---|
241 | gr = gr * ( newvap/zqsat(ig,l) - seq ) |
---|
242 | dzq = min( max( Cste * gr,-zq(ig,l,nq-1) ) |
---|
243 | * , zq(ig,l,nq) ) |
---|
244 | |
---|
245 | c Nucleation (sat ratio must be larger than a critical value) |
---|
246 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
247 | if (satu.gt.1.) then |
---|
248 | if (satu.le.1.4.and.zq(ig,l,nq-1).lt.1.e-8) |
---|
249 | * dzq = 0. |
---|
250 | endif |
---|
251 | |
---|
252 | Else |
---|
253 | c Old version |
---|
254 | c ~~~~~~~~~~~ |
---|
255 | if (zq(ig,l,nq).ge.zqsat(ig,l))then ! Condensation |
---|
256 | dzq=zq(ig,l,nq)-zqsat(ig,l) |
---|
257 | elseif(zq(ig,l,nq).lt.zqsat(ig,l))then ! Sublimation |
---|
258 | dzq=-min(zqsat(ig,l)-zq(ig,l,nq),zq(ig,l,nq-1)) |
---|
259 | endif |
---|
260 | |
---|
261 | Endif |
---|
262 | |
---|
263 | c Water Mass change |
---|
264 | c ~~~~~~~~~~~~~~~~~ |
---|
265 | zq(ig,l,nq-1)=zq(ig,l,nq-1)+dzq |
---|
266 | zq(ig,l,nq)=zq(ig,l,nq)-dzq |
---|
267 | |
---|
268 | rice(ig,l)=max( CBRT ( (zq(ig,l,nq-1)/rho_ice |
---|
269 | & +dustcores(ig,l)*(4./3.)*pi*rdust(ig,l)**3.) |
---|
270 | & /(dustcores(ig,l)*4./3.*pi) ), rdust(ig,l)) |
---|
271 | |
---|
272 | if(activice.and.pclc(ig).gt.0) |
---|
273 | & rice(ig,l)=rice(ig,l)/CBRT(pclc(ig)) |
---|
274 | |
---|
275 | enddo |
---|
276 | enddo |
---|
277 | |
---|
278 | else ! if not iceparty |
---|
279 | |
---|
280 | c Saturation couche nlay a 2 : |
---|
281 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
282 | do l=nlay,2, -1 |
---|
283 | do ig=1,ngrid |
---|
284 | if (zq(ig,l,nq).gt.zqsat(ig,l))then |
---|
285 | zq(ig,l-1,nq) = zq(ig,l-1,nq)+(zq(ig,l,nq)-zqsat(ig,l)) |
---|
286 | & *(pplev(ig,l)-pplev(ig,l+1))/(pplev(ig,l-1)-pplev(ig,l)) |
---|
287 | zq(ig,l,nq)=zqsat(ig,l) |
---|
288 | endif |
---|
289 | enddo |
---|
290 | enddo |
---|
291 | |
---|
292 | c Saturation couche l=1 si pas iceparty |
---|
293 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
294 | do ig=1,ngridmx |
---|
295 | if (zq(ig,1,nq).gt.zqsat(ig,1))then |
---|
296 | pdqscloud(ig,nq)= (zq(ig,1,nq)-zqsat(ig,1)) |
---|
297 | & *(pplev(ig,1)-pplev(ig,2))/(g*ptimestep) |
---|
298 | zq(ig,1,nq)=zqsat(ig,1) |
---|
299 | endif |
---|
300 | end do |
---|
301 | |
---|
302 | endif ! (iceparty) |
---|
303 | |
---|
304 | c Tendance finale |
---|
305 | c ~~~~~~~~~~~~~~~ |
---|
306 | do l=1, nlay |
---|
307 | do ig=1,ngridmx |
---|
308 | pdqcloud(ig,l,nq) = (zq(ig,l,nq) - zq0(ig,l,nq))/ptimestep |
---|
309 | if(iceparty) then |
---|
310 | pdqcloud(ig,l,nq-1) = |
---|
311 | & (zq(ig,l,nq-1) - zq0(ig,l,nq-1))/ptimestep |
---|
312 | endif |
---|
313 | lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)**2)*1.e+3 |
---|
314 | pdtcloud(ig,l)=-pdqcloud(ig,l,nq)*lw/cpp |
---|
315 | end do |
---|
316 | end do |
---|
317 | |
---|
318 | c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005 |
---|
319 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
320 | c Then that should not affect the ice particle radius |
---|
321 | do ig=1,ngridmx |
---|
322 | if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then |
---|
323 | if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3))) |
---|
324 | & rice(ig,2)=rice(ig,3) |
---|
325 | rice(ig,1)=rice(ig,2) |
---|
326 | end if |
---|
327 | end do |
---|
328 | |
---|
329 | c************************************************** |
---|
330 | c Output |
---|
331 | c************************************************** |
---|
332 | do ig=1,ngridmx |
---|
333 | do l=1 ,nlay |
---|
334 | c masse de vapeur d'eau dans la couche l |
---|
335 | mtot(ig)=mtot(ig)+pq(ig,l,nq)*masse(ig,l) |
---|
336 | if(iceparty) then |
---|
337 | c masse de glace d'eau dans la couche l |
---|
338 | icetot(ig)=icetot(ig)+masse(ig,l)*pq(ig,l,nq-1) |
---|
339 | c rayon moyen des cristaux dans la colonne ig |
---|
340 | rave(ig)=rave(ig)+masse(ig,l)*pq(ig,l,nq-1)*rice(ig,l) |
---|
341 | cloudmass(ig,l)=masse(ig,l)*pq(ig,l,nq-1) |
---|
342 | endif |
---|
343 | enddo |
---|
344 | if(iceparty) then |
---|
345 | rave(ig)=rave(ig)/max(icetot(ig),1.e-30) |
---|
346 | if (icetot(ig)*1000.lt.0.01) rave(ig)=0. |
---|
347 | endif |
---|
348 | enddo ! (ngridmx) |
---|
349 | |
---|
350 | c Computing abs optical depth at 825 cm-1 in each layer |
---|
351 | c to simulate NEW TES retrieval |
---|
352 | do ig=1,ngridmx |
---|
353 | tau825(ig)=0. |
---|
354 | do l=1 ,nlay |
---|
355 | Qabs=min(max(0.4e6*rice(ig,l)-0.05 ,0.),1.2) |
---|
356 | op825(ig,l)= 0.75*Qabs*pq(ig,l,nq-1)*masse(ig,l) |
---|
357 | & / (rho_ice*rice(ig,l)) |
---|
358 | tau825(ig)=tau825(ig)+ op825(ig,l) |
---|
359 | end do |
---|
360 | end do |
---|
361 | |
---|
362 | c ****************************************** |
---|
363 | c Output in diagfi |
---|
364 | c ****************************************** |
---|
365 | |
---|
366 | if (ngrid.gt.1) then |
---|
367 | |
---|
368 | CALL WRITEDIAGFI(ngridmx,'mtot', |
---|
369 | & 'total mass of water vapor','kg/m2',2,mtot) |
---|
370 | |
---|
371 | if (callstats) then |
---|
372 | call wstats(ngrid,"mtot", |
---|
373 | . "total mass of water vapor","kg/m2",2,mtot) |
---|
374 | endif |
---|
375 | |
---|
376 | c if(.not.iceparty) |
---|
377 | c & call WRITEDIAGFI(ngridmx,'taucond', |
---|
378 | c & 'taux cond H2O ice','kg/kg/s',3,taucond) |
---|
379 | c |
---|
380 | if(iceparty) then |
---|
381 | CALL WRITEDIAGFI(ngridmx,'icetot', |
---|
382 | & 'total mass of water ice','kg/m2',2,icetot) |
---|
383 | if (callstats) then |
---|
384 | call wstats(ngrid,"icetot", |
---|
385 | . "total mass of water ice","kg/m2",2,icetot) |
---|
386 | endif |
---|
387 | c CALL WRITEDIAGFI(ngridmx,'cloudmass', |
---|
388 | c & 'mass of water ice/layer','kg/m2',3,cloudmass) |
---|
389 | c CALL WRITEDIAGFI(ngridmx,'rice','ice radius', |
---|
390 | c & 'meter',3,rice) |
---|
391 | CALL WRITEDIAGFI(ngridmx,'rave','Mean ice radius', |
---|
392 | & 'meter',2,rave) |
---|
393 | CALL WRITEDIAGFI(ngridmx,'tauice','tau abs 825 cm-1', |
---|
394 | & '',2,tau825) |
---|
395 | endif |
---|
396 | |
---|
397 | else |
---|
398 | |
---|
399 | c CALL WRITEG1D(ngridmx,1,mtot,'gas','kg/m2') |
---|
400 | c CALL WRITEG1D(ngridmx,1,icetot,'ice','kg/m2') |
---|
401 | c call WRITEG1D(ngridmx,nlay,rice,'rice','um') |
---|
402 | c call WRITEG1D(ngridmx,nlay,dustcores,'pouss','#/layer') |
---|
403 | |
---|
404 | endif |
---|
405 | 99 continue |
---|
406 | |
---|
407 | RETURN |
---|
408 | END |
---|
409 | |
---|