Context Navigation

improvedclouds.F @ 640

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

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

Download in other formats:

Original Format