Context Navigation

improvedclouds_mod.F @ 3026

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

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

Download in other formats:

Original Format