Context Navigation

improvedclouds.F @ 706

Last change on this file since 706 was 689, checked in by flefevre, 13 years ago
Minor change. Better handling of very small water mixing ratio at saturation (for very low temperatures).
File size: 20.0 KB

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

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

Download in other formats:

Original Format