Context Navigation

improvedclouds.F @ 1655

Last change on this file since 1655 was 1307, checked in by tnavarro, 10 years ago
Changed two microphysics thresholds for timesteps used in Mesoscale and LES
File size: 20.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	USE updaterad
8	use tracer_mod, only: rho_ice, nuice_sed, igcm_h2o_vap,
9	& igcm_h2o_ice, igcm_dust_mass,
10	& igcm_dust_number, igcm_ccn_mass,
11	& igcm_ccn_number
12	use conc_mod, only: mmean
13	USE comcstfi_h
14	implicit none
15
16
17	c------------------------------------------------------------------
18	c This routine is used to form clouds when a parcel of the GCM is
19	c saturated. It includes the ability to have supersaturation, a
20	c computation of the nucleation rates, growthrates and the
21	c scavenging of dust particles by clouds.
22	c It is worth noting that the amount of dust is computed using the
23	c dust optical depth computed in aeropacity.F. That's why
24	c the variable called "tauscaling" is used to convert
25	c pq(dust_mass) and pq(dust_number), which are relative
26	c quantities, to absolute and realistic quantities stored in zq.
27	c This has to be done to convert the inputs into absolute
28	c values, but also to convert the outputs back into relative
29	c values which are then used by the sedimentation and advection
30	c schemes.
31
32	c Authors: J.-B. Madeleine, based on the work by Franck Montmessin
33	c (October 2011)
34	c T. Navarro, debug,correction, new scheme (October-April 2011)
35	c A. Spiga, optimization (February 2012)
36	c------------------------------------------------------------------
37	#include "callkeys.h"
38	#include "microphys.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(ngrid) ! 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	INTEGER ig,l,i
73
74	REAL zq(ngrid,nlay,nq) ! local value of tracers
75	REAL zq0(ngrid,nlay,nq) ! local initial value of tracers
76	REAL zt(ngrid,nlay) ! local value of temperature
77	REAL zqsat(ngrid,nlay) ! saturation
78	REAL lw !Latent heat of sublimation (J.kg-1)
79	REAL cste
80	REAL dMice ! mass of condensed ice
81	! REAL sumcheck
82	REAL*8 ph2o ! Water vapor partial pressure (Pa)
83	REAL*8 satu ! Water vapor saturation ratio over ice
84	REAL*8 Mo,No
85	REAL*8 Rn, Rm, dev2, n_derf, m_derf
86	REAL*8 n_aer(nbin_cld) ! number conc. of particle/each size bin
87	REAL*8 m_aer(nbin_cld) ! mass mixing ratio of particle/each size bin
88
89	REAL*8 sig ! Water-ice/air surface tension (N.m)
90	EXTERNAL sig
91
92	REAL dN,dM
93	REAL rate(nbin_cld) ! nucleation rate
94	REAL seq
95
96	REAL rice(ngrid,nlay) ! Ice mass mean radius (m)
97	! (r_c in montmessin_2004)
98	REAL rhocloud(ngrid,nlay) ! Cloud density (kg.m-3)
99	REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m)
100
101	REAL res ! Resistance growth
102
103
104	c Parameters of the size discretization
105	c used by the microphysical scheme
106	DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m)
107	DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m)
108	DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.0001e-6
109	! Minimum boundary radius (m)
110	DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m)
111	DOUBLE PRECISION vrat_cld ! Volume ratio
112	DOUBLE PRECISION rb_cld(nbin_cld+1)! boundary values of each rad_cld bin (m)
113	SAVE rb_cld
114	DOUBLE PRECISION dr_cld(nbin_cld) ! width of each rad_cld bin (m)
115	DOUBLE PRECISION vol_cld(nbin_cld) ! particle volume for each bin (m3)
116
117
118	REAL sigma_ice ! Variance of the ice and CCN distributions
119	SAVE sigma_ice
120
121
122
123	c----------------------------------
124	c TESTS
125
126	INTEGER countcells
127
128	LOGICAL test_flag ! flag for test/debuging outputs
129	SAVE test_flag
130
131
132	REAL satubf(ngrid,nlay),satuaf(ngrid,nlay)
133	REAL res_out(ngrid,nlay)
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	res_out(:,:) = 0
220	rice(:,:) = 1.e-8
221
222	c Initialize the tendencies
223	pdqcloud(1:ngrid,1:nlay,1:nq)=0
224	pdtcloud(1:ngrid,1:nlay)=0
225
226
227	zt(1:ngrid,1:nlay) =
228	& pt(1:ngrid,1:nlay) +
229	& pdt(1:ngrid,1:nlay) * ptimestep
230
231	zq(1:ngrid,1:nlay,1:nq) =
232	& pq(1:ngrid,1:nlay,1:nq) +
233	& pdq(1:ngrid,1:nlay,1:nq) * ptimestep
234
235
236	WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
237	& zq(1:ngrid,1:nlay,1:nq) = 1.e-30
238
239	zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
240
241	!=============================================================
242	! 2. Compute saturation
243	!=============================================================
244
245
246	dev2 = 1. / ( sqrt(2.) * sigma_ice )
247
248	call watersat(ngrid*nlay,zt,pplay,zqsat)
249
250	countcells = 0
251
252	c Main loop over the GCM's grid
253	DO l=1,nlay
254	DO ig=1,ngrid
255
256	c Get the partial pressure of water vapor and its saturation ratio
257	ph2o = zq(ig,l,igcm_h2o_vap) * (mmean(ig,l)/18.) * pplay(ig,l)
258	satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
259
260	!=============================================================
261	! 3. Nucleation
262	!=============================================================
263
264	IF ( satu .ge. 1. ) THEN ! if there is condensation
265
266	call updaterccn(zq(ig,l,igcm_dust_mass),
267	& zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
268
269
270	c Expand the dust moments into a binned distribution
271	Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig) + 1.e-30
272	No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30
273	Rn = rdust(ig,l)
274	Rn = -dlog(Rn)
275	Rm = Rn - 3. * sigma_ice*sigma_ice
276	n_derf = derf( (rb_cld(1)+Rn) *dev2)
277	m_derf = derf( (rb_cld(1)+Rm) *dev2)
278	do i = 1, nbin_cld
279	n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
280	m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
281	n_derf = derf( (rb_cld(i+1)+Rn) *dev2)
282	m_derf = derf( (rb_cld(i+1)+Rm) *dev2)
283	n_aer(i) = n_aer(i) + 0.5 * No * n_derf
284	m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
285	enddo
286
287	! sumcheck = 0
288	! do i = 1, nbin_cld
289	! sumcheck = sumcheck + n_aer(i)
290	! enddo
291	! sumcheck = abs(sumcheck/No - 1)
292	! if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
293	! print*, "WARNING, No sumcheck PROBLEM"
294	! print*, "sumcheck, No",sumcheck, No
295	! print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
296	! print*, "Dust binned distribution", n_aer
297	! endif
298	!
299	! sumcheck = 0
300	! do i = 1, nbin_cld
301	! sumcheck = sumcheck + m_aer(i)
302	! enddo
303	! sumcheck = abs(sumcheck/Mo - 1)
304	! if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) then
305	! print*, "WARNING, Mo sumcheck PROBLEM"
306	! print*, "sumcheck, Mo",sumcheck, Mo
307	! print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
308	! print*, "Dust binned distribution", m_aer
309	! endif
310
311
312	c Get the rates of nucleation
313	call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
314
315	dN = 0.
316	dM = 0.
317	do i = 1, nbin_cld
318	n_aer(i) = n_aer(i)/( 1. + rate(i)*ptimestep)
319	m_aer(i) = m_aer(i)/( 1. + rate(i)*ptimestep)
320	dN = dN + n_aer(i) * rate(i) * ptimestep
321	dM = dM + m_aer(i) * rate(i) * ptimestep
322	enddo
323
324
325	c Update Dust particles
326	zq(ig,l,igcm_dust_mass) =
327	& zq(ig,l,igcm_dust_mass) - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
328	zq(ig,l,igcm_dust_number) =
329	& zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
330	c Update CCNs
331	zq(ig,l,igcm_ccn_mass) =
332	& zq(ig,l,igcm_ccn_mass) + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
333	zq(ig,l,igcm_ccn_number) =
334	& zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
335
336	ENDIF ! of is satu >1
337
338	!=============================================================
339	! 4. Ice growth: scheme for radius evolution
340	!=============================================================
341
342	c We trigger crystal growth if and only if there is at least one nuclei (N>1).
343	c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
344	c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
345
346	IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.) THEN ! we trigger crystal growth
347
348
349	call updaterice_micro(zq(ig,l,igcm_h2o_ice),
350	& zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
351	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
352
353	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
354
355	c saturation at equilibrium
356	c rice should not be too small, otherwise seq value is not valid
357	seq = exp(2.sig(zt(ig,l))mh2o / (rho_icergpzt(ig,l)*
358	& max(rice(ig,l),1.e-7)))
359
360	c get resistance growth
361	call growthrate(zt(ig,l),pplay(ig,l),
362	& real(ph2o/satu),rice(ig,l),res)
363
364	res_out(ig,l) = res
365
366	ccccccc implicit scheme of mass growth
367
368	dMice =
369	& (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))
370	& /(reszqsat(ig,l)/(csteNo*rice(ig,l)) + 1.)
371
372
373	! With the above scheme, dMice cannot be bigger than vapor,
374	! but can be bigger than all available ice.
375	dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
376	dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless...
377
378	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice
379	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice
380
381
382	countcells = countcells + 1
383
384	! latent heat release
385	lw=(2834.3-0.28*(zt(ig,l)-To)-
386	& 0.004(zt(ig,l)-To)(zt(ig,l)-To))*1.e+3
387	pdtcloud(ig,l)= dMice*lw/cpp/ptimestep
388
389
390
391	!=============================================================
392	! 5. Dust cores released, tendancies, latent heat, etc ...
393	!=============================================================
394
395	c If all the ice particles sublimate, all the condensation
396	c nuclei are released:
397	if (zq(ig,l,igcm_h2o_ice).le.1.e-28) then
398
399	c Water
400	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)
401	& + zq(ig,l,igcm_h2o_ice)
402	zq(ig,l,igcm_h2o_ice) = 0.
403	c Dust particles
404	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass)
405	& + zq(ig,l,igcm_ccn_mass)
406	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number)
407	& + zq(ig,l,igcm_ccn_number)
408	c CCNs
409	zq(ig,l,igcm_ccn_mass) = 0.
410	zq(ig,l,igcm_ccn_number) = 0.
411
412	endif
413
414	ENDIF !of if Nccn>1
415
416	ENDDO ! of ig loop
417	ENDDO ! of nlayer loop
418
419
420	! Get cloud tendencies
421	pdqcloud(1:ngrid,1:nlay,igcm_h2o_vap) =
422	& (zq(1:ngrid,1:nlay,igcm_h2o_vap) -
423	& zq0(1:ngrid,1:nlay,igcm_h2o_vap))/ptimestep
424	pdqcloud(1:ngrid,1:nlay,igcm_h2o_ice) =
425	& (zq(1:ngrid,1:nlay,igcm_h2o_ice) -
426	& zq0(1:ngrid,1:nlay,igcm_h2o_ice))/ptimestep
427	pdqcloud(1:ngrid,1:nlay,igcm_ccn_mass) =
428	& (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
429	& zq0(1:ngrid,1:nlay,igcm_ccn_mass))/ptimestep
430	pdqcloud(1:ngrid,1:nlay,igcm_ccn_number) =
431	& (zq(1:ngrid,1:nlay,igcm_ccn_number) -
432	& zq0(1:ngrid,1:nlay,igcm_ccn_number))/ptimestep
433
434	if (scavenging) then
435
436	pdqcloud(1:ngrid,1:nlay,igcm_dust_mass) =
437	& (zq(1:ngrid,1:nlay,igcm_dust_mass) -
438	& zq0(1:ngrid,1:nlay,igcm_dust_mass))/ptimestep
439	pdqcloud(1:ngrid,1:nlay,igcm_dust_number) =
440	& (zq(1:ngrid,1:nlay,igcm_dust_number) -
441	& zq0(1:ngrid,1:nlay,igcm_dust_number))/ptimestep
442
443	endif
444
445
446
447	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
448	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
449	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
450	IF (test_flag) then
451
452	! error2d(:) = 0.
453	DO l=1,nlay
454	DO ig=1,ngrid
455	! error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
456	satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l)
457	satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
458	ENDDO
459	ENDDO
460
461	print*, 'count is ',countcells, ' i.e. ',
462	& countcells100/(nlayngrid), '% for microphys computation'
463
464	#ifndef MESOSCALE
465	! IF (ngrid.ne.1) THEN ! 3D
466	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
467	! & satu_out)
468	! call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
469	! & dM_out)
470	! call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
471	! & dN_out)
472	! call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2,
473	! & error2d)
474	! call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3,
475	! & zqsat)
476	! ENDIF
477
478	! IF (ngrid.eq.1) THEN ! 1D
479	! call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1,
480	! & error_out)
481	call WRITEdiagfi(ngrid,"resist","resistance","s/m2",1,
482	& res_out)
483	call WRITEdiagfi(ngrid,"satu_bf","satu before","kg/kg",1,
484	& satubf)
485	call WRITEdiagfi(ngrid,"satu_af","satu after","kg/kg",1,
486	& satuaf)
487	call WRITEdiagfi(ngrid,"vapbf","h2ovap before","kg/kg",1,
488	& zq0(1,1,igcm_h2o_vap))
489	call WRITEdiagfi(ngrid,"vapaf","h2ovap after","kg/kg",1,
490	& zq(1,1,igcm_h2o_vap))
491	call WRITEdiagfi(ngrid,"icebf","h2oice before","kg/kg",1,
492	& zq0(1,1,igcm_h2o_ice))
493	call WRITEdiagfi(ngrid,"iceaf","h2oice after","kg/kg",1,
494	& zq(1,1,igcm_h2o_ice))
495	call WRITEdiagfi(ngrid,"ccnbf","ccn before","/kg",1,
496	& zq0(1,1,igcm_ccn_number))
497	call WRITEdiagfi(ngrid,"ccnaf","ccn after","/kg",1,
498	& zq(1,1,igcm_ccn_number))
499	c call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
500	c & gr_out)
501	c call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
502	c & rate_out)
503	c call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
504	c & dM_out)
505	c call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
506	c & dN_out)
507	call WRITEdiagfi(ngrid,"zqsat","p vap sat","kg/kg",1,
508	& zqsat)
509	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
510	! & satu_out)
511	call WRITEdiagfi(ngrid,"rice","ice radius","m",1,
512	& rice)
513	! call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
514	! & rdust)
515	! call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
516	! & rsedcloud)
517	! call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
518	! & rhocloud)
519	! ENDIF
520	#endif
521
522	ENDIF ! endif test_flag
523	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
524	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
525	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
526
527	return
528	end
529
530
531
532	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
533	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
534	c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0
535	c It is an analytical solution to the ice radius growth equation,
536	c with the approximation of a constant 'reduced' cunningham correction factor
537	c (lambda in growthrate.F) taken at radius req instead of rice
538	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
539	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
540
541	c subroutine phi(rice,req,coeff1,coeff2,time)
542	c
543	c implicit none
544	c
545	c ! inputs
546	c real rice ! ice radius
547	c real req ! ice radius at equilibirum
548	c real coeff1 ! coeff for the log
549	c real coeff2 ! coeff for the arctan
550	c
551	c ! output
552	c real time
553	c
554	c !local
555	c real var
556	c
557	c ! 1.73205 is sqrt(3)
558	c
559	c var = max(
560	c & abs(rice-req) / sqrt(ricerice + ricereq + req*req),1e-30)
561	c
562	c time =
563	c & coeff1 *
564	c & log( var )
565	c & + coeff2 * 1.73205 *
566	c & atan( (2rice+req) / (1.73205req) )
567	c
568	c return
569	c end
570
571
572

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

Download in other formats:

Original Format