Context Navigation

improvedclouds.F @ 358

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

LMDZ.MARS:

A FIRST VERSION WITH SCAVENGING
TRANSPARENT TO CASUAL USER
FOR WATER CYCLE FOLKS, SEE BELOW

[TN and AS solved conflicts due to r330 >> r357]

07/11/11 == JBM

Changed watercloud.F to call two separate routines,

simpleclouds.F or improvedclouds.F, which are a simplified and
full microphysical scheme for cloud formation, respectively.
Removed the tag called "improved" in watercloud.F, and added
another tag called "microphys" which is defined in callphys.F
instead. Changed routines: callkeys, inifis, physiq, watercloud.

Reimplemented the use of typical profiles of dust particle sizes

and CCNs in simpleclouds.F. Corrected the previously used
analytical CCN profile. Moved ccn_factor to simpleclouds.F,
which won't be used in the new microphysical scheme. Changed
routines: aeropacity, physiq, simpleclouds, watercloud.

Computed rdust at the end of callsedim instead of updatereffrad,

except at firstcall. Renamed rfall into rsedcloud and moved it
in simpleclouds. Moved nuice_sed to initracer.F and added it to
tracer.h. Changed routines: callsedim, physiq, tracer.h,
watercloud, initracer, simpleclouds, updatereffrad.

Added two tracers for the CCNs. Added the different tests in

initracer.F to make sure that, for example, microphys is not
used without doubleq, etc. Corrected an inconsistency in
callsedim.F, and changed the way r0 is updated. Changes
routines: callsedim, inifis, initracer, physiq, testphys1d,
tracer.h.

Added the ability to have a spatially variable density in

newsedim (same method as that used for the radius of
sedimentation). Required the addition of one input to newsedim,
which is the size of the density variable "rho". Changed
routines: callsedim, newsedim.

Added an output to aeropacity called "tauscaling", which is a

factor to convert dust_mass and dust_number to true quantities,
based on tauvis. Removed ccn and qdust from aeropacity, which
are obsolete now.

Wrote improvedcloud.F which includes all the microphysical

scheme, and connected it to the sedimentation scheme. Added and
changed routines: callsedim, physiq, growthrate, nuclea,
improvedclouds, initracer, watercloud, microphys.h.

07/11/11 == TN

Added CCN & dust tracers updates in physiq.F

Corrected a bug that can give negative CCN numbers, due to the
use of single precision values (only 6 decimals) whereas up to 10E+10
CCN can be formed or disappears...

Corrected physical bug that causes h2o_vap or h2o_ice

to be negative in improvedclouds.F.

Corrected physical bug that causes CCN & dust tracers

to be negative in improvedclouds.F when all ice sublimates and
releases dust

Added parameter contact mteta in callphys.def

Default value is still 0.95, see inifis.F

Changed tendancy computation for dust_number in improvedclouds.F

that was not the right one. Indeed, the scavenged dust_number tracer
is computed from the dust_mass one, and its tendancy before scavenging
must be taken into account to compute its scavenging's tendancy.

File size: 22.5 KB

