Context Navigation

improvedclouds.F @ 756

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

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

Download in other formats:

Original Format