Context Navigation

improvedclouds_mod.F @ 2599

Last change on this file since 2599 was 2437, checked in by mvals, 5 years ago

Mars GCM:

improvedclouds_mod.F: update of the nucleation equation with its analytical resolution
writediagmicrofi.F: makes it possible to get outputs from the microphysics (call example in watercloud_mod.F)

File size: 24.0 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	dN = dN + n_aer(i)(exp(-rate(i)microtimestep)-1.)
329	dM = dM + m_aer(i)(exp(-rate(i)microtimestep)-1.)
330	enddo
331
332
333	c Update Dust particles
334	zq(ig,l,igcm_dust_mass) =
335	& zq(ig,l,igcm_dust_mass) + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
336	zq(ig,l,igcm_dust_number) =
337	& zq(ig,l,igcm_dust_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
338	c Update CCNs
339	zq(ig,l,igcm_ccn_mass) =
340	& zq(ig,l,igcm_ccn_mass) - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
341	zq(ig,l,igcm_ccn_number) =
342	& zq(ig,l,igcm_ccn_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
343
344	ENDIF ! of is satu >1
345
346	!=============================================================
347	! 4. Ice growth: scheme for radius evolution
348	!=============================================================
349
350	c We trigger crystal growth if and only if there is at least one nuclei (N>1).
351	c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
352	c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
353
354	IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.) THEN ! we trigger crystal growth
355
356
357	call updaterice_micro(zq(ig,l,igcm_h2o_ice),
358	& zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
359	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
360
361	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
362
363	c saturation at equilibrium
364	c rice should not be too small, otherwise seq value is not valid
365	seq = exp(2.sig(zt(ig,l))mh2o / (rho_icergpzt(ig,l)*
366	& max(rice(ig,l),1.e-7)))
367
368	c get resistance growth
369	call growthrate(zt(ig,l),pplay(ig,l),
370	& real(ph2o/satu),rice(ig,l),res,Dv)
371
372	res_out(ig,l) = res
373
374	ccccccc implicit scheme of mass growth
375
376	dMice =
377	& (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))
378	& /(reszqsat(ig,l)/(csteNo*rice(ig,l)) + 1.)
379
380
381	! With the above scheme, dMice cannot be bigger than vapor,
382	! but can be bigger than all available ice.
383	dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
384	dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless...
385
386	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice
387	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice
388
389
390	countcells = countcells + 1
391
392	! latent heat release
393	lw=(2834.3-0.28*(zt(ig,l)-To)-
394	& 0.004(zt(ig,l)-To)(zt(ig,l)-To))*1.e+3
395	subpdtcloud(ig,l)= dMice*lw/cpp/microtimestep
396
397	c Special case of the isotope of water HDO
398	if (hdo) then
399	!! condensation
400	if (dMice.gt.0.0) then
401	!! do we use fractionation?
402	if (hdofrac) then
403	!! Calculation of the HDO vapor coefficient
404	Dv_hdo = 1./3. * sqrt( 8kbzzt(ig,l)/(pi*mhdo/nav) )
405	& * kbz * zt(ig,l) /
406	& ( pi * pplay(ig,l) * (molco2+molhdo)*(molco2+molhdo)
407	& * sqrt(1.+mhdo/mco2) )
408	!! Calculation of the fractionnation coefficient at equilibrium
409	alpha(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2)
410	c alpha = exp(13525./zt(ig,l)**2.-5.59d-2) !Lamb
411	!! Calculation of the 'real' fractionnation coefficient
412	alpha_c(ig,l) = (alpha(ig,l)*satu)/
413	& ( (alpha(ig,l)(Dv/Dv_hdo)(satu-1.)) + 1.)
414	c alpha_c(ig,l) = alpha(ig,l) ! to test without the effect of cinetics
415	else
416	alpha_c(ig,l) = 1.d0
417	endif
418	if (zq0(ig,l,igcm_h2o_vap).gt.qperemin) then
419	dMice_hdo=
420	& dMicealpha_c(ig,l)
421	& ( zq0(ig,l,igcm_hdo_vap)
422	& /zq0(ig,l,igcm_h2o_vap) )
423	else
424	dMice_hdo=0.
425	endif
426	!! sublimation
427	else
428	if (zq0(ig,l,igcm_h2o_ice).gt.qperemin) then
429	dMice_hdo=
430	& dMice*
431	& ( zq0(ig,l,igcm_hdo_ice)
432	& /zq0(ig,l,igcm_h2o_ice) )
433	else
434	dMice_hdo=0.
435	endif
436	endif !if (dMice.gt.0.0)
437
438	dMice_hdo = max(dMice_hdo,-zq(ig,l,igcm_hdo_ice))
439	dMice_hdo = min(dMice_hdo,zq(ig,l,igcm_hdo_vap))
440
441	zq(ig,l,igcm_hdo_ice) = zq(ig,l,igcm_hdo_ice)+dMice_hdo
442	zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)-dMice_hdo
443
444	endif ! if (hdo)
445
446	!=============================================================
447	! 5. Dust cores released, tendancies, latent heat, etc ...
448	!=============================================================
449
450	c If all the ice particles sublimate, all the condensation
451	c nuclei are released:
452	if (zq(ig,l,igcm_h2o_ice).le.1.e-28) then
453
454	c Water
455	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)
456	& + zq(ig,l,igcm_h2o_ice)
457	zq(ig,l,igcm_h2o_ice) = 0.
458	if (hdo) then
459	zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)
460	& + zq(ig,l,igcm_hdo_ice)
461	zq(ig,l,igcm_hdo_ice) = 0.
462	endif
463	c Dust particles
464	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass)
465	& + zq(ig,l,igcm_ccn_mass)
466	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number)
467	& + zq(ig,l,igcm_ccn_number)
468	c CCNs
469	zq(ig,l,igcm_ccn_mass) = 0.
470	zq(ig,l,igcm_ccn_number) = 0.
471
472	endif
473
474	ENDIF !of if Nccn>1
475
476	ENDDO ! of ig loop
477	ENDDO ! of nlayer loop
478
479
480	! Get cloud tendencies
481	subpdqcloud(1:ngrid,1:nlay,igcm_h2o_vap) =
482	& (zq(1:ngrid,1:nlay,igcm_h2o_vap) -
483	& zq0(1:ngrid,1:nlay,igcm_h2o_vap))/microtimestep
484	subpdqcloud(1:ngrid,1:nlay,igcm_h2o_ice) =
485	& (zq(1:ngrid,1:nlay,igcm_h2o_ice) -
486	& zq0(1:ngrid,1:nlay,igcm_h2o_ice))/microtimestep
487	if (hdo) then
488	subpdqcloud(1:ngrid,1:nlay,igcm_hdo_vap) =
489	& (zq(1:ngrid,1:nlay,igcm_hdo_vap) -
490	& zq0(1:ngrid,1:nlay,igcm_hdo_vap))/microtimestep
491	subpdqcloud(1:ngrid,1:nlay,igcm_hdo_ice) =
492	& (zq(1:ngrid,1:nlay,igcm_hdo_ice) -
493	& zq0(1:ngrid,1:nlay,igcm_hdo_ice))/microtimestep
494	endif
495	subpdqcloud(1:ngrid,1:nlay,igcm_ccn_mass) =
496	& (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
497	& zq0(1:ngrid,1:nlay,igcm_ccn_mass))/microtimestep
498	subpdqcloud(1:ngrid,1:nlay,igcm_ccn_number) =
499	& (zq(1:ngrid,1:nlay,igcm_ccn_number) -
500	& zq0(1:ngrid,1:nlay,igcm_ccn_number))/microtimestep
501
502	if (scavenging) then
503
504	subpdqcloud(1:ngrid,1:nlay,igcm_dust_mass) =
505	& (zq(1:ngrid,1:nlay,igcm_dust_mass) -
506	& zq0(1:ngrid,1:nlay,igcm_dust_mass))/microtimestep
507	subpdqcloud(1:ngrid,1:nlay,igcm_dust_number) =
508	& (zq(1:ngrid,1:nlay,igcm_dust_number) -
509	& zq0(1:ngrid,1:nlay,igcm_dust_number))/microtimestep
510
511	endif
512
513
514
515	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
516	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
517	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
518	IF (test_flag) then
519
520	! error2d(:) = 0.
521	DO l=1,nlay
522	DO ig=1,ngrid
523	! error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
524	satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l)
525	satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
526	ENDDO
527	ENDDO
528
529	print*, 'count is ',countcells, ' i.e. ',
530	& countcells100/(nlayngrid), '% for microphys computation'
531
532	#ifndef MESOSCALE
533	! IF (ngrid.ne.1) THEN ! 3D
534	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
535	! & satu_out)
536	! call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
537	! & dM_out)
538	! call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
539	! & dN_out)
540	! call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2,
541	! & error2d)
542	! call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3,
543	! & zqsat)
544	! ENDIF
545
546	! IF (ngrid.eq.1) THEN ! 1D
547	! call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1,
548	! & error_out)
549	call WRITEdiagfi(ngrid,"resist","resistance","s/m2",1,
550	& res_out)
551	call WRITEdiagfi(ngrid,"satu_bf","satu before","kg/kg",1,
552	& satubf)
553	call WRITEdiagfi(ngrid,"satu_af","satu after","kg/kg",1,
554	& satuaf)
555	call WRITEdiagfi(ngrid,"vapbf","h2ovap before","kg/kg",1,
556	& zq0(1,1,igcm_h2o_vap))
557	call WRITEdiagfi(ngrid,"vapaf","h2ovap after","kg/kg",1,
558	& zq(1,1,igcm_h2o_vap))
559	call WRITEdiagfi(ngrid,"icebf","h2oice before","kg/kg",1,
560	& zq0(1,1,igcm_h2o_ice))
561	call WRITEdiagfi(ngrid,"iceaf","h2oice after","kg/kg",1,
562	& zq(1,1,igcm_h2o_ice))
563	call WRITEdiagfi(ngrid,"ccnbf","ccn before","/kg",1,
564	& zq0(1,1,igcm_ccn_number))
565	call WRITEdiagfi(ngrid,"ccnaf","ccn after","/kg",1,
566	& zq(1,1,igcm_ccn_number))
567	c call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
568	c & gr_out)
569	c call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
570	c & rate_out)
571	c call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
572	c & dM_out)
573	c call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
574	c & dN_out)
575	call WRITEdiagfi(ngrid,"zqsat","p vap sat","kg/kg",1,
576	& zqsat)
577	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
578	! & satu_out)
579	call WRITEdiagfi(ngrid,"rice","ice radius","m",1,
580	& rice)
581	! call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
582	! & rdust)
583	! call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
584	! & rsedcloud)
585	! call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
586	! & rhocloud)
587	! ENDIF
588	#endif
589
590	ENDIF ! endif test_flag
591	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
592	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
593	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
594
595	return
596
597
598
599
600	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
601	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
602	c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0
603	c It is an analytical solution to the ice radius growth equation,
604	c with the approximation of a constant 'reduced' cunningham correction factor
605	c (lambda in growthrate.F) taken at radius req instead of rice
606	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
607	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
608
609	c subroutine phi(rice,req,coeff1,coeff2,time)
610	c
611	c implicit none
612	c
613	c ! inputs
614	c real rice ! ice radius
615	c real req ! ice radius at equilibirum
616	c real coeff1 ! coeff for the log
617	c real coeff2 ! coeff for the arctan
618	c
619	c ! output
620	c real time
621	c
622	c !local
623	c real var
624	c
625	c ! 1.73205 is sqrt(3)
626	c
627	c var = max(
628	c & abs(rice-req) / sqrt(ricerice + ricereq + req*req),1e-30)
629	c
630	c time =
631	c & coeff1 *
632	c & log( var )
633	c & + coeff2 * 1.73205 *
634	c & atan( (2rice+req) / (1.73205req) )
635	c
636	c return
637	c end
638
639
640
641	END SUBROUTINE improvedclouds
642
643	END MODULE improvedclouds_mod

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

Download in other formats:

Original Format