Context Navigation

improvedclouds_mod.F @ 3011

Last change on this file since 3011 was 3008, checked in by emillour, 2 years ago
Mars PCM: Some code cleanup around microphysics. Turn microphys.h into module microphys_h.F90, and while at it also turn nuclea.F, growthrate.F90 and massflowrateco2.F90 into modules. EM
File size: 29.8 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
400	call updaterice_micro(zq(ig,l,igcm_h2o_ice),
401	& zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
402	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
403
404	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
405
406	c saturation at equilibrium
407	c rice should not be too small, otherwise seq value is not valid
408	seq = exp(2.sig(zt(ig,l))mh2o / (rho_icergpzt(ig,l)*
409	& max(rice(ig,l),1.e-7)))
410
411	c get resistance growth
412	call growthrate(zt(ig,l),pplay(ig,l),
413	& real(ph2o/satu),rice(ig,l),res,Dv)
414
415	res_out(ig,l) = res
416
417	ccccccc implicit scheme of mass growth
418	c cste here must be computed at each step
419	cste = 4pirho_ice*microtimestep
420	c print*, 'ig',ig
421	c print*, 'l',l
422	c print*, 'cste',cste
423	c print*, 'seq',seq
424	c print*, 'zqsat',zqsat(ig,l)
425	c print*, 'zq vap',zq(ig,l,igcm_h2o_vap)
426	c print*, 'res',res
427	c print*, 'No',No
428	c print*, 'rice',rice(ig,l)
429
430	dMice =
431	& (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))
432	& /(reszqsat(ig,l)/(csteNo*rice(ig,l)) + 1.)
433
434
435	! With the above scheme, dMice cannot be bigger than vapor,
436	! but can be bigger than all available ice.
437	dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
438	dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless...
439
440	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice
441	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice
442
443
444	countcells = countcells + 1
445
446	! latent heat release
447	lw=(2834.3-0.28*(zt(ig,l)-To)-
448	& 0.004(zt(ig,l)-To)(zt(ig,l)-To))*1.e+3
449	subpdtcloud(ig,l)= dMice*lw/cpp/microtimestep
450
451	c Special case of the isotope of water HDO
452	if (hdo) then
453	!! condensation
454	if (dMice.gt.0.0) then
455	!! do we use fractionation?
456	if (hdofrac) then
457	!! Calculation of the HDO vapor coefficient
458	Dv_hdo = 1./3. * sqrt( 8kbzzt(ig,l)/(pi*mhdo/nav) )
459	& * kbz * zt(ig,l) /
460	& ( pi * pplay(ig,l) * (molco2+molhdo)*(molco2+molhdo)
461	& * sqrt(1.+mhdo/mco2) )
462	!! Calculation of the fractionnation coefficient at equilibrium
463	alpha(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2)
464	c alpha = exp(13525./zt(ig,l)**2.-5.59d-2) !Lamb
465	!! Calculation of the 'real' fractionnation coefficient
466	alpha_c(ig,l) = (alpha(ig,l)*satu)/
467	& ( (alpha(ig,l)(Dv/Dv_hdo)(satu-1.)) + 1.)
468	c alpha_c(ig,l) = alpha(ig,l) ! to test without the effect of cinetics
469	else
470	alpha_c(ig,l) = 1.d0
471	endif
472	if (zq0(ig,l,igcm_h2o_vap).gt.qperemin) then
473	dMice_hdo=
474	& dMicealpha_c(ig,l)
475	& ( zq0(ig,l,igcm_hdo_vap)
476	& /zq0(ig,l,igcm_h2o_vap) )
477	else
478	dMice_hdo=0.
479	endif
480	!! sublimation
481	else
482	if (zq0(ig,l,igcm_h2o_ice).gt.qperemin) then
483	dMice_hdo=
484	& dMice*
485	& ( zq0(ig,l,igcm_hdo_ice)
486	& /zq0(ig,l,igcm_h2o_ice) )
487	else
488	dMice_hdo=0.
489	endif
490	endif !if (dMice.gt.0.0)
491
492	dMice_hdo = max(dMice_hdo,-zq(ig,l,igcm_hdo_ice))
493	dMice_hdo = min(dMice_hdo,zq(ig,l,igcm_hdo_vap))
494
495	zq(ig,l,igcm_hdo_ice) = zq(ig,l,igcm_hdo_ice)+dMice_hdo
496	zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)-dMice_hdo
497
498	endif ! if (hdo)
499
500	!=============================================================
501	! 5. Dust cores released, tendancies, latent heat, etc ...
502	!=============================================================
503
504	c If all the ice particles sublimate, all the condensation
505	c nuclei are released:
506	if (zq(ig,l,igcm_h2o_ice).le.1.e-28) then
507
508	c Water
509	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)
510	& + zq(ig,l,igcm_h2o_ice)
511	zq(ig,l,igcm_h2o_ice) = 0.
512	if (hdo) then
513	zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)
514	& + zq(ig,l,igcm_hdo_ice)
515	zq(ig,l,igcm_hdo_ice) = 0.
516	endif
517	c Dust particles
518	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass)
519	& + zq(ig,l,igcm_ccn_mass)
520	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number)
521	& + zq(ig,l,igcm_ccn_number)
522	c CCNs
523	zq(ig,l,igcm_ccn_mass) = 0.
524	zq(ig,l,igcm_ccn_number) = 0.
525
526	endif
527
528	ENDIF !of if Nccn>1
529
530
531	! No more getting tendency : we increment tracers & temp directly
532
533	! Increment tracers & temp for following microtimestep
534	! Dust :
535	! Special treatment of dust_mass / number for scavenging ?
536	IF (scavenging) THEN
537	zq(ig,l,igcm_dust_mass) =zq(ig,l,igcm_dust_mass)+
538	& pdq(ig,l,igcm_dust_mass)*microtimestep
539	zq(ig,l,igcm_dust_number) =zq(ig,l,igcm_dust_number)+
540	& pdq(ig,l,igcm_dust_number)*microtimestep
541	ELSE
542	zq(ig,l,igcm_dust_mass) =zq0(ig,l,igcm_dust_mass)
543	zq(ig,l,igcm_dust_number) =zq0(ig,l,igcm_dust_number)
544	ENDIF !(scavenging) THEN
545	zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) +
546	& pdq(ig,l,igcm_ccn_mass)*microtimestep
547	zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) +
548	& pdq(ig,l,igcm_ccn_number)*microtimestep
549
550	! Water :
551	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+
552	& pdq(ig,l,igcm_h2o_ice)*microtimestep
553	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)+
554	& pdq(ig,l,igcm_h2o_vap)*microtimestep
555
556	! HDO (if computed) :
557	IF (hdo) THEN
558	zq(ig,l,igcm_hdo_ice) =
559	& zq(ig,l,igcm_hdo_ice)+
560	& pdq(ig,l,igcm_hdo_ice)*microtimestep
561	zq(ig,l,igcm_hdo_vap) =
562	& zq(ig,l,igcm_hdo_vap)+
563	& pdq(ig,l,igcm_hdo_vap)*microtimestep
564	ENDIF ! hdo
565
566	c Could also set subpdtcloud to 0 if not activice to make it simpler
567	c or change name of the flag
568	IF (.not.activice) THEN
569	subpdtcloud(ig,l)=0.
570	ENDIF
571	! Temperature :
572	zt(ig,l) = zt(ig,l)+(pdt(ig,l)+
573	& subpdtcloud(ig,l))*microtimestep
574
575	c Prevent negative tracers ! JN
576	DO i=1,nq
577	IF(zq(ig,l,i).lt.1.e-30) zq(ig,l,i)=1.e-30
578	ENDDO !i=1,nq
579
580	c Increment time spent in here
581	spenttime=spenttime+microtimestep
582	count_micro(ig,l)=count_micro(ig,l)+1
583	IF (cloud_adapt_ts) then
584	c Estimation of how much is actually condensing/subliming
585	c print*, 'zq',zq(ig,l,igcm_h2o_vap)
586	c print*, 'zq0',zq0(ig,l,igcm_h2o_vap)
587	c print*, 'ptimestep',ptimestep
588	zdq=(zq(ig,l,igcm_h2o_vap)-zq0(ig,l,igcm_h2o_vap))
589	& *ptimestep/microtimestep
590	c Estimation of how much is actually condensing/subliming
591	c zdq=min(abs(dMice)*ptimestep,zpotcond(ig,l))
592	call adapt_imicro(ptimestep,zdq,
593	& zimicro(ig,l))
594	ENDIF! (cloud_adapt_ts) then
595	ENDDO ! while (.not. ending_ts)
596	ENDDO ! of ig loop
597	ENDDO ! of nlayer loop
598
599
600	c------ Useful outputs to check how it went
601	call write_output("zpotcond_inst","zpotcond_inst "//
602	& "microphysics","(kg/kg)",zpotcond_inst(:,:))
603	call write_output("zpotcond_full","zpotcond_full "//
604	& "microphysics","(kg/kg)",zpotcond_full(:,:))
605	call write_output("zpotcond","zpotcond "//
606	& "microphysics","(kg/kg)",zpotcond(:,:))
607	call write_output("count_micro","count_micro "//
608	& "after microphysics","integer",count_micro(:,:))
609
610
611
612	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
613	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
614	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
615	IF (test_flag) then
616
617	! error2d(:) = 0.
618	DO l=1,nlay
619	DO ig=1,ngrid
620	! error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
621	satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l)
622	satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
623	ENDDO
624	ENDDO
625
626	print*, 'count is ',countcells, ' i.e. ',
627	& countcells100/(nlayngrid), '% for microphys computation'
628
629	#ifndef MESOSCALE
630	! IF (ngrid.ne.1) THEN ! 3D
631	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
632	! & satu_out)
633	! call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
634	! & dM_out)
635	! call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
636	! & dN_out)
637	! call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2,
638	! & error2d)
639	! call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3,
640	! & zqsat)
641	! ENDIF
642
643	! IF (ngrid.eq.1) THEN ! 1D
644	! call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1,
645	! & error_out)
646	call WRITEdiagfi(ngrid,"resist","resistance","s/m2",1,
647	& res_out)
648	call WRITEdiagfi(ngrid,"satu_bf","satu before","kg/kg",1,
649	& satubf)
650	call WRITEdiagfi(ngrid,"satu_af","satu after","kg/kg",1,
651	& satuaf)
652	call WRITEdiagfi(ngrid,"vapbf","h2ovap before","kg/kg",1,
653	& zq0(1,1,igcm_h2o_vap))
654	call WRITEdiagfi(ngrid,"vapaf","h2ovap after","kg/kg",1,
655	& zq(1,1,igcm_h2o_vap))
656	call WRITEdiagfi(ngrid,"icebf","h2oice before","kg/kg",1,
657	& zq0(1,1,igcm_h2o_ice))
658	call WRITEdiagfi(ngrid,"iceaf","h2oice after","kg/kg",1,
659	& zq(1,1,igcm_h2o_ice))
660	call WRITEdiagfi(ngrid,"ccnbf","ccn before","/kg",1,
661	& zq0(1,1,igcm_ccn_number))
662	call WRITEdiagfi(ngrid,"ccnaf","ccn after","/kg",1,
663	& zq(1,1,igcm_ccn_number))
664	c call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
665	c & gr_out)
666	c call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
667	c & rate_out)
668	c call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
669	c & dM_out)
670	c call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
671	c & dN_out)
672	call WRITEdiagfi(ngrid,"zqsat","p vap sat","kg/kg",1,
673	& zqsat)
674	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
675	! & satu_out)
676	call WRITEdiagfi(ngrid,"rice","ice radius","m",1,
677	& rice)
678	! call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
679	! & rdust)
680	! call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
681	! & rsedcloud)
682	! call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
683	! & rhocloud)
684	! ENDIF
685	#endif
686
687	ENDIF ! endif test_flag
688	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
689	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
690	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
691
692	return
693
694
695
696
697	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
698	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
699	c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0
700	c It is an analytical solution to the ice radius growth equation,
701	c with the approximation of a constant 'reduced' cunningham correction factor
702	c (lambda in growthrate.F) taken at radius req instead of rice
703	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
704	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
705
706	c subroutine phi(rice,req,coeff1,coeff2,time)
707	c
708	c implicit none
709	c
710	c ! inputs
711	c real rice ! ice radius
712	c real req ! ice radius at equilibirum
713	c real coeff1 ! coeff for the log
714	c real coeff2 ! coeff for the arctan
715	c
716	c ! output
717	c real time
718	c
719	c !local
720	c real var
721	c
722	c ! 1.73205 is sqrt(3)
723	c
724	c var = max(
725	c & abs(rice-req) / sqrt(ricerice + ricereq + req*req),1e-30)
726	c
727	c time =
728	c & coeff1 *
729	c & log( var )
730	c & + coeff2 * 1.73205 *
731	c & atan( (2rice+req) / (1.73205req) )
732	c
733	c return
734	c end
735
736
737
738	END SUBROUTINE improvedclouds
739
740	SUBROUTINE adapt_imicro(ptimestep,potcond,
741	$ zimicro)
742
743	c Adaptative timestep for water ice clouds.
744	c Works using a powerlaw to compute the minimal duration of subtimestep
745	c (in s) should all the avalaible vapor (resp. ice) condenses (resp.sublimates)
746	c Then, we use the instantaneous vap (ice) gradient extrapolated to the
747	c rest of duration to increase subtimestep duration, for computing
748	c efficiency.
749
750	real,intent(in) :: ptimestep ! total duration of physics (sec)
751	real,intent(in) :: potcond ! condensible vapor / sublimable ice(kg/kg)
752	real :: alpha, beta ! Coefficients for power law
753	integer,intent(out) :: zimicro ! number of ptimestep division
754
755	c zimicro = 30
756	c Coefficients good enough for present-day Mars :
757	c alpha = 1.87485684e+09
758	c beta = 1.45655856e+00
759	c Coefficients covering high obliquity scenarios :
760	alpha=3.36711332e+15
761	beta=1.98636414e+00
762	c nimicro=min(max(alphaabs(potcond)*beta,5.),7000.)
763	c zimicro=ceiling((ptimestep/timeleft)*nimicro)
764	c zimicro=ceiling((timeleft/ptimestep)*nimicro)
765	zimicro=ceiling(min(max(alphaabs(potcond)*beta,5.),7000.))
766
767	END SUBROUTINE adapt_imicro
768
769	END MODULE improvedclouds_mod

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

Download in other formats:

Original Format