Context Navigation

improvedclouds_mod.F @ 2417

Last change on this file since 2417 was 2407, checked in by mvals, 5 years ago

Mars GCM:
Implementation of HDO in the microphysics of water ice clouds:

improvedclouds_mod.F: calculation of the HDO flux
growthrate.F, microphys.h: addings to take into account the effect of diffusion kinetics on fractionation
callsedim_mod.F: sedimentation of HDO as an isotope of water in the microphysics case

File size: 24.2 KB

Line
1	MODULE improvedclouds_mod
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	subroutine improvedclouds(ngrid,nlay,microtimestep,
8	& pplay,pteff,sum_subpdt,
9	& pqeff,sum_subpdq,subpdqcloud,subpdtcloud,
10	& nq,tauscaling)
11	USE updaterad, ONLY: updaterice_micro, updaterccn
12	USE watersat_mod, ONLY: watersat
13	use tracer_mod, only: rho_ice, nuice_sed, igcm_h2o_vap,
14	& igcm_h2o_ice, igcm_dust_mass,
15	& igcm_dust_number, igcm_ccn_mass,
16	& igcm_ccn_number,
17	& igcm_hdo_vap,igcm_hdo_ice,
18	& qperemin
19	use conc_mod, only: mmean
20	use comcstfi_h, only: pi, cpp
21	implicit none
22
23
24	c------------------------------------------------------------------
25	c This routine is used to form clouds when a parcel of the GCM is
26	c saturated. It includes the ability to have supersaturation, a
27	c computation of the nucleation rates, growthrates and the
28	c scavenging of dust particles by clouds.
29	c It is worth noting that the amount of dust is computed using the
30	c dust optical depth computed in aeropacity.F. That's why
31	c the variable called "tauscaling" is used to convert
32	c pq(dust_mass) and pq(dust_number), which are relative
33	c quantities, to absolute and realistic quantities stored in zq.
34	c This has to be done to convert the inputs into absolute
35	c values, but also to convert the outputs back into relative
36	c values which are then used by the sedimentation and advection
37	c schemes.
38
39	c Authors: J.-B. Madeleine, based on the work by Franck Montmessin
40	c (October 2011)
41	c T. Navarro, debug,correction, new scheme (October-April 2011)
42	c A. Spiga, optimization (February 2012)
43	c------------------------------------------------------------------
44	#include "callkeys.h"
45	#include "microphys.h"
46	c------------------------------------------------------------------
47	c Inputs/outputs:
48
49	INTEGER, INTENT(IN) :: ngrid,nlay
50	INTEGER, INTENT(IN) :: nq ! nombre de traceurs
51	REAL, INTENT(IN) :: microtimestep ! pas de temps physique (s)
52	REAL, INTENT(IN) :: pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
53	REAL, INTENT(IN) :: pteff(ngrid,nlay) ! temperature at the middle of the
54	! layers (K)
55	REAL, INTENT(IN) :: sum_subpdt(ngrid,nlay)! tendance temperature des autres
56	! param.
57	REAL, INTENT(IN) :: pqeff(ngrid,nlay,nq) ! traceur (kg/kg)
58	REAL, INTENT(IN) :: sum_subpdq(ngrid,nlay,nq) ! tendance avant condensation
59	! (kg/kg.s-1)
60	REAL, INTENT(IN) :: tauscaling(ngrid) ! Convertion factor for qdust and Ndust
61
62	REAL, INTENT(OUT) :: subpdqcloud(ngrid,nlay,nq) ! tendance de la condensation
63	! H2O(kg/kg.s-1)
64	REAL, INTENT(OUT) :: subpdtcloud(ngrid,nlay) ! tendance temperature due
65	! a la chaleur latente
66
67	c------------------------------------------------------------------
68	c Local variables:
69
70	LOGICAL firstcall
71	DATA firstcall/.true./
72	SAVE firstcall
73
74	REAL*8 derf ! Error function
75	!external derf
76
77	INTEGER ig,l,i
78
79	REAL zq(ngrid,nlay,nq) ! local value of tracers
80	REAL zq0(ngrid,nlay,nq) ! local initial value of tracers
81	REAL zt(ngrid,nlay) ! local value of temperature
82	REAL zqsat(ngrid,nlay) ! saturation
83	REAL lw !Latent heat of sublimation (J.kg-1)
84	REAL cste
85	REAL dMice ! mass of condensed ice
86	REAL dMice_hdo ! mass of condensed HDO ice
87	REAL alpha(ngrid,nlay) ! HDO equilibrium fractionation coefficient (Saturation=1)
88	REAL alpha_c(ngrid,nlay) ! HDO real fractionation coefficient
89	! REAL sumcheck
90	REAL*8 ph2o ! Water vapor partial pressure (Pa)
91	REAL*8 satu ! Water vapor saturation ratio over ice
92	REAL*8 Mo,No
93	REAL*8 Rn, Rm, dev2, n_derf, m_derf
94	REAL*8 n_aer(nbin_cld) ! number conc. of particle/each size bin
95	REAL*8 m_aer(nbin_cld) ! mass mixing ratio of particle/each size bin
96
97	REAL*8 sig ! Water-ice/air surface tension (N.m)
98	EXTERNAL sig
99
100	REAL dN,dM
101	REAL rate(nbin_cld) ! nucleation rate
102	REAL seq
103
104	REAL rice(ngrid,nlay) ! Ice mass mean radius (m)
105	! (r_c in montmessin_2004)
106	REAL rhocloud(ngrid,nlay) ! Cloud density (kg.m-3)
107	REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m)
108
109	REAL res ! Resistance growth
110	REAL Dv,Dv_hdo ! Water/HDO vapor diffusion coefficient
111
112
113	c Parameters of the size discretization
114	c used by the microphysical scheme
115	DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m)
116	DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m)
117	DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.0001e-6
118	! Minimum boundary radius (m)
119	DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m)
120	DOUBLE PRECISION vrat_cld ! Volume ratio
121	DOUBLE PRECISION rb_cld(nbin_cld+1)! boundary values of each rad_cld bin (m)
122	SAVE rb_cld
123	DOUBLE PRECISION dr_cld(nbin_cld) ! width of each rad_cld bin (m)
124	DOUBLE PRECISION vol_cld(nbin_cld) ! particle volume for each bin (m3)
125
126
127	REAL sigma_ice ! Variance of the ice and CCN distributions
128	SAVE sigma_ice
129
130
131
132	c----------------------------------
133	c TESTS
134
135	INTEGER countcells
136
137	LOGICAL test_flag ! flag for test/debuging outputs
138	SAVE test_flag
139
140
141	REAL satubf(ngrid,nlay),satuaf(ngrid,nlay)
142	REAL res_out(ngrid,nlay)
143
144
145	c------------------------------------------------------------------
146
147	! AS: firstcall OK absolute
148	IF (firstcall) THEN
149	!=============================================================
150	! 0. Definition of the size grid
151	!=============================================================
152	c rad_cld is the primary radius grid used for microphysics computation.
153	c The grid spacing is computed assuming a constant volume ratio
154	c between two consecutive bins; i.e. vrat_cld.
155	c vrat_cld is determined from the boundary values of the size grid:
156	c rmin_cld and rmax_cld.
157	c The rb_cld array contains the boundary values of each rad_cld bin.
158	c dr_cld is the width of each rad_cld bin.
159
160	c Volume ratio between two adjacent bins
161	! vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
162	! vrat_cld = exp(vrat_cld)
163	vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
164	vrat_cld = exp(vrat_cld)
165	write(,) "vrat_cld", vrat_cld
166
167	rb_cld(1) = rbmin_cld
168	rad_cld(1) = rmin_cld
169	vol_cld(1) = 4./3. * dble(pi) * rmin_cldrmin_cldrmin_cld
170	! vol_cld(1) = 4./3. * pi * rmin_cldrmin_cldrmin_cld
171
172	do i=1,nbin_cld-1
173	rad_cld(i+1) = rad_cld(i) * vrat_cld**(1./3.)
174	vol_cld(i+1) = vol_cld(i) * vrat_cld
175	enddo
176
177	do i=1,nbin_cld
178	rb_cld(i+1)= ( (2.vrat_cld) / (vrat_cld+1.) )(1./3.)
179	& rad_cld(i)
180	dr_cld(i) = rb_cld(i+1) - rb_cld(i)
181	enddo
182	rb_cld(nbin_cld+1) = rbmax_cld
183	dr_cld(nbin_cld) = rb_cld(nbin_cld+1) - rb_cld(nbin_cld)
184
185	print*, ' '
186	print*,'Microphysics: size bin information:'
187	print*,'i,rb_cld(i), rad_cld(i),dr_cld(i)'
188	print*,'-----------------------------------'
189	do i=1,nbin_cld
190	write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i),
191	& dr_cld(i)
192	enddo
193	write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1)
194	print*,'-----------------------------------'
195
196	do i=1,nbin_cld+1
197	! rb_cld(i) = log(rb_cld(i))
198	rb_cld(i) = log(rb_cld(i)) !! we save that so that it is not computed
199	!! at each timestep and gridpoint
200	enddo
201
202	c Contact parameter of water ice on dust ( m=cos(theta) )
203	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
204	! mteta = 0.95
205	write(,) 'water_param contact parameter:', mteta
206
207	c Volume of a water molecule (m3)
208	vo1 = mh2o / dble(rho_ice)
209	c Variance of the ice and CCN distributions
210	sigma_ice = sqrt(log(1.+nuice_sed))
211
212
213	write(,) 'Variance of ice & CCN distribs :', sigma_ice
214	write(,) 'nuice for sedimentation:', nuice_sed
215	write(,) 'Volume of a water molecule:', vo1
216
217
218	test_flag = .false.
219
220	firstcall=.false.
221	END IF
222
223
224	!=============================================================
225	! 1. Initialisation
226	!=============================================================
227	cste = 4pirho_ice*microtimestep
228
229	res_out(:,:) = 0
230	rice(:,:) = 1.e-8
231
232	c Initialize the tendencies
233	subpdqcloud(1:ngrid,1:nlay,1:nq)=0
234	subpdtcloud(1:ngrid,1:nlay)=0
235
236
237	zt(1:ngrid,1:nlay) =
238	& pteff(1:ngrid,1:nlay) +
239	& sum_subpdt(1:ngrid,1:nlay) * microtimestep
240
241	zq(1:ngrid,1:nlay,1:nq) =
242	& pqeff(1:ngrid,1:nlay,1:nq) +
243	& sum_subpdq(1:ngrid,1:nlay,1:nq) * microtimestep
244
245
246	WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
247	& zq(1:ngrid,1:nlay,1:nq) = 1.e-30
248
249	zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
250
251	!=============================================================
252	! 2. Compute saturation
253	!=============================================================
254
255
256	dev2 = 1. / ( sqrt(2.) * sigma_ice )
257
258	call watersat(ngrid*nlay,zt,pplay,zqsat)
259
260	countcells = 0
261
262	c Main loop over the GCM's grid
263	DO l=1,nlay
264	DO ig=1,ngrid
265
266	c Get the partial pressure of water vapor and its saturation ratio
267	ph2o = zq(ig,l,igcm_h2o_vap) * (mmean(ig,l)/18.) * pplay(ig,l)
268	satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
269
270	!=============================================================
271	! 3. Nucleation
272	!=============================================================
273
274	IF ( satu .ge. 1. ) THEN ! if there is condensation
275
276	call updaterccn(zq(ig,l,igcm_dust_mass),
277	& zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
278
279
280	c Expand the dust moments into a binned distribution
281	Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig) + 1.e-30
282	No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30
283	Rn = rdust(ig,l)
284	Rn = -log(Rn)
285	Rm = Rn - 3. * sigma_ice*sigma_ice
286	n_derf = derf( (rb_cld(1)+Rn) *dev2)
287	m_derf = derf( (rb_cld(1)+Rm) *dev2)
288	do i = 1, nbin_cld
289	n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
290	m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
291	n_derf = derf( (rb_cld(i+1)+Rn) *dev2)
292	m_derf = derf( (rb_cld(i+1)+Rm) *dev2)
293	n_aer(i) = n_aer(i) + 0.5 * No * n_derf
294	m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
295	enddo
296
297	! sumcheck = 0
298	! do i = 1, nbin_cld
299	! sumcheck = sumcheck + n_aer(i)
300	! enddo
301	! sumcheck = abs(sumcheck/No - 1)
302	! if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
303	! print*, "WARNING, No sumcheck PROBLEM"
304	! print*, "sumcheck, No",sumcheck, No
305	! print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
306	! print*, "Dust binned distribution", n_aer
307	! endif
308	!
309	! sumcheck = 0
310	! do i = 1, nbin_cld
311	! sumcheck = sumcheck + m_aer(i)
312	! enddo
313	! sumcheck = abs(sumcheck/Mo - 1)
314	! if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) then
315	! print*, "WARNING, Mo sumcheck PROBLEM"
316	! print*, "sumcheck, Mo",sumcheck, Mo
317	! print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
318	! print*, "Dust binned distribution", m_aer
319	! endif
320
321
322	c Get the rates of nucleation
323	call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
324
325	dN = 0.
326	dM = 0.
327	do i = 1, nbin_cld
328	n_aer(i) = n_aer(i)/( 1. + rate(i)*microtimestep)
329	m_aer(i) = m_aer(i)/( 1. + rate(i)*microtimestep)
330	dN = dN + n_aer(i) * rate(i) * microtimestep
331	dM = dM + m_aer(i) * rate(i) * microtimestep
332	enddo
333
334
335	c Update Dust particles
336	zq(ig,l,igcm_dust_mass) =
337	& zq(ig,l,igcm_dust_mass) - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
338	zq(ig,l,igcm_dust_number) =
339	& zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
340	c Update CCNs
341	zq(ig,l,igcm_ccn_mass) =
342	& zq(ig,l,igcm_ccn_mass) + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
343	zq(ig,l,igcm_ccn_number) =
344	& zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
345
346	ENDIF ! of is satu >1
347
348	!=============================================================
349	! 4. Ice growth: scheme for radius evolution
350	!=============================================================
351
352	c We trigger crystal growth if and only if there is at least one nuclei (N>1).
353	c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
354	c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
355
356	IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.) THEN ! we trigger crystal growth
357
358
359	call updaterice_micro(zq(ig,l,igcm_h2o_ice),
360	& zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
361	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
362
363	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
364
365	c saturation at equilibrium
366	c rice should not be too small, otherwise seq value is not valid
367	seq = exp(2.sig(zt(ig,l))mh2o / (rho_icergpzt(ig,l)*
368	& max(rice(ig,l),1.e-7)))
369
370	c get resistance growth
371	call growthrate(zt(ig,l),pplay(ig,l),
372	& real(ph2o/satu),rice(ig,l),res,Dv)
373
374	res_out(ig,l) = res
375
376	ccccccc implicit scheme of mass growth
377
378	dMice =
379	& (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))
380	& /(reszqsat(ig,l)/(csteNo*rice(ig,l)) + 1.)
381
382
383	! With the above scheme, dMice cannot be bigger than vapor,
384	! but can be bigger than all available ice.
385	dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
386	dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless...
387
388	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice
389	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice
390
391
392	countcells = countcells + 1
393
394	! latent heat release
395	lw=(2834.3-0.28*(zt(ig,l)-To)-
396	& 0.004(zt(ig,l)-To)(zt(ig,l)-To))*1.e+3
397	subpdtcloud(ig,l)= dMice*lw/cpp/microtimestep
398
399	c Special case of the isotope of water HDO
400	if (hdo) then
401	!! condensation
402	if (dMice.gt.0.0) then
403	!! do we use fractionation?
404	if (hdofrac) then
405	!! Calculation of the HDO vapor coefficient
406	Dv_hdo = 1./3. * sqrt( 8kbzzt(ig,l)/(pi*mhdo/nav) )
407	& * kbz * zt(ig,l) /
408	& ( pi * pplay(ig,l) * (molco2+molhdo)*(molco2+molhdo)
409	& * sqrt(1.+mhdo/mco2) )
410	!! Calculation of the fractionnation coefficient at equilibrium
411	alpha(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2)
412	c alpha = exp(13525./zt(ig,l)**2.-5.59d-2) !Lamb
413	!! Calculation of the 'real' fractionnation coefficient
414	alpha_c(ig,l) = (alpha(ig,l)*satu)/
415	& ( (alpha(ig,l)(Dv/Dv_hdo)(satu-1.)) + 1.)
416	c alpha_c(ig,l) = alpha(ig,l) ! to test without the effect of cinetics
417	else
418	alpha_c(ig,l) = 1.d0
419	endif
420	if (zq0(ig,l,igcm_h2o_vap).gt.qperemin) then
421	dMice_hdo=
422	& dMicealpha_c(ig,l)
423	& ( zq0(ig,l,igcm_hdo_vap)
424	& /zq0(ig,l,igcm_h2o_vap) )
425	else
426	dMice_hdo=0.
427	endif
428	!! sublimation
429	else
430	if (zq0(ig,l,igcm_h2o_ice).gt.qperemin) then
431	dMice_hdo=
432	& dMice*
433	& ( zq0(ig,l,igcm_hdo_ice)
434	& /zq0(ig,l,igcm_h2o_ice) )
435	else
436	dMice_hdo=0.
437	endif
438	endif !if (dMice.gt.0.0)
439
440	dMice_hdo = max(dMice_hdo,-zq(ig,l,igcm_hdo_ice))
441	dMice_hdo = min(dMice_hdo,zq(ig,l,igcm_hdo_vap))
442
443	zq(ig,l,igcm_hdo_ice) = zq(ig,l,igcm_hdo_ice)+dMice_hdo
444	zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)-dMice_hdo
445
446	endif ! if (hdo)
447
448	!=============================================================
449	! 5. Dust cores released, tendancies, latent heat, etc ...
450	!=============================================================
451
452	c If all the ice particles sublimate, all the condensation
453	c nuclei are released:
454	if (zq(ig,l,igcm_h2o_ice).le.1.e-28) then
455
456	c Water
457	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)
458	& + zq(ig,l,igcm_h2o_ice)
459	zq(ig,l,igcm_h2o_ice) = 0.
460	if (hdo) then
461	zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)
462	& + zq(ig,l,igcm_hdo_ice)
463	zq(ig,l,igcm_hdo_ice) = 0.
464	endif
465	c Dust particles
466	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass)
467	& + zq(ig,l,igcm_ccn_mass)
468	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number)
469	& + zq(ig,l,igcm_ccn_number)
470	c CCNs
471	zq(ig,l,igcm_ccn_mass) = 0.
472	zq(ig,l,igcm_ccn_number) = 0.
473
474	endif
475
476	ENDIF !of if Nccn>1
477
478	ENDDO ! of ig loop
479	ENDDO ! of nlayer loop
480
481
482	! Get cloud tendencies
483	subpdqcloud(1:ngrid,1:nlay,igcm_h2o_vap) =
484	& (zq(1:ngrid,1:nlay,igcm_h2o_vap) -
485	& zq0(1:ngrid,1:nlay,igcm_h2o_vap))/microtimestep
486	subpdqcloud(1:ngrid,1:nlay,igcm_h2o_ice) =
487	& (zq(1:ngrid,1:nlay,igcm_h2o_ice) -
488	& zq0(1:ngrid,1:nlay,igcm_h2o_ice))/microtimestep
489	if (hdo) then
490	subpdqcloud(1:ngrid,1:nlay,igcm_hdo_vap) =
491	& (zq(1:ngrid,1:nlay,igcm_hdo_vap) -
492	& zq0(1:ngrid,1:nlay,igcm_hdo_vap))/microtimestep
493	subpdqcloud(1:ngrid,1:nlay,igcm_hdo_ice) =
494	& (zq(1:ngrid,1:nlay,igcm_hdo_ice) -
495	& zq0(1:ngrid,1:nlay,igcm_hdo_ice))/microtimestep
496	endif
497	subpdqcloud(1:ngrid,1:nlay,igcm_ccn_mass) =
498	& (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
499	& zq0(1:ngrid,1:nlay,igcm_ccn_mass))/microtimestep
500	subpdqcloud(1:ngrid,1:nlay,igcm_ccn_number) =
501	& (zq(1:ngrid,1:nlay,igcm_ccn_number) -
502	& zq0(1:ngrid,1:nlay,igcm_ccn_number))/microtimestep
503
504	if (scavenging) then
505
506	subpdqcloud(1:ngrid,1:nlay,igcm_dust_mass) =
507	& (zq(1:ngrid,1:nlay,igcm_dust_mass) -
508	& zq0(1:ngrid,1:nlay,igcm_dust_mass))/microtimestep
509	subpdqcloud(1:ngrid,1:nlay,igcm_dust_number) =
510	& (zq(1:ngrid,1:nlay,igcm_dust_number) -
511	& zq0(1:ngrid,1:nlay,igcm_dust_number))/microtimestep
512
513	endif
514
515
516
517	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
518	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
519	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
520	IF (test_flag) then
521
522	! error2d(:) = 0.
523	DO l=1,nlay
524	DO ig=1,ngrid
525	! error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
526	satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l)
527	satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
528	ENDDO
529	ENDDO
530
531	print*, 'count is ',countcells, ' i.e. ',
532	& countcells100/(nlayngrid), '% for microphys computation'
533
534	#ifndef MESOSCALE
535	! IF (ngrid.ne.1) THEN ! 3D
536	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
537	! & satu_out)
538	! call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
539	! & dM_out)
540	! call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
541	! & dN_out)
542	! call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2,
543	! & error2d)
544	! call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3,
545	! & zqsat)
546	! ENDIF
547
548	! IF (ngrid.eq.1) THEN ! 1D
549	! call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1,
550	! & error_out)
551	call WRITEdiagfi(ngrid,"resist","resistance","s/m2",1,
552	& res_out)
553	call WRITEdiagfi(ngrid,"satu_bf","satu before","kg/kg",1,
554	& satubf)
555	call WRITEdiagfi(ngrid,"satu_af","satu after","kg/kg",1,
556	& satuaf)
557	call WRITEdiagfi(ngrid,"vapbf","h2ovap before","kg/kg",1,
558	& zq0(1,1,igcm_h2o_vap))
559	call WRITEdiagfi(ngrid,"vapaf","h2ovap after","kg/kg",1,
560	& zq(1,1,igcm_h2o_vap))
561	call WRITEdiagfi(ngrid,"icebf","h2oice before","kg/kg",1,
562	& zq0(1,1,igcm_h2o_ice))
563	call WRITEdiagfi(ngrid,"iceaf","h2oice after","kg/kg",1,
564	& zq(1,1,igcm_h2o_ice))
565	call WRITEdiagfi(ngrid,"ccnbf","ccn before","/kg",1,
566	& zq0(1,1,igcm_ccn_number))
567	call WRITEdiagfi(ngrid,"ccnaf","ccn after","/kg",1,
568	& zq(1,1,igcm_ccn_number))
569	c call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
570	c & gr_out)
571	c call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
572	c & rate_out)
573	c call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
574	c & dM_out)
575	c call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
576	c & dN_out)
577	call WRITEdiagfi(ngrid,"zqsat","p vap sat","kg/kg",1,
578	& zqsat)
579	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
580	! & satu_out)
581	call WRITEdiagfi(ngrid,"rice","ice radius","m",1,
582	& rice)
583	! call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
584	! & rdust)
585	! call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
586	! & rsedcloud)
587	! call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
588	! & rhocloud)
589	! ENDIF
590	#endif
591
592	ENDIF ! endif test_flag
593	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
594	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
595	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
596
597	return
598
599
600
601
602	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
603	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
604	c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0
605	c It is an analytical solution to the ice radius growth equation,
606	c with the approximation of a constant 'reduced' cunningham correction factor
607	c (lambda in growthrate.F) taken at radius req instead of rice
608	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
609	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
610
611	c subroutine phi(rice,req,coeff1,coeff2,time)
612	c
613	c implicit none
614	c
615	c ! inputs
616	c real rice ! ice radius
617	c real req ! ice radius at equilibirum
618	c real coeff1 ! coeff for the log
619	c real coeff2 ! coeff for the arctan
620	c
621	c ! output
622	c real time
623	c
624	c !local
625	c real var
626	c
627	c ! 1.73205 is sqrt(3)
628	c
629	c var = max(
630	c & abs(rice-req) / sqrt(ricerice + ricereq + req*req),1e-30)
631	c
632	c time =
633	c & coeff1 *
634	c & log( var )
635	c & + coeff2 * 1.73205 *
636	c & atan( (2rice+req) / (1.73205req) )
637	c
638	c return
639	c end
640
641
642
643	END SUBROUTINE improvedclouds
644
645	END MODULE improvedclouds_mod

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

Download in other formats:

Original Format