Context Navigation

improvedclouds_mod.F @ 2997

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

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

Download in other formats:

Original Format