Context Navigation

improvedclouds_mod.F @ 3804

Last change on this file since 3804 was 3740, checked in by emillour, 10 months ago
Mars PCM: Code tidying; removed unused mufract.F (mucorr.F does the job); turned ambiguously named sig.F to sig_h2o.F90 and make it a module. EM
File size: 26.6 KB

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

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

Download in other formats:

Original Format