Context Navigation

watercloud.F @ 117

Last change on this file since 117 was 83, checked in by aslmd, 14 years ago

mars + LMD_MM_MARS : modifications pour cycle de l'eau, valeurs tunees JBM

used settings reached by JBM to obtain his PhD results

alb_surfice = 0.45 --- in physiq.F and meso_physiq.F
ccn_factor = 4.5 --- in watercloud.F
nuice_sed = 0.45 --- in callsedim.F

NB: this is supposed to be further refined in the future

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,iq,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 seq ! Equilibrium saturation ration (accounting for curvature effect)
89	REAL 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 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	INTEGER i
107
108	! indexes of water vapour, water ice and dust tracers:
109	INTEGER,SAVE :: i_h2o=0 ! water vapour
110	INTEGER,SAVE :: i_ice=0 ! water ice
111	CHARACTER(LEN=20) :: tracername ! to temporarly store text
112
113	c ** un petit test de coherence
114	c --------------------------
115
116	IF (firstcall) THEN
117	IF(ngrid.NE.ngridmx) THEN
118	PRINT*,'STOP dans watercloud'
119	PRINT*,'probleme de dimensions :'
120	PRINT*,'ngrid =',ngrid
121	PRINT*,'ngridmx =',ngridmx
122	STOP
123	ENDIF
124
125	if (nq.gt.nqmx) then
126	write(,) 'stop in watercloud (nq.gt.nqmx)!'
127	write(,) 'nq=',nq,' nqmx=',nqmx
128	stop
129	endif
130
131	i_h2o=igcm_h2o_vap
132	i_ice=igcm_h2o_ice
133
134	write(,) "watercloud: i_h2o=",i_h2o
135	write(,) " i_ice=",i_ice
136
137	firstcall=.false.
138	ENDIF ! of IF (firstcall)
139
140
141	c-----------------------------------------------------------------------
142	c 1. initialisation
143	c -----------------
144
145	c On "update" la valeur de q(nqmx) (water vapor) et temperature.
146	c On effectue qqes calculs preliminaires sur les couches :
147	c masse (kg.m-2), epaisseur(m).
148
149	do l=1,nlay
150	do ig=1,ngrid
151	zq(ig,l,i_h2o)=pq(ig,l,i_h2o)+pdq(ig,l,i_h2o)*ptimestep
152	zq(ig,l,i_h2o)=max(zq(ig,l,i_h2o),1.E-30) ! FF 12/2004
153	zq0(ig,l,i_h2o)=zq(ig,l,i_h2o)
154	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
155	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
156	epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l)
157
158	zq(ig,l,i_ice)=pq(ig,l,i_ice)+pdq(ig,l,i_ice)*ptimestep
159	zq(ig,l,i_ice)=max(zq(ig,l,i_ice),0.) ! FF 12/2004
160	zq0(ig,l,i_ice)=zq(ig,l,i_ice)
161
162	c This typical profile is not used anymore; rdust is
163	c set up in updatereffrad.F.
164	c rdust(ig,l)= max(.8e-6*exp(-pzlay(ig,l)/18000.),1.e-9)
165	enddo
166	enddo
167
168	do l=1,nlay
169	do ig=1,ngrid
170	c Calcul du rayon moyen des particules de glace.
171	c Hypothese : Dans une couche, la glace presente se
172	c repartit uniformement autour du nbre de poussieres
173	c dont le rayon moyen est prescrit par rdust.
174	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
175	rice(ig,l)=CBRT( ( zq(ig,l,i_ice)/rho_ice+
176	& ccn(ig,l)(4./3.)pirdust(ig,l)*3. )
177	& / (ccn(ig,l)4./3.pi) )
178	rice(ig,l)=max(rice(ig,l),rdust(ig,l))
179	c Effective variance of the size distribution
180	nuice(ig,l)=nuice_ref
181	enddo ! of do ig=1,ngrid
182	enddo ! of dol=1,nlay
183
184	pdqscloud(1:ngrid,1:nq)=0
185	pdqcloud(1:ngrid,1:nlay,1:nq)=0
186	pdtcloud(1:ngrid,1:nlay)=0
187
188	c icetot(1:ngrid)=0
189	c rave(1:ngrid)=0
190
191	c ----------------------------------------------
192	c
193	c
194	c Rapport de melange a saturation dans la couche l : -------
195	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
196
197	call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
198
199	c taux de condensation (kg/kg/s-1) dans les differentes couches
200	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
201
202	c Iceparty is not used anymore: water=>iceparty (JBM).
203	c if(iceparty) then
204
205	do l=1,nlay
206	do ig=1,ngrid
207
208	IF (improved) then
209	c Improved microphysics scheme
210	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
211
212	Ctot = zq(ig,l,i_h2o) + zq(ig,l,i_ice)
213	ph2o = zq(ig,l,i_h2o) * 44. / 18. * pplay(ig,l)
214	satu = zq(ig,l,i_h2o) / zqsat(ig,l)
215
216	call growthrate(ptimestep,zt(ig,l),pplay(ig,l),
217	& ph2o,ph2o/satu,seq,rice(ig,l),gr)
218	Cste = ptimestep * 4. * pi * rice(ig,l)
219	* * rho_ice * ccn(ig,l)
220	up = zq(ig,l,i_h2o) + Cste * gr * seq
221	dwn = 1. + Cste * gr / zqsat(ig,l)
222	newvap = min(up/dwn,Ctot)
223
224	gr = gr * ( newvap/zqsat(ig,l) - seq )
225	dzq = min( max( Cste * gr,-zq(ig,l,i_ice) )
226	* , zq(ig,l,i_h2o) )
227
228	c Nucleation (sat ratio must be larger than a critical value)
229	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
230	if (satu.gt.1.) then
231	if (satu.le.1.4.and.zq(ig,l,i_ice).lt.1.e-8)
232	* dzq = 0.
233	endif
234
235	ELSE
236	c Old version
237	c ~~~~~~~~~~~
238	if (zq(ig,l,i_h2o).ge.zqsat(ig,l))then ! Condensation
239	dzq=zq(ig,l,i_h2o)-zqsat(ig,l)
240	elseif(zq(ig,l,i_h2o).lt.zqsat(ig,l))then ! Sublimation
241	dzq=-min(zqsat(ig,l)-zq(ig,l,i_h2o),zq(ig,l,i_ice))
242	endif
243
244	ENDIF ! of IF (improved)
245
246	c Water Mass change
247	c ~~~~~~~~~~~~~~~~~
248	zq(ig,l,i_ice)=zq(ig,l,i_ice)+dzq
249	zq(ig,l,i_h2o)=zq(ig,l,i_h2o)-dzq
250
251	rice(ig,l)=max( CBRT ( (zq(ig,l,i_ice)/rho_ice
252	& +ccn(ig,l)(4./3.)pirdust(ig,l)*3.)
253	& /(ccn(ig,l)4./3.pi) ), rdust(ig,l))
254
255	enddo ! of do ig=1,ngrid
256	enddo ! of do l=1,nlay
257
258	c The following part have been commented because iceparty
259	c is not used anymore: water=>iceparty (JBM).
260
261	c else ! if not iceparty
262
263	c Saturation couche nlay a 2 :
264	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
265	c do l=nlay,2, -1
266	c do ig=1,ngrid
267	c if (zq(ig,l,i_h2o).gt.zqsat(ig,l))then
268	c zq(ig,l-1,i_h2o)= zq(ig,l-1,i_h2o)+
269	c & (zq(ig,l,i_h2o)-zqsat(ig,l))
270	c & *(pplev(ig,l)-pplev(ig,l+1))/(pplev(ig,l-1)-pplev(ig,l))
271	c zq(ig,l,i_h2o)=zqsat(ig,l)
272	c endif
273	c enddo
274	c enddo
275
276	c Saturation couche l=1 si pas iceparty
277	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
278	c do ig=1,ngridmx
279	c if (zq(ig,1,i_h2o).gt.zqsat(ig,1))then
280	c pdqscloud(ig,i_ice)=(zq(ig,1,i_h2o)-zqsat(ig,1))
281	c & (pplev(ig,1)-pplev(ig,2))/(gptimestep)
282	c zq(ig,1,i_h2o)=zqsat(ig,1)
283	c endif
284	c enddo
285
286	c endif ! of if (iceparty)
287
288	c Tendance finale
289	c ~~~~~~~~~~~~~~~
290	do l=1, nlay
291	do ig=1,ngridmx
292	pdqcloud(ig,l,i_h2o)=(zq(ig,l,i_h2o)
293	& -zq0(ig,l,i_h2o))/ptimestep
294	pdqcloud(ig,l,i_ice) =
295	& (zq(ig,l,i_ice) - zq0(ig,l,i_ice))/ptimestep
296	lw=(2834.3-0.28(zt(ig,l)-To)-0.004(zt(ig,l)-To)*2)1.e+3
297	pdtcloud(ig,l)=-pdqcloud(ig,l,i_h2o)*lw/cpp
298	end do
299	end do
300
301	c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
302	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
303	c Then that should not affect the ice particle radius
304	do ig=1,ngridmx
305	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,2)))then
306	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,3)))
307	& rice(ig,2)=rice(ig,3)
308	rice(ig,1)=rice(ig,2)
309	end if
310	end do
311
312	c**************************************************
313	c Output
314	c**************************************************
315	! NB: for diagnostics use zq(), the updated value of tracers
316
317	c do ig=1,ngridmx
318	c do l=1 ,nlay
319	c masse de glace d'eau dans la couche l
320	c icetot(ig)=icetot(ig)+masse(ig,l)*zq(ig,l,i_ice)
321	c rayon moyen des cristaux dans la colonne ig
322	c rave(ig)=rave(ig)+masse(ig,l)zq(ig,l,i_ice)rice(ig,l)
323	c enddo
324	c rave(ig)=rave(ig)/max(icetot(ig),1.e-30)
325	c if (icetot(ig)*1000.lt.0.01) rave(ig)=0.
326	c enddo ! (ngridmx)
327	c**************************************************
328
329	RETURN
330	END
331

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

Download in other formats:

Original Format