Context Navigation

improvedclouds_mod.F @ 3525

Last change on this file since 3525 was 3127, checked in by csegonne, 2 years ago

MARS PCM:

Bug fixes for improvedcloud_mod.F :
Initialization of dMice when low ccn number
Initialization of zdq at first sub time step (spenttime=0)
Move 'coef' to fix the formula for zimicro

File size: 26.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	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	ELSE
519	! Initialization of dMice when it's not computed
520	dMice=0 ! no condensation/sublimation to account for
521	ENDIF !of if Nccn>1
522
523
524	! No more getting tendency : we increment tracers & temp directly
525
526	! Increment tracers & temp for following microtimestep
527	! Dust :
528	! Special treatment of dust_mass / number for scavenging ?
529	IF (scavenging) THEN
530	zq(ig,l,igcm_dust_mass) =zq(ig,l,igcm_dust_mass)+
531	& pdq(ig,l,igcm_dust_mass)*microtimestep
532	zq(ig,l,igcm_dust_number) =zq(ig,l,igcm_dust_number)+
533	& pdq(ig,l,igcm_dust_number)*microtimestep
534	ELSE
535	zq(ig,l,igcm_dust_mass) =zq0(ig,l,igcm_dust_mass)
536	zq(ig,l,igcm_dust_number) =zq0(ig,l,igcm_dust_number)
537	ENDIF !(scavenging) THEN
538	zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) +
539	& pdq(ig,l,igcm_ccn_mass)*microtimestep
540	zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) +
541	& pdq(ig,l,igcm_ccn_number)*microtimestep
542
543	! Water :
544	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+
545	& pdq(ig,l,igcm_h2o_ice)*microtimestep
546	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)+
547	& pdq(ig,l,igcm_h2o_vap)*microtimestep
548
549	! HDO (if computed) :
550	IF (hdo) THEN
551	zq(ig,l,igcm_hdo_ice) =
552	& zq(ig,l,igcm_hdo_ice)+
553	& pdq(ig,l,igcm_hdo_ice)*microtimestep
554	zq(ig,l,igcm_hdo_vap) =
555	& zq(ig,l,igcm_hdo_vap)+
556	& pdq(ig,l,igcm_hdo_vap)*microtimestep
557	ENDIF ! hdo
558
559	c Could also set subpdtcloud to 0 if not activice to make it simpler
560	c or change name of the flag
561	IF (.not.activice) THEN
562	subpdtcloud(ig,l)=0.
563	ENDIF
564	! Temperature :
565	zt(ig,l) = zt(ig,l)+(pdt(ig,l)+
566	& subpdtcloud(ig,l))*microtimestep
567
568	c Prevent negative tracers ! JN
569	DO i=1,nq
570	IF(zq(ig,l,i).lt.1.e-30) zq(ig,l,i)=1.e-30
571	ENDDO !i=1,nq
572
573	IF (cloud_adapt_ts) then
574	c Estimation of how much is actually condensing/subliming
575	IF (spenttime.ne.0) then
576	zdq=(dMice/spenttime)*(ptimestep-spenttime)
577	ELSE
578	!Initialization for spenttime=0
579	zdq=zpotcond(ig,l)*((ptimestep-spenttime)/
580	& ptimestep)
581	ENDIF
582	c Threshold with powerlaw (sanity check)
583	zdq=min(abs(zdq),
584	& abs(zpotcond(ig,l)*((ptimestep-spenttime)/ptimestep)))
585	call adapt_imicro(ptimestep,zdq,
586	& zimicro(ig,l))
587	ENDIF! (cloud_adapt_ts) then
588	c Increment time spent in here
589	spenttime=spenttime+microtimestep
590	count_micro(ig,l)=count_micro(ig,l)+1
591	ENDDO ! while (.not. ending_ts)
592	ENDDO ! of ig loop
593	ENDDO ! of nlayer loop
594
595
596	c------ Useful outputs to check how it went
597	call write_output("zpotcond_inst","zpotcond_inst "//
598	& "microphysics","(kg/kg)",zpotcond_inst(:,:))
599	call write_output("zpotcond_full","zpotcond_full "//
600	& "microphysics","(kg/kg)",zpotcond_full(:,:))
601	call write_output("zpotcond","zpotcond "//
602	& "microphysics","(kg/kg)",zpotcond(:,:))
603	call write_output("count_micro","count_micro "//
604	& "after microphysics","integer",count_micro(:,:))
605
606
607
608	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
609	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
610	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
611	IF (test_flag) then
612
613	! error2d(:) = 0.
614	DO l=1,nlay
615	DO ig=1,ngrid
616	! error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
617	satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l)
618	satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
619	ENDDO
620	ENDDO
621
622	print*, 'count is ',countcells, ' i.e. ',
623	& countcells100/(nlayngrid), '% for microphys computation'
624
625	ENDIF ! endif test_flag
626	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
627	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
628	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
629
630	return
631
632
633
634
635	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
636	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
637	c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0
638	c It is an analytical solution to the ice radius growth equation,
639	c with the approximation of a constant 'reduced' cunningham correction factor
640	c (lambda in growthrate.F) taken at radius req instead of rice
641	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
642	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
643
644	c subroutine phi(rice,req,coeff1,coeff2,time)
645	c
646	c implicit none
647	c
648	c ! inputs
649	c real rice ! ice radius
650	c real req ! ice radius at equilibirum
651	c real coeff1 ! coeff for the log
652	c real coeff2 ! coeff for the arctan
653	c
654	c ! output
655	c real time
656	c
657	c !local
658	c real var
659	c
660	c ! 1.73205 is sqrt(3)
661	c
662	c var = max(
663	c & abs(rice-req) / sqrt(ricerice + ricereq + req*req),1e-30)
664	c
665	c time =
666	c & coeff1 *
667	c & log( var )
668	c & + coeff2 * 1.73205 *
669	c & atan( (2rice+req) / (1.73205req) )
670	c
671	c return
672	c end
673
674
675
676	END SUBROUTINE improvedclouds
677
678	SUBROUTINE adapt_imicro(ptimestep,potcond,
679	$ zimicro)
680
681	c Adaptative timestep for water ice clouds.
682	c Works using a powerlaw to compute the minimal duration of subtimestep
683	c (in s) should all the avalaible vapor (resp. ice) condenses (resp.sublimates)
684	c Then, we use the instantaneous vap (ice) gradient extrapolated to the
685	c rest of duration to increase subtimestep duration, for computing
686	c efficiency.
687
688	real,intent(in) :: ptimestep ! total duration of physics (sec)
689	real,intent(in) :: potcond ! condensible vapor / sublimable ice(kg/kg)
690	real :: alpha, beta ! Coefficients for power law
691	real :: defstep,coef ! Default ptimestep of 7.5 mins (iphysiq=5)
692	integer,intent(out) :: zimicro ! number of ptimestep division
693
694	c Default ptimestep : defstep (7.5 mins)
695	defstep=88775.*5./960.
696	coef=ptimestep/defstep
697	c Conservative coefficients :
698	alpha=1.81846504e+11
699	beta=1.54550140e+00
700	zimicro=ceiling(coefmin(max(alphaabs(potcond)**beta,5.),7000.))
701
702	END SUBROUTINE adapt_imicro
703
704	END MODULE improvedclouds_mod

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

Download in other formats:

Original Format