Context Navigation

improvedclouds_mod.F @ 2937

Last change on this file since 2937 was 2932, checked in by romain.vande, 3 years ago

Mars PCM:
Add a new routine to write output called write_output.
It needs to be imported (for example : use write_output_mod, only: write_output)
Then, it requires only 4 arguments : the name of the variable, its title, its units and the variable itself.
It detects the dimension of the variable and decide to output either in diagfi or diagsoil.
It is also compatible with XIOS (xml file needs to be adapted)
Writediagfi and writediagsoil routines are still available in case.
RV

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

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

Download in other formats:

Original Format