Line
1	subroutine improvedclouds(ngrid,nlay,ptimestep,
2	& pplay,pt,pdt,
3	& pq,pdq,pdqcloud,pdqscloud,pdtcloud,
4	& nq,tauscaling,rdust,rice,nuice,
5	& rsedcloud,rhocloud)
6	implicit none
7	c------------------------------------------------------------------
8	c This routine is used to form clouds when a parcel of the GCM is
9	c saturated. It includes the ability to have supersaturation, a
10	c computation of the nucleation rates, growthrates and the
11	c scavenging of dust particles by clouds.
12	c It is worth noting that the amount of dust is computed using the
13	c dust optical depth computed in aeropacity.F. That's why
14	c the variable called "tauscaling" is used to convert
15	c pq(dust_mass) and pq(dust_number), which are relative
16	c quantities, to absolute and realistic quantities stored in zq.
17	c This has to be done to convert the inputs into absolute
18	c values, but also to convert the outputs back into relative
19	c values which are then used by the sedimentation and advection
20	c schemes.
21
22	c Authors: J.-B. Madeleine, based on the work by Franck Montmessin
23	c (October 2011)
24	c------------------------------------------------------------------
25	#include "dimensions.h"
26	#include "dimphys.h"
27	#include "comcstfi.h"
28	#include "callkeys.h"
29	#include "tracer.h"
30	#include "comgeomfi.h"
31	#include "dimradmars.h"
32	#include "microphys.h"
33	c------------------------------------------------------------------
34	c Inputs:
35
36	INTEGER ngrid,nlay
37	integer nq ! nombre de traceurs
38	REAL ptimestep ! pas de temps physique (s)
39	REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
40	REAL pt(ngrid,nlay) ! temperature at the middle of the
41	! layers (K)
42	REAL pdt(ngrid,nlay) ! tendance temperature des autres
43	! param.
44	REAL pq(ngrid,nlay,nq) ! traceur (kg/kg)
45	REAL pdq(ngrid,nlay,nq) ! tendance avant condensation
46	! (kg/kg.s-1)
47	REAL tauscaling(ngridmx) ! Convertion factor for qdust and Ndust
48	REAL rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)
49
50	c Outputs:
51	REAL rice(ngrid,nlay) ! Ice mass mean radius (m)
52	! (r_c in montmessin_2004)
53	REAL nuice(ngrid,nlay) ! Estimated effective variance
54	! of the size distribution
55	REAL rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius
56	REAL rhocloud(ngridmx,nlayermx) ! Cloud density (kg.m-3)
57	REAL pdqcloud(ngrid,nlay,nq) ! tendance de la condensation
58	! H2O(kg/kg.s-1)
59	REAL pdqscloud(ngrid,nq) ! flux en surface (kg.m-2.s-1)
60	REAL pdtcloud(ngrid,nlay) ! tendance temperature due
61	! a la chaleur latente
62
63	c------------------------------------------------------------------
64	c Local variables:
65
66	LOGICAL firstcall
67	DATA firstcall/.true./
68	SAVE firstcall
69
70	REAL*8 derf ! Error function
71	!external derf
72
73	REAL CBRT
74	EXTERNAL CBRT
75
76	INTEGER ig,l,i
77
78	REAL zq(ngridmx,nlayermx,nqmx) ! local value of tracers
79	REAL zq0(ngridmx,nlayermx,nqmx) ! local initial value of tracers
80	REAL zt(ngridmx,nlayermx) ! local value of temperature
81	REAL zqsat(ngridmx,nlayermx) ! saturation
82	REAL lw !Latent heat of sublimation (J.kg-1)
83	REAL Cste
84	REAL dMice ! mass of condensated ice
85	REAL sumcheck
86	REAL*8 ph2o ! Water vapor partial pressure (Pa)
87	REAL*8 satu ! Water vapor saturation ratio over ice
88	REAL*8 Mo,No
89	REAL*8 dN,dM,newvap
90	REAL*8 Rn, Rm, dev2
91	REAL*8 n_aer(nbin_cld) ! number conc. of particle/each size bin
92	REAL*8 m_aer(nbin_cld) ! mass mixing ratio of particle/each size bin
93	REAL*8 rate(nbin_cld) ! nucleation rate
94	REAL*8 up,dwn,Ctot,gr,seq
95	REAL*8 sig ! Water-ice/air surface tension (N.m)
96	EXTERNAL sig
97
98	c Parameters of the size discretization
99	c used by the microphysical scheme
100	DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m)
101	DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m)
102	DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.0001e-6
103	! Minimum boundary radius (m)
104	DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m)
105	DOUBLE PRECISION vrat_cld ! Volume ratio
106	DOUBLE PRECISION rb_cld(nbin_cld+1)! boundary values of each rad_cld bin (m)
107	SAVE rb_cld
108	DOUBLE PRECISION dr_cld(nbin_cld)! width of each rad_cld bin (m)
109	DOUBLE PRECISION vol_cld(nbin_cld) ! particle volume for each bin (m3)
110
111	REAL sigma_ice ! Variance of the ice and CCN distributions
112	SAVE sigma_ice
113
114	c some outputs for 1D
115	REAL satu_out(ngridmx,nlayermx) ! satu ratio for output
116	REAL dN_out(ngridmx,nlayermx) ! mass variation for output
117	REAL dM_out(ngridmx,nlayermx) ! number variation for output
118	REAL Mcon_out(ngridmx,nlayermx) ! mass to be condensed (not dMice !!)
119	REAL gr_out(ngridmx,nlayermx) ! for 1d output
120	REAL newvap_out(ngridmx,nlayermx) ! for 1d output
121
122	c------------------------------------------------------------------
123
124	IF (firstcall) THEN
125
126	c Definition of the size grid
127	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~
128	c rad_cld is the primary radius grid used for microphysics computation.
129	c The grid spacing is computed assuming a constant volume ratio
130	c between two consecutive bins; i.e. vrat_cld.
131	c vrat_cld is determined from the boundary values of the size grid:
132	c rmin_cld and rmax_cld.
133	c The rb_cld array contains the boundary values of each rad_cld bin.
134	c dr_cld is the width of each rad_cld bin.
135
136	c Volume ratio between two adjacent bins
137	vrat_cld = dlog(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
138	vrat_cld = dexp(vrat_cld)
139	write(,) "vrat_cld", vrat_cld
140
141	rb_cld(1) = rbmin_cld
142	rad_cld(1) = rmin_cld
143	vol_cld(1) = 4./3. * dble(pi) * rmin_cld**3.
144
145	do i=1,nbin_cld-1
146	rad_cld(i+1) = rad_cld(i) * vrat_cld**(1./3.)
147	vol_cld(i+1) = vol_cld(i) * vrat_cld
148	enddo
149
150	do i=1,nbin_cld
151	rb_cld(i+1)= ( (2.vrat_cld) / (vrat_cld+1.) )(1./3.)
152	& rad_cld(i)
153	dr_cld(i) = rb_cld(i+1) - rb_cld(i)
154	enddo
155	rb_cld(nbin_cld+1) = rbmax_cld
156	dr_cld(nbin_cld) = rb_cld(nbin_cld+1) - rb_cld(nbin_cld)
157
158	print*, ' '
159	print*,'Microphysics: size bin information:'
160	print*,'i,rb_cld(i), rad_cld(i),dr_cld(i)'
161	print*,'-----------------------------------'
162	do i=1,nbin_cld
163	write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i),
164	& dr_cld(i)
165	enddo
166	write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1)
167	print*,'-----------------------------------'
168
169	c Contact parameter of water ice on dust ( m=cos(theta) )
170	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
171	! mteta = 0.95
172	write(,) 'water_param contact parameter:', mteta
173
174	c Volume of a water molecule (m3)
175	vo1 = mh2o / dble(rho_ice)
176	c Variance of the ice and CCN distributions
177	sigma_ice = sqrt(log(1.+nuice_sed))
178
179	write(,) 'Variance of ice & CCN distribs :', sigma_ice
180	write(,) 'Volume of a water molecule:', vo1
181
182	firstcall=.false.
183	END IF
184
185	c-----------------------------------------------------------------------
186	c 1. Initialization
187	c-----------------------------------------------------------------------
188	c Update the needed variables
189
190	c write(,) "tauscaling", tauscaling
191
192	!write(,) "pq ccn_mass", pq(ig,:,igcm_ccn_mass)
193	!write(,) "pdq ccn_mass", pdq(ig,:,igcm_ccn_mass)
194	!write(,) "pq ccn_number", pq(ig,:,igcm_ccn_number)
195	!write(,) "pdq ccn_number", pdq(ig,:,igcm_ccn_number)
196
197	! print*, "improvedcloud debut pdq",
198	! & pdq(:,:,igcm_ccn_number)*ptimestep
199	! print*, "improvedcloud debut pq", pq(:,:,igcm_ccn_number)
200
201
202	do l=1,nlay
203	do ig=1,ngrid
204	c Temperature
205	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
206	c Dust mass mixing ratio
207	c (converted to the true value using tauscaling)
208	zq(ig,l,igcm_dust_mass) =
209	& pq(ig,l,igcm_dust_mass) +
210	& pdq(ig,l,igcm_dust_mass) * ptimestep
211	zq(ig,l,igcm_dust_mass) =
212	& zq(ig,l,igcm_dust_mass) * tauscaling(ig)
213	zq(ig,l,igcm_dust_mass)=max(zq(ig,l,igcm_dust_mass),1.E-30)
214	zq0(ig,l,igcm_dust_mass)=zq(ig,l,igcm_dust_mass)
215	c Dust particle number
216	c (converted to the true value using rdust and tauscaling)
217	! zq(ig,l,igcm_dust_number) =
218	! & pq(ig,l,igcm_dust_number) +
219	! & pdq(ig,l,igcm_dust_number) * ptimestep
220	zq(ig,l,igcm_dust_number) =
221	& (1.e0/rdust(ig,l))*3. r3n_q * zq(ig,l,igcm_dust_mass)
222	zq(ig,l,igcm_dust_number)=
223	& max(zq(ig,l,igcm_dust_number),1.E-30)
224	zq0(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number)
225	c CCN mass mixing ratio
226	zq(ig,l,igcm_ccn_mass)=
227	& pq(ig,l,igcm_ccn_mass)+pdq(ig,l,igcm_ccn_mass)*ptimestep
228	zq(ig,l,igcm_ccn_mass)=max(zq(ig,l,igcm_ccn_mass),1.E-30)
229	zq0(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass)
230	c write(,) "pq,zq ccn_mass", pq(ig,l,igcm_ccn_mass),
231	c & zq(ig,l,igcm_ccn_mass)
232	c CCN particle number
233	zq(ig,l,igcm_ccn_number)=
234	& pq(ig,l,igcm_ccn_number)+pdq(ig,l,igcm_ccn_number)*ptimestep
235	zq(ig,l,igcm_ccn_number)=max(zq(ig,l,igcm_ccn_number),1.E-30)
236	zq0(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)
237	c write(,) "pq,zq ccn_number", pq(ig,l,igcm_ccn_number),
238	c & zq(ig,l,igcm_ccn_number)
239	c Water vapor
240	zq(ig,l,igcm_h2o_vap)=
241	& pq(ig,l,igcm_h2o_vap)+pdq(ig,l,igcm_h2o_vap)*ptimestep
242	zq(ig,l,igcm_h2o_vap)=max(zq(ig,l,igcm_h2o_vap),1.E-30) ! FF 12/2004
243	zq0(ig,l,igcm_h2o_vap)=zq(ig,l,igcm_h2o_vap)
244	c Water ice
245	zq(ig,l,igcm_h2o_ice)=
246	& pq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep
247	zq(ig,l,igcm_h2o_ice)=max(zq(ig,l,igcm_h2o_ice),0.) ! FF 12/2004
248	zq0(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
249	enddo
250	enddo
251
252	!print*, "improvedcloud debut pq", pq(1,:,igcm_dust_number)
253
254	c------------------------------------------------------------------
255	c Cloud microphysical scheme
256	c------------------------------------------------------------------
257
258	Cste = ptimestep * 4. * pi * rho_ice
259
260	call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
261
262	c write(,) "ccn_number avant loop phy",
263	c & zq(ig,:,igcm_ccn_number)
264	c write(,) "ccn_mass avant loop phy",
265	c & zq(ig,:,igcm_ccn_mass)
266
267	c Main loop over the GCM's grid
268	DO l=1,nlay
269	!ig = 1
270	DO ig=1,ngrid
271
272	c Get the partial pressure of water vapor and its saturation ratio
273	ph2o = zq(ig,l,igcm_h2o_vap) * (44./18.) * pplay(ig,l)
274	satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
275	satu_out(ig,l) = satu
276	!write(,) "l \| h2o_vap \| zqsat \| satu \| ph2o"
277	!write(,) l,zq(ig,l,igcm_h2o_vap), zqsat(ig,l), satu, ph2o
278
279	c Expand the dust moments into a binned distribution
280	Mo = zq(ig,l,igcm_dust_mass)
281	No = zq(ig,l,igcm_dust_number)+ 1.e-30
282	Rn = rdust(ig,l)
283	Rm = Rn * exp( 3. * sigma_ice**2. )
284	Rn = 1. / Rn
285	Rm = 1. / Rm
286	dev2 = 1. / ( sqrt(2.) * sigma_ice )
287	do i = 1, nbin_cld
288	n_aer(i) = 0.5 * No * ( derf( dlog(rb_cld(i+1)Rn) dev2 )
289	& -derf( dlog(rb_cld(i) * Rn) * dev2 ) )
290	m_aer(i) = 0.5 * Mo * ( derf( dlog(rb_cld(i+1)Rm) dev2 )
291	& -derf( dlog(rb_cld(i) * Rm) * dev2 ) )
292	enddo
293
294	! sumcheck = 0
295	! do i = 1, nbin_cld
296	! sumcheck = sumcheck + n_aer(i)
297	! enddo
298	! sumcheck = abs(sumcheck/No - 1)
299	! if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
300	! print*, "WARNING, No sumcheck PROBLEM"
301	! print*, "sumcheck, No",sumcheck, No
302	! print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
303	! print*, "Dust binned distribution", n_aer
304	! endif
305	!
306	! sumcheck = 0
307	! do i = 1, nbin_cld
308	! sumcheck = sumcheck + M_aer(i)
309	! enddo
310	! sumcheck = abs(sumcheck/Mo - 1)
311	! if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) then
312	! print*, "WARNING, Mo sumcheck PROBLEM"
313	! print*, "sumcheck, No",sumcheck, No
314	! print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
315	! print*, "Dust binned distribution", m_aer
316	! endif
317
318	c Get the rates of nucleation
319	call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
320
321	dN = 0.
322	dM = 0.
323	do i = 1, nbin_cld
324	n_aer(i) = n_aer(i) / ( 1. + rate(i)*ptimestep )
325	m_aer(i) = m_aer(i) / ( 1. + rate(i)*ptimestep )
326	dN = dN + n_aer(i) * rate(i) * ptimestep
327	dM = dM + m_aer(i) * rate(i) * ptimestep
328	enddo
329
330	! dN = min( max(dN,-zq(ig,l,igcm_ccn_number) ),
331	! & zq(ig,l,igcm_dust_number) )
332	!
333	! dM = min( max(dM,-zq(ig,l,igcm_ccn_mass) ),
334	! & zq(ig,l,igcm_dust_mass) )
335
336
337	c IF (zq(ig,l,igcm_ccn_number).ge.1.e-20) THEN
338
339	Mo = zq0(ig,l,igcm_h2o_ice) +
340	& zq0(ig,l,igcm_ccn_mass) + 1.e-30
341	No = zq0(ig,l,igcm_ccn_number)
342	!write(,) "l,cloud particles,cloud mass",l, No, Mo
343	rhocloud(ig,l) = zq0(ig,l,igcm_h2o_ice) / Mo * rho_ice
344	& +zq0(ig,l,igcm_ccn_mass) / Mo * rho_dust
345	rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
346	rice(ig,l) =
347	& ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.)
348	nuice(ig,l)=nuice_ref ! used for rad. transfer calculations
349	if (Mo.lt.1.e-20) rice(ig,l) = 1.e-8
350	seq = exp( 2.sig(zt(ig,l))mh2o /
351	& (rho_icergpzt(ig,l)*rice(ig,l)) )
352
353	call growthrate(ptimestep,zt(ig,l),pplay(ig,l),
354	& ph2o,ph2o/satu,seq,rice(ig,l),gr)
355
356	up = Cste * gr * rice(ig,l) * No * seq +
357	& zq(ig,l,igcm_h2o_vap)
358	dwn = Cste * gr * rice(ig,l) * No / zqsat(ig,l)+ 1.
359
360	Ctot = zq0(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_h2o_vap)
361
362	newvap = min(up/dwn,Ctot)
363
364	gr = gr * ( newvap/zqsat(ig,l) - seq )
365
366
367	dMice = min( max(Cste * No * rice(ig,l) * gr,
368	& -zq(ig,l,igcm_h2o_ice) ),
369	& zq(ig,l,igcm_h2o_vap) )
370
371	c----------- TESTS 1D output ---------
372	if (ngrid.eq.1) then
373	Mcon_out(ig,l) = dMice
374	newvap_out(ig,l) = newvap
375	gr_out(ig,l) = gr
376	dN_out(ig,l) = dN
377	dM_out(ig,l) = dM
378	endif
379	c-------------------------------------
380
381	c Water ice
382	zq(ig,l,igcm_h2o_ice) = zq0(ig,l,igcm_h2o_ice) +
383	& dMice
384
385	c Water vapor
386	zq(ig,l,igcm_h2o_vap) = zq0(ig,l,igcm_h2o_vap) -
387	& dMice
388
389
390	c If all the ice particles sublimate, all the condensation
391	c nuclei are released:
392	if (zq(ig,l,igcm_h2o_ice).le.1e-30) then
393	c Water ice particles
394	zq(ig,l,igcm_h2o_ice) = 0.
395	c Dust particles
396	zq(ig,l,igcm_dust_mass ) = zq0(ig,l,igcm_dust_mass) +
397	& zq0(ig,l,igcm_ccn_mass)
398	zq(ig,l,igcm_dust_number ) = zq0(ig,l,igcm_dust_number) +
399	& zq0(ig,l,igcm_ccn_number)
400	c CCNs
401	zq(ig,l,igcm_ccn_mass) = 0.
402	zq(ig,l,igcm_ccn_number) = 0.
403	c for coherence
404	dM = 0
405	dN = 0
406	endif
407	c ELSE
408	cc Initialize rhocloud and rice to avoid divisions by 0
409	c rhocloud(ig,l) = 1.e-10
410	c rice(ig,l) = 1.e-8
411	c dM = 0
412	c dN = 0
413	c ENDIF ! of if (zq(ig,l,igcm_ccn_number).ge.1.e-20)
414
415	c Dust particles
416	zq(ig,l,igcm_dust_mass ) = zq(ig,l,igcm_dust_mass ) - dM
417	zq(ig,l,igcm_dust_number ) = zq(ig,l,igcm_dust_number ) - dN
418	c CCNs
419	zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + dM
420	zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + dN
421
422	ENDDO
423	ENDDO
424
425
426	! print*, "improvedclouds zq0 abs.", zq0(:,:,igcm_ccn_number)
427	! print*, "improvedclouds zq abs.", zq(:,:,igcm_ccn_number)
428	c------------------------------------------------------------------
429	c Convert the initial values back into relative values
430	c (has to be done before updating rdust!)
431	c------------------------------------------------------------------
432
433	do l=1, nlay
434	do ig=1,ngrid
435	zq0(ig,l,igcm_dust_mass) =
436	& zq0(ig,l,igcm_dust_mass) / tauscaling(ig)
437	zq0(ig,l,igcm_dust_number) =
438	& (1.e0/rdust(ig,l))*3. r3n_q *
439	& zq0(ig,l,igcm_dust_mass)
440	enddo ! of do ig=1,ngrid
441	enddo ! of do l=1,nlay
442
443	c------------------------------------------------------------------
444	c Update the dust radius
445	c------------------------------------------------------------------
446
447	DO l=1,nlay
448	DO ig=1,ngrid
449	rdust(ig,l)=
450	& CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
451	& max(zq(ig,l,igcm_dust_number),0.01))
452	rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
453	ENDDO
454	ENDDO
455
456	c------------------------------------------------------------------
457	c Convert zq back into relative values (only applies to dust)
458	c------------------------------------------------------------------
459
460	do l=1, nlay
461	do ig=1,ngrid
462	c Dust mass mixing ratio
463	c (converted back into relative value using tauscaling)
464	zq(ig,l,igcm_dust_mass) =
465	& zq(ig,l,igcm_dust_mass) / tauscaling(ig)
466	zq(ig,l,igcm_dust_mass)=max(zq(ig,l,igcm_dust_mass),1.E-30)
467	c Dust particle number
468	c (converted back into relative value)
469	zq(ig,l,igcm_dust_number) =
470	& (1.e0/rdust(ig,l))*3. r3n_q * zq(ig,l,igcm_dust_mass)
471	zq(ig,l,igcm_dust_number)=
472	& max(zq(ig,l,igcm_dust_number),1.E-30)
473	enddo ! of do ig=1,ngrid
474	enddo ! of do l=1,nlay
475
476	! print*, "improvedclouds zq0 rel.", zq0(1,:,igcm_ccn_number)
477	! print*, "improvedclouds zq rel.", zq(1,:,igcm_ccn_number)
478
479	c------------------------------------------------------------------
480	c Compute the sedimentation radius
481	c------------------------------------------------------------------
482
483	do l=1, nlay
484	do ig=1,ngrid
485	rsedcloud(ig,l)=max( rice(ig,l)(1.+nuice_sed)*3.,
486	& rdust(ig,l) )
487	rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
488	enddo ! of do ig=1,ngrid
489	enddo ! of do l=1,nlay
490
491	c------------------------------------------------------------------
492	c Force positive values
493	c------------------------------------------------------------------
494
495	! do l=1, nlay
496	! do ig=1,ngrid
497	! zq(ig,l,igcm_ccn_mass)=
498	! & max(zq(ig,l,igcm_ccn_mass),1e-30)
499	! zq(ig,l,igcm_ccn_number)=
500	! & max(zq(ig,l,igcm_ccn_number),1e-30)
501	! zq(ig,l,igcm_h2o_vap)=
502	! & max(zq(ig,l,igcm_h2o_vap),1e-30)
503	! zq(ig,l,igcm_h2o_ice)=
504	! & max(zq(ig,l,igcm_h2o_ice),1e-30)
505	! enddo ! of do ig=1,ngrid
506	! enddo ! of do l=1,nlay
507
508
509	c------------------------------------------------------------------
510	c Compute the final tendencies
511	c------------------------------------------------------------------
512
513	c Initialize the tendencies
514	pdqscloud(1:ngrid,1:nq)=0
515	pdqcloud(1:ngrid,1:nlay,1:nq)=0
516	pdtcloud(1:ngrid,1:nlay)=0
517
518	c Update the tendencies
519	do l=1, nlay
520	do ig=1,ngrid
521	pdqcloud(ig,l,igcm_dust_mass)=(zq(ig,l,igcm_dust_mass)
522	& -zq0(ig,l,igcm_dust_mass))/ptimestep
523	pdqcloud(ig,l,igcm_dust_number)=(zq(ig,l,igcm_dust_number)
524	& -pq(ig,l,igcm_dust_number))/ptimestep
525	& - pdq(ig,l,igcm_dust_number) !!! AJOUT TN
526	pdqcloud(ig,l,igcm_ccn_mass)=(zq(ig,l,igcm_ccn_mass)
527	& -zq0(ig,l,igcm_ccn_mass))/ptimestep
528	pdqcloud(ig,l,igcm_ccn_number)=(zq(ig,l,igcm_ccn_number)
529	& -zq0(ig,l,igcm_ccn_number))/ptimestep
530	pdqcloud(ig,l,igcm_h2o_vap)=(zq(ig,l,igcm_h2o_vap)
531	& -zq0(ig,l,igcm_h2o_vap))/ptimestep
532	pdqcloud(ig,l,igcm_h2o_ice)=(zq(ig,l,igcm_h2o_ice)
533	& -zq0(ig,l,igcm_h2o_ice))/ptimestep
534	lw=(2834.3-0.28(zt(ig,l)-To)-0.004(zt(ig,l)-To)*2)1.e+3
535	pdtcloud(ig,l)=-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp
536	end do
537	end do
538	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
539	! & satu_out)
540
541	! print, "improvedclouds pdqdt",
542	! & pdqcloud(:,:,igcm_ccn_number)*ptimestep
543
544
545	c------------------------------------------------------------------
546	c TEST_JBM
547	IF (ngrid.eq.1) THEN
548	c call WRITEDIAGFI(ngrid,"tausca","tauscaling","",0,
549	c & tauscaling)
550	call WRITEDIAGFI(ngrid,"newvap","h2o newvap","kg",1,
551	& newvap_out)
552	call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
553	& gr_out)
554	call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
555	& dM_out)
556	call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
557	& dN_out)
558	call WRITEDIAGFI(ngrid,"mcond","h2o condensed mass","kg",1,
559	& Mcon_out)
560	call WRITEDIAGFI(ngrid,"zqsat","p vap sat","kg/kg",1,
561	& zqsat)
562	call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
563	& satu_out)
564	call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
565	& rice)
566	call WRITEDIAGFI(ngrid,"rdust","rdust","m",1,
567	& rdust)
568	call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
569	& rsedcloud)
570	call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
571	& rhocloud)
572	ENDIF
573	c------------------------------------------------------------------
574	return
575	end

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

Download in other formats:

Original Format