Context Navigation

improvedclouds_mod.F @ 2616

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

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

Download in other formats:

Original Format