Context Navigation

watercloud.F @ 171

Last change on this file since 171 was 120, checked in by emillour, 14 years ago

-Mars GCM:

set internal computations using double precision in growthrate.F and

watercloud.F (otherwise we sometimes end up with Nans).

add extra checks in newcondens.F to avoid possibility of out of bounds

evaluation of array masse()

File size: 11.5 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	real rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)
58
59	c Outputs:
60	c -------
61
62	real pdqcloud(ngrid,nlay,nq) ! tendance de la condensation H2O(kg/kg.s-1)
63	real pdqscloud(ngrid,nq) ! flux en surface (kg.m-2.s-1)
64	REAL pdtcloud(ngrid,nlay) ! tendance temperature due
65	! a la chaleur latente
66
67	REAL rice(ngrid,nlay) ! Ice mass mean radius (m)
68	! (r_c in montmessin_2004)
69	REAL nuice(ngrid,nlay) ! Estimated effective variance
70	! of the size distribution
71
72	c local:
73	c ------
74
75	REAL CBRT
76	EXTERNAL CBRT
77	INTEGER ig,l
78
79
80	REAL zq(ngridmx,nlayermx,nqmx) ! local value of tracers
81	REAL zq0(ngridmx,nlayermx,nqmx) ! local initial value of tracers
82	REAL zqsat(ngridmx,nlayermx) ! saturation
83	REAL zt(ngridmx,nlayermx) ! local value of temperature
84
85	REAL masse (ngridmx,nlayermx)
86	REAL epaisseur (ngridmx,nlayermx)
87	! REAL rfinal ! Ice crystal radius after condensation(m)
88	REAL*8 seq ! Equilibrium saturation ration (accounting for curvature effect)
89	REAL*8 dzq ! masse de glace echangee (kg/kg)
90	REAL lw !Latent heat of sublimation (J.kg-1)
91	REAL,PARAMETER :: To=273.15 ! reference temperature, T=273.15 K
92
93	REAL Ctot
94	REAL*8 ph2o,satu
95	REAL*8 gr,Cste,up,dwn,newvap
96
97	LOGICAL,SAVE :: firstcall=.true.
98	! To use more refined microphysics, set improved to .true.
99	LOGICAL,PARAMETER :: improved=.true.
100
101	c Pour diagnostique :
102	c ~~~~~~~~~~~~~~~~~
103	c REAL icetot(ngridmx) ! Total mass of water ice (kg/m2)
104	c REAL rave(ngridmx) ! Mean crystal radius in a column (m)
105
106	! indexes of water vapour, water ice and dust tracers:
107	INTEGER,SAVE :: i_h2o=0 ! water vapour
108	INTEGER,SAVE :: i_ice=0 ! water ice
109
110	c ** un petit test de coherence
111	c --------------------------
112
113	IF (firstcall) THEN
114	IF(ngrid.NE.ngridmx) THEN
115	PRINT*,'STOP dans watercloud'
116	PRINT*,'probleme de dimensions :'
117	PRINT*,'ngrid =',ngrid
118	PRINT*,'ngridmx =',ngridmx
119	STOP
120	ENDIF
121
122	if (nq.gt.nqmx) then
123	write(,) 'stop in watercloud (nq.gt.nqmx)!'
124	write(,) 'nq=',nq,' nqmx=',nqmx
125	stop
126	endif
127
128	i_h2o=igcm_h2o_vap
129	i_ice=igcm_h2o_ice
130
131	write(,) "watercloud: i_h2o=",i_h2o
132	write(,) " i_ice=",i_ice
133
134	firstcall=.false.
135	ENDIF ! of IF (firstcall)
136
137
138	c-----------------------------------------------------------------------
139	c 1. initialisation
140	c -----------------
141
142	c On "update" la valeur de q(nqmx) (water vapor) et temperature.
143	c On effectue qqes calculs preliminaires sur les couches :
144	c masse (kg.m-2), epaisseur(m).
145
146	do l=1,nlay
147	do ig=1,ngrid
148	zq(ig,l,i_h2o)=pq(ig,l,i_h2o)+pdq(ig,l,i_h2o)*ptimestep
149	zq(ig,l,i_h2o)=max(zq(ig,l,i_h2o),1.E-30) ! FF 12/2004
150	zq0(ig,l,i_h2o)=zq(ig,l,i_h2o)
151	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
152	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
153	epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l)
154
155	zq(ig,l,i_ice)=pq(ig,l,i_ice)+pdq(ig,l,i_ice)*ptimestep
156	zq(ig,l,i_ice)=max(zq(ig,l,i_ice),0.) ! FF 12/2004
157	zq0(ig,l,i_ice)=zq(ig,l,i_ice)
158
159	c This typical profile is not used anymore; rdust is
160	c set up in updatereffrad.F.
161	c rdust(ig,l)= max(.8e-6*exp(-pzlay(ig,l)/18000.),1.e-9)
162	enddo
163	enddo
164
165	do l=1,nlay
166	do ig=1,ngrid
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	rice(ig,l)=CBRT( ( zq(ig,l,i_ice)/rho_ice+
173	& ccn(ig,l)(4./3.)pirdust(ig,l)*3. )
174	& / (ccn(ig,l)4./3.pi) )
175	rice(ig,l)=max(rice(ig,l),rdust(ig,l))
176	c Effective variance of the size distribution
177	nuice(ig,l)=nuice_ref
178	enddo ! of do ig=1,ngrid
179	enddo ! of dol=1,nlay
180
181	pdqscloud(1:ngrid,1:nq)=0
182	pdqcloud(1:ngrid,1:nlay,1:nq)=0
183	pdtcloud(1:ngrid,1:nlay)=0
184
185	c icetot(1:ngrid)=0
186	c rave(1:ngrid)=0
187
188	c ----------------------------------------------
189	c
190	c
191	c Rapport de melange a saturation dans la couche l : -------
192	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
193
194	call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
195
196	c taux de condensation (kg/kg/s-1) dans les differentes couches
197	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
198
199	c Iceparty is not used anymore: water=>iceparty (JBM).
200	c if(iceparty) then
201
202	do l=1,nlay
203	do ig=1,ngrid
204
205	IF (improved) then
206	c Improved microphysics scheme
207	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
208	Ctot = zq(ig,l,i_h2o) + zq(ig,l,i_ice)
209	ph2o = zq(ig,l,i_h2o) * 44. / 18. * pplay(ig,l)
210	satu = zq(ig,l,i_h2o) / zqsat(ig,l)
211
212	call growthrate(ptimestep,zt(ig,l),pplay(ig,l),
213	& ph2o,ph2o/satu,seq,rice(ig,l),gr)
214	Cste = ptimestep * 4. * pi * rice(ig,l)
215	* * rho_ice * ccn(ig,l)
216	up = zq(ig,l,i_h2o) + Cste * gr * seq
217	dwn = 1. + Cste * gr / zqsat(ig,l)
218	newvap = min(up/dwn,Ctot)
219
220	gr = gr * ( newvap/zqsat(ig,l) - seq )
221	dzq = min( max( Cste * gr,-zq(ig,l,i_ice) )
222	* , zq(ig,l,i_h2o) )
223
224	c Nucleation (sat ratio must be larger than a critical value)
225	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
226	if (satu.gt.1.) then
227	if (satu.le.1.4.and.zq(ig,l,i_ice).lt.1.e-8)
228	* dzq = 0.
229	endif
230
231	ELSE
232	c Old version
233	c ~~~~~~~~~~~
234	if (zq(ig,l,i_h2o).ge.zqsat(ig,l))then ! Condensation
235	dzq=zq(ig,l,i_h2o)-zqsat(ig,l)
236	elseif(zq(ig,l,i_h2o).lt.zqsat(ig,l))then ! Sublimation
237	dzq=-min(zqsat(ig,l)-zq(ig,l,i_h2o),zq(ig,l,i_ice))
238	endif
239
240	ENDIF ! of IF (improved)
241
242	c Water Mass change
243	c ~~~~~~~~~~~~~~~~~
244	zq(ig,l,i_ice)=zq(ig,l,i_ice)+dzq
245	zq(ig,l,i_h2o)=zq(ig,l,i_h2o)-dzq
246
247	rice(ig,l)=max( CBRT ( (zq(ig,l,i_ice)/rho_ice
248	& +ccn(ig,l)(4./3.)pirdust(ig,l)*3.)
249	& /(ccn(ig,l)4./3.pi) ), rdust(ig,l))
250
251	enddo ! of do ig=1,ngrid
252	enddo ! of do l=1,nlay
253
254	c The following part have been commented because iceparty
255	c is not used anymore: water=>iceparty (JBM).
256
257	c else ! if not iceparty
258
259	c Saturation couche nlay a 2 :
260	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
261	c do l=nlay,2, -1
262	c do ig=1,ngrid
263	c if (zq(ig,l,i_h2o).gt.zqsat(ig,l))then
264	c zq(ig,l-1,i_h2o)= zq(ig,l-1,i_h2o)+
265	c & (zq(ig,l,i_h2o)-zqsat(ig,l))
266	c & *(pplev(ig,l)-pplev(ig,l+1))/(pplev(ig,l-1)-pplev(ig,l))
267	c zq(ig,l,i_h2o)=zqsat(ig,l)
268	c endif
269	c enddo
270	c enddo
271
272	c Saturation couche l=1 si pas iceparty
273	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
274	c do ig=1,ngridmx
275	c if (zq(ig,1,i_h2o).gt.zqsat(ig,1))then
276	c pdqscloud(ig,i_ice)=(zq(ig,1,i_h2o)-zqsat(ig,1))
277	c & (pplev(ig,1)-pplev(ig,2))/(gptimestep)
278	c zq(ig,1,i_h2o)=zqsat(ig,1)
279	c endif
280	c enddo
281
282	c endif ! of if (iceparty)
283
284	c Tendance finale
285	c ~~~~~~~~~~~~~~~
286	do l=1, nlay
287	do ig=1,ngridmx
288	pdqcloud(ig,l,i_h2o)=(zq(ig,l,i_h2o)
289	& -zq0(ig,l,i_h2o))/ptimestep
290	pdqcloud(ig,l,i_ice) =
291	& (zq(ig,l,i_ice) - zq0(ig,l,i_ice))/ptimestep
292	lw=(2834.3-0.28(zt(ig,l)-To)-0.004(zt(ig,l)-To)*2)1.e+3
293	pdtcloud(ig,l)=-pdqcloud(ig,l,i_h2o)*lw/cpp
294	end do
295	end do
296
297	c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
298	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
299	c Then that should not affect the ice particle radius
300	do ig=1,ngridmx
301	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,2)))then
302	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,3)))
303	& rice(ig,2)=rice(ig,3)
304	rice(ig,1)=rice(ig,2)
305	end if
306	end do
307
308	c**************************************************
309	c Output
310	c**************************************************
311	! NB: for diagnostics use zq(), the updated value of tracers
312
313	c do ig=1,ngridmx
314	c do l=1 ,nlay
315	c masse de glace d'eau dans la couche l
316	c icetot(ig)=icetot(ig)+masse(ig,l)*zq(ig,l,i_ice)
317	c rayon moyen des cristaux dans la colonne ig
318	c rave(ig)=rave(ig)+masse(ig,l)zq(ig,l,i_ice)rice(ig,l)
319	c enddo
320	c rave(ig)=rave(ig)/max(icetot(ig),1.e-30)
321	c if (icetot(ig)*1000.lt.0.01) rave(ig)=0.
322	c enddo ! (ngridmx)
323	c**************************************************
324
325	RETURN
326	END
327

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

Download in other formats:

Original Format