Context Navigation

improvedclouds_mod.F @ 2332

Last change on this file since 2332 was 2304, checked in by emillour, 5 years ago
Mars GCM: Some code tidying: use getin_p() instead of getin() and use "call abort_physic" instead of "stop" or "call abort" EM
File size: 21.5 KB

Line
1	MODULE improvedclouds_mod
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	subroutine improvedclouds(ngrid,nlay,microtimestep,
8	& pplay,pteff,sum_subpdt,
9	& pqeff,sum_subpdq,subpdqcloud,subpdtcloud,
10	& nq,tauscaling)
11	USE updaterad, ONLY: updaterice_micro, updaterccn
12	USE watersat_mod, ONLY: watersat
13	use tracer_mod, only: rho_ice, nuice_sed, igcm_h2o_vap,
14	& igcm_h2o_ice, igcm_dust_mass,
15	& igcm_dust_number, igcm_ccn_mass,
16	& igcm_ccn_number
17	use conc_mod, only: mmean
18	use comcstfi_h, only: pi, cpp
19	implicit none
20
21
22	c------------------------------------------------------------------
23	c This routine is used to form clouds when a parcel of the GCM is
24	c saturated. It includes the ability to have supersaturation, a
25	c computation of the nucleation rates, growthrates and the
26	c scavenging of dust particles by clouds.
27	c It is worth noting that the amount of dust is computed using the
28	c dust optical depth computed in aeropacity.F. That's why
29	c the variable called "tauscaling" is used to convert
30	c pq(dust_mass) and pq(dust_number), which are relative
31	c quantities, to absolute and realistic quantities stored in zq.
32	c This has to be done to convert the inputs into absolute
33	c values, but also to convert the outputs back into relative
34	c values which are then used by the sedimentation and advection
35	c schemes.
36
37	c Authors: J.-B. Madeleine, based on the work by Franck Montmessin
38	c (October 2011)
39	c T. Navarro, debug,correction, new scheme (October-April 2011)
40	c A. Spiga, optimization (February 2012)
41	c------------------------------------------------------------------
42	#include "callkeys.h"
43	#include "microphys.h"
44	c------------------------------------------------------------------
45	c Inputs/outputs:
46
47	INTEGER, INTENT(IN) :: ngrid,nlay
48	INTEGER, INTENT(IN) :: nq ! nombre de traceurs
49	REAL, INTENT(IN) :: microtimestep ! pas de temps physique (s)
50	REAL, INTENT(IN) :: pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
51	REAL, INTENT(IN) :: pteff(ngrid,nlay) ! temperature at the middle of the
52	! layers (K)
53	REAL, INTENT(IN) :: sum_subpdt(ngrid,nlay)! tendance temperature des autres
54	! param.
55	REAL, INTENT(IN) :: pqeff(ngrid,nlay,nq) ! traceur (kg/kg)
56	REAL, INTENT(IN) :: sum_subpdq(ngrid,nlay,nq) ! tendance avant condensation
57	! (kg/kg.s-1)
58	REAL, INTENT(IN) :: tauscaling(ngrid) ! Convertion factor for qdust and Ndust
59
60	REAL, INTENT(OUT) :: subpdqcloud(ngrid,nlay,nq) ! tendance de la condensation
61	! H2O(kg/kg.s-1)
62	REAL, INTENT(OUT) :: subpdtcloud(ngrid,nlay) ! tendance temperature due
63	! a la chaleur latente
64
65	c------------------------------------------------------------------
66	c Local variables:
67
68	LOGICAL firstcall
69	DATA firstcall/.true./
70	SAVE firstcall
71
72	REAL*8 derf ! Error function
73	!external derf
74
75	INTEGER ig,l,i
76
77	REAL zq(ngrid,nlay,nq) ! local value of tracers
78	REAL zq0(ngrid,nlay,nq) ! local initial value of tracers
79	REAL zt(ngrid,nlay) ! local value of temperature
80	REAL zqsat(ngrid,nlay) ! saturation
81	REAL lw !Latent heat of sublimation (J.kg-1)
82	REAL cste
83	REAL dMice ! mass of condensed ice
84	! REAL sumcheck
85	REAL*8 ph2o ! Water vapor partial pressure (Pa)
86	REAL*8 satu ! Water vapor saturation ratio over ice
87	REAL*8 Mo,No
88	REAL*8 Rn, Rm, dev2, n_derf, m_derf
89	REAL*8 n_aer(nbin_cld) ! number conc. of particle/each size bin
90	REAL*8 m_aer(nbin_cld) ! mass mixing ratio of particle/each size bin
91
92	REAL*8 sig ! Water-ice/air surface tension (N.m)
93	EXTERNAL sig
94
95	REAL dN,dM
96	REAL rate(nbin_cld) ! nucleation rate
97	REAL seq
98
99	REAL rice(ngrid,nlay) ! Ice mass mean radius (m)
100	! (r_c in montmessin_2004)
101	REAL rhocloud(ngrid,nlay) ! Cloud density (kg.m-3)
102	REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m)
103
104	REAL res ! Resistance growth
105
106
107	c Parameters of the size discretization
108	c used by the microphysical scheme
109	DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m)
110	DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m)
111	DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.0001e-6
112	! Minimum boundary radius (m)
113	DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m)
114	DOUBLE PRECISION vrat_cld ! Volume ratio
115	DOUBLE PRECISION rb_cld(nbin_cld+1)! boundary values of each rad_cld bin (m)
116	SAVE rb_cld
117	DOUBLE PRECISION dr_cld(nbin_cld) ! width of each rad_cld bin (m)
118	DOUBLE PRECISION vol_cld(nbin_cld) ! particle volume for each bin (m3)
119
120
121	REAL sigma_ice ! Variance of the ice and CCN distributions
122	SAVE sigma_ice
123
124
125
126	c----------------------------------
127	c TESTS
128
129	INTEGER countcells
130
131	LOGICAL test_flag ! flag for test/debuging outputs
132	SAVE test_flag
133
134
135	REAL satubf(ngrid,nlay),satuaf(ngrid,nlay)
136	REAL res_out(ngrid,nlay)
137
138
139	c------------------------------------------------------------------
140
141	! AS: firstcall OK absolute
142	IF (firstcall) THEN
143	!=============================================================
144	! 0. Definition of the size grid
145	!=============================================================
146	c rad_cld is the primary radius grid used for microphysics computation.
147	c The grid spacing is computed assuming a constant volume ratio
148	c between two consecutive bins; i.e. vrat_cld.
149	c vrat_cld is determined from the boundary values of the size grid:
150	c rmin_cld and rmax_cld.
151	c The rb_cld array contains the boundary values of each rad_cld bin.
152	c dr_cld is the width of each rad_cld bin.
153
154	c Volume ratio between two adjacent bins
155	! vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
156	! vrat_cld = exp(vrat_cld)
157	vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
158	vrat_cld = exp(vrat_cld)
159	write(,) "vrat_cld", vrat_cld
160
161	rb_cld(1) = rbmin_cld
162	rad_cld(1) = rmin_cld
163	vol_cld(1) = 4./3. * dble(pi) * rmin_cldrmin_cldrmin_cld
164	! vol_cld(1) = 4./3. * pi * rmin_cldrmin_cldrmin_cld
165
166	do i=1,nbin_cld-1
167	rad_cld(i+1) = rad_cld(i) * vrat_cld**(1./3.)
168	vol_cld(i+1) = vol_cld(i) * vrat_cld
169	enddo
170
171	do i=1,nbin_cld
172	rb_cld(i+1)= ( (2.vrat_cld) / (vrat_cld+1.) )(1./3.)
173	& rad_cld(i)
174	dr_cld(i) = rb_cld(i+1) - rb_cld(i)
175	enddo
176	rb_cld(nbin_cld+1) = rbmax_cld
177	dr_cld(nbin_cld) = rb_cld(nbin_cld+1) - rb_cld(nbin_cld)
178
179	print*, ' '
180	print*,'Microphysics: size bin information:'
181	print*,'i,rb_cld(i), rad_cld(i),dr_cld(i)'
182	print*,'-----------------------------------'
183	do i=1,nbin_cld
184	write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i),
185	& dr_cld(i)
186	enddo
187	write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1)
188	print*,'-----------------------------------'
189
190	do i=1,nbin_cld+1
191	! rb_cld(i) = log(rb_cld(i))
192	rb_cld(i) = log(rb_cld(i)) !! we save that so that it is not computed
193	!! at each timestep and gridpoint
194	enddo
195
196	c Contact parameter of water ice on dust ( m=cos(theta) )
197	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
198	! mteta = 0.95
199	write(,) 'water_param contact parameter:', mteta
200
201	c Volume of a water molecule (m3)
202	vo1 = mh2o / dble(rho_ice)
203	c Variance of the ice and CCN distributions
204	sigma_ice = sqrt(log(1.+nuice_sed))
205
206
207	write(,) 'Variance of ice & CCN distribs :', sigma_ice
208	write(,) 'nuice for sedimentation:', nuice_sed
209	write(,) 'Volume of a water molecule:', vo1
210
211
212	test_flag = .false.
213
214	firstcall=.false.
215	END IF
216
217
218	!=============================================================
219	! 1. Initialisation
220	!=============================================================
221	cste = 4pirho_ice*microtimestep
222
223	res_out(:,:) = 0
224	rice(:,:) = 1.e-8
225
226	c Initialize the tendencies
227	subpdqcloud(1:ngrid,1:nlay,1:nq)=0
228	subpdtcloud(1:ngrid,1:nlay)=0
229
230
231	zt(1:ngrid,1:nlay) =
232	& pteff(1:ngrid,1:nlay) +
233	& sum_subpdt(1:ngrid,1:nlay) * microtimestep
234
235	zq(1:ngrid,1:nlay,1:nq) =
236	& pqeff(1:ngrid,1:nlay,1:nq) +
237	& sum_subpdq(1:ngrid,1:nlay,1:nq) * microtimestep
238
239
240	WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
241	& zq(1:ngrid,1:nlay,1:nq) = 1.e-30
242
243	zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
244
245	!=============================================================
246	! 2. Compute saturation
247	!=============================================================
248
249
250	dev2 = 1. / ( sqrt(2.) * sigma_ice )
251
252	call watersat(ngrid*nlay,zt,pplay,zqsat)
253
254	countcells = 0
255
256	c Main loop over the GCM's grid
257	DO l=1,nlay
258	DO ig=1,ngrid
259
260	c Get the partial pressure of water vapor and its saturation ratio
261	ph2o = zq(ig,l,igcm_h2o_vap) * (mmean(ig,l)/18.) * pplay(ig,l)
262	satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
263
264	!=============================================================
265	! 3. Nucleation
266	!=============================================================
267
268	IF ( satu .ge. 1. ) THEN ! if there is condensation
269
270	call updaterccn(zq(ig,l,igcm_dust_mass),
271	& zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
272
273
274	c Expand the dust moments into a binned distribution
275	Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig) + 1.e-30
276	No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30
277	Rn = rdust(ig,l)
278	Rn = -log(Rn)
279	Rm = Rn - 3. * sigma_ice*sigma_ice
280	n_derf = derf( (rb_cld(1)+Rn) *dev2)
281	m_derf = derf( (rb_cld(1)+Rm) *dev2)
282	do i = 1, nbin_cld
283	n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
284	m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
285	n_derf = derf( (rb_cld(i+1)+Rn) *dev2)
286	m_derf = derf( (rb_cld(i+1)+Rm) *dev2)
287	n_aer(i) = n_aer(i) + 0.5 * No * n_derf
288	m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
289	enddo
290
291	! sumcheck = 0
292	! do i = 1, nbin_cld
293	! sumcheck = sumcheck + n_aer(i)
294	! enddo
295	! sumcheck = abs(sumcheck/No - 1)
296	! if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
297	! print*, "WARNING, No sumcheck PROBLEM"
298	! print*, "sumcheck, No",sumcheck, No
299	! print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
300	! print*, "Dust binned distribution", n_aer
301	! endif
302	!
303	! sumcheck = 0
304	! do i = 1, nbin_cld
305	! sumcheck = sumcheck + m_aer(i)
306	! enddo
307	! sumcheck = abs(sumcheck/Mo - 1)
308	! if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) then
309	! print*, "WARNING, Mo sumcheck PROBLEM"
310	! print*, "sumcheck, Mo",sumcheck, Mo
311	! print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
312	! print*, "Dust binned distribution", m_aer
313	! endif
314
315
316	c Get the rates of nucleation
317	call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
318
319	dN = 0.
320	dM = 0.
321	do i = 1, nbin_cld
322	n_aer(i) = n_aer(i)/( 1. + rate(i)*microtimestep)
323	m_aer(i) = m_aer(i)/( 1. + rate(i)*microtimestep)
324	dN = dN + n_aer(i) * rate(i) * microtimestep
325	dM = dM + m_aer(i) * rate(i) * microtimestep
326	enddo
327
328
329	c Update Dust particles
330	zq(ig,l,igcm_dust_mass) =
331	& zq(ig,l,igcm_dust_mass) - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
332	zq(ig,l,igcm_dust_number) =
333	& zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
334	c Update CCNs
335	zq(ig,l,igcm_ccn_mass) =
336	& zq(ig,l,igcm_ccn_mass) + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
337	zq(ig,l,igcm_ccn_number) =
338	& zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
339
340	ENDIF ! of is satu >1
341
342	!=============================================================
343	! 4. Ice growth: scheme for radius evolution
344	!=============================================================
345
346	c We trigger crystal growth if and only if there is at least one nuclei (N>1).
347	c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
348	c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
349
350	IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.) THEN ! we trigger crystal growth
351
352
353	call updaterice_micro(zq(ig,l,igcm_h2o_ice),
354	& zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
355	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
356
357	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
358
359	c saturation at equilibrium
360	c rice should not be too small, otherwise seq value is not valid
361	seq = exp(2.sig(zt(ig,l))mh2o / (rho_icergpzt(ig,l)*
362	& max(rice(ig,l),1.e-7)))
363
364	c get resistance growth
365	call growthrate(zt(ig,l),pplay(ig,l),
366	& real(ph2o/satu),rice(ig,l),res)
367
368	res_out(ig,l) = res
369
370	ccccccc implicit scheme of mass growth
371
372	dMice =
373	& (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))
374	& /(reszqsat(ig,l)/(csteNo*rice(ig,l)) + 1.)
375
376
377	! With the above scheme, dMice cannot be bigger than vapor,
378	! but can be bigger than all available ice.
379	dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
380	dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless...
381
382	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice
383	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice
384
385
386	countcells = countcells + 1
387
388	! latent heat release
389	lw=(2834.3-0.28*(zt(ig,l)-To)-
390	& 0.004(zt(ig,l)-To)(zt(ig,l)-To))*1.e+3
391	subpdtcloud(ig,l)= dMice*lw/cpp/microtimestep
392
393
394
395	!=============================================================
396	! 5. Dust cores released, tendancies, latent heat, etc ...
397	!=============================================================
398
399	c If all the ice particles sublimate, all the condensation
400	c nuclei are released:
401	if (zq(ig,l,igcm_h2o_ice).le.1.e-28) then
402
403	c Water
404	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)
405	& + zq(ig,l,igcm_h2o_ice)
406	zq(ig,l,igcm_h2o_ice) = 0.
407	c Dust particles
408	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass)
409	& + zq(ig,l,igcm_ccn_mass)
410	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number)
411	& + zq(ig,l,igcm_ccn_number)
412	c CCNs
413	zq(ig,l,igcm_ccn_mass) = 0.
414	zq(ig,l,igcm_ccn_number) = 0.
415
416	endif
417
418	ENDIF !of if Nccn>1
419
420	ENDDO ! of ig loop
421	ENDDO ! of nlayer loop
422
423
424	! Get cloud tendencies
425	subpdqcloud(1:ngrid,1:nlay,igcm_h2o_vap) =
426	& (zq(1:ngrid,1:nlay,igcm_h2o_vap) -
427	& zq0(1:ngrid,1:nlay,igcm_h2o_vap))/microtimestep
428	subpdqcloud(1:ngrid,1:nlay,igcm_h2o_ice) =
429	& (zq(1:ngrid,1:nlay,igcm_h2o_ice) -
430	& zq0(1:ngrid,1:nlay,igcm_h2o_ice))/microtimestep
431	subpdqcloud(1:ngrid,1:nlay,igcm_ccn_mass) =
432	& (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
433	& zq0(1:ngrid,1:nlay,igcm_ccn_mass))/microtimestep
434	subpdqcloud(1:ngrid,1:nlay,igcm_ccn_number) =
435	& (zq(1:ngrid,1:nlay,igcm_ccn_number) -
436	& zq0(1:ngrid,1:nlay,igcm_ccn_number))/microtimestep
437
438	if (scavenging) then
439
440	subpdqcloud(1:ngrid,1:nlay,igcm_dust_mass) =
441	& (zq(1:ngrid,1:nlay,igcm_dust_mass) -
442	& zq0(1:ngrid,1:nlay,igcm_dust_mass))/microtimestep
443	subpdqcloud(1:ngrid,1:nlay,igcm_dust_number) =
444	& (zq(1:ngrid,1:nlay,igcm_dust_number) -
445	& zq0(1:ngrid,1:nlay,igcm_dust_number))/microtimestep
446
447	endif
448
449
450
451	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
452	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
453	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
454	IF (test_flag) then
455
456	! error2d(:) = 0.
457	DO l=1,nlay
458	DO ig=1,ngrid
459	! error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
460	satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l)
461	satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
462	ENDDO
463	ENDDO
464
465	print*, 'count is ',countcells, ' i.e. ',
466	& countcells100/(nlayngrid), '% for microphys computation'
467
468	#ifndef MESOSCALE
469	! IF (ngrid.ne.1) THEN ! 3D
470	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
471	! & satu_out)
472	! call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
473	! & dM_out)
474	! call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
475	! & dN_out)
476	! call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2,
477	! & error2d)
478	! call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3,
479	! & zqsat)
480	! ENDIF
481
482	! IF (ngrid.eq.1) THEN ! 1D
483	! call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1,
484	! & error_out)
485	call WRITEdiagfi(ngrid,"resist","resistance","s/m2",1,
486	& res_out)
487	call WRITEdiagfi(ngrid,"satu_bf","satu before","kg/kg",1,
488	& satubf)
489	call WRITEdiagfi(ngrid,"satu_af","satu after","kg/kg",1,
490	& satuaf)
491	call WRITEdiagfi(ngrid,"vapbf","h2ovap before","kg/kg",1,
492	& zq0(1,1,igcm_h2o_vap))
493	call WRITEdiagfi(ngrid,"vapaf","h2ovap after","kg/kg",1,
494	& zq(1,1,igcm_h2o_vap))
495	call WRITEdiagfi(ngrid,"icebf","h2oice before","kg/kg",1,
496	& zq0(1,1,igcm_h2o_ice))
497	call WRITEdiagfi(ngrid,"iceaf","h2oice after","kg/kg",1,
498	& zq(1,1,igcm_h2o_ice))
499	call WRITEdiagfi(ngrid,"ccnbf","ccn before","/kg",1,
500	& zq0(1,1,igcm_ccn_number))
501	call WRITEdiagfi(ngrid,"ccnaf","ccn after","/kg",1,
502	& zq(1,1,igcm_ccn_number))
503	c call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
504	c & gr_out)
505	c call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
506	c & rate_out)
507	c call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
508	c & dM_out)
509	c call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
510	c & dN_out)
511	call WRITEdiagfi(ngrid,"zqsat","p vap sat","kg/kg",1,
512	& zqsat)
513	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
514	! & satu_out)
515	call WRITEdiagfi(ngrid,"rice","ice radius","m",1,
516	& rice)
517	! call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
518	! & rdust)
519	! call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
520	! & rsedcloud)
521	! call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
522	! & rhocloud)
523	! ENDIF
524	#endif
525
526	ENDIF ! endif test_flag
527	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
528	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
529	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
530
531	return
532
533
534
535
536	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
537	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
538	c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0
539	c It is an analytical solution to the ice radius growth equation,
540	c with the approximation of a constant 'reduced' cunningham correction factor
541	c (lambda in growthrate.F) taken at radius req instead of rice
542	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
543	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
544
545	c subroutine phi(rice,req,coeff1,coeff2,time)
546	c
547	c implicit none
548	c
549	c ! inputs
550	c real rice ! ice radius
551	c real req ! ice radius at equilibirum
552	c real coeff1 ! coeff for the log
553	c real coeff2 ! coeff for the arctan
554	c
555	c ! output
556	c real time
557	c
558	c !local
559	c real var
560	c
561	c ! 1.73205 is sqrt(3)
562	c
563	c var = max(
564	c & abs(rice-req) / sqrt(ricerice + ricereq + req*req),1e-30)
565	c
566	c time =
567	c & coeff1 *
568	c & log( var )
569	c & + coeff2 * 1.73205 *
570	c & atan( (2rice+req) / (1.73205req) )
571	c
572	c return
573	c end
574
575
576
577	END SUBROUTINE improvedclouds
578
579	END MODULE improvedclouds_mod

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

Download in other formats:

Original Format