Context Navigation

improvedclouds.F @ 1009

Last change on this file since 1009 was 768, checked in by tnavarro, 12 years ago
water conservation bug
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
245	WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
246	& zq(1:ngrid,1:nlay,1:nq) = 1.e-30
247
248	zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
249
250	!=============================================================
251	! 2. Compute saturation
252	!=============================================================
253
254
255	dev2 = 1. / ( sqrt(2.) * sigma_ice )
256
257	call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
258
259	countcells = 0
260
261	c Main loop over the GCM's grid
262	DO l=1,nlay
263	DO ig=1,ngrid
264
265	c Get the partial pressure of water vapor and its saturation ratio
266	ph2o = zq(ig,l,igcm_h2o_vap) * (mmean(ig,l)/18.) * pplay(ig,l)
267	satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
268
269	!=============================================================
270	! 3. Nucleation
271	!=============================================================
272
273	IF ( satu .ge. 1. ) THEN ! if there is condensation
274
275	call updaterccn(zq(ig,l,igcm_dust_mass),
276	& zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
277
278
279	c Expand the dust moments into a binned distribution
280	Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig) + 1.e-30
281	No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30
282	Rn = rdust(ig,l)
283	Rn = -dlog(Rn)
284	Rm = Rn - 3. * sigma_ice*sigma_ice
285	n_derf = derf( (rb_cld(1)+Rn) *dev2)
286	m_derf = derf( (rb_cld(1)+Rm) *dev2)
287	do i = 1, nbin_cld
288	n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
289	m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
290	n_derf = derf( (rb_cld(i+1)+Rn) *dev2)
291	m_derf = derf( (rb_cld(i+1)+Rm) *dev2)
292	n_aer(i) = n_aer(i) + 0.5 * No * n_derf
293	m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
294	enddo
295
296	! sumcheck = 0
297	! do i = 1, nbin_cld
298	! sumcheck = sumcheck + n_aer(i)
299	! enddo
300	! sumcheck = abs(sumcheck/No - 1)
301	! if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
302	! print*, "WARNING, No sumcheck PROBLEM"
303	! print*, "sumcheck, No",sumcheck, No
304	! print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
305	! print*, "Dust binned distribution", n_aer
306	! endif
307	!
308	! sumcheck = 0
309	! do i = 1, nbin_cld
310	! sumcheck = sumcheck + m_aer(i)
311	! enddo
312	! sumcheck = abs(sumcheck/Mo - 1)
313	! if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) then
314	! print*, "WARNING, Mo sumcheck PROBLEM"
315	! print*, "sumcheck, Mo",sumcheck, Mo
316	! print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
317	! print*, "Dust binned distribution", m_aer
318	! endif
319
320
321	c Get the rates of nucleation
322	call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
323
324	dN = 0.
325	dM = 0.
326	do i = 1, nbin_cld
327	n_aer(i) = n_aer(i)/( 1. + rate(i)*ptimestep)
328	m_aer(i) = m_aer(i)/( 1. + rate(i)*ptimestep)
329	dN = dN + n_aer(i) * rate(i) * ptimestep
330	dM = dM + m_aer(i) * rate(i) * ptimestep
331	enddo
332
333
334	c Update Dust particles
335	zq(ig,l,igcm_dust_mass) =
336	& zq(ig,l,igcm_dust_mass) - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
337	zq(ig,l,igcm_dust_number) =
338	& zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
339	c Update CCNs
340	zq(ig,l,igcm_ccn_mass) =
341	& zq(ig,l,igcm_ccn_mass) + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
342	zq(ig,l,igcm_ccn_number) =
343	& zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
344
345	ENDIF ! of is satu >1
346
347	!=============================================================
348	! 4. Ice growth: scheme for radius evolution
349	!=============================================================
350
351	c We trigger crystal growth if and only if there is at least one nuclei (N>1).
352	c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
353	c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
354
355	IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.) THEN ! we trigger crystal growth
356
357
358	call updaterice_micro(zq(ig,l,igcm_h2o_ice),
359	& zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
360	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
361
362	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
363
364	c saturation at equilibrium
365	c rice should not be too small, otherwise seq value is not valid
366	seq = exp(2.sig(zt(ig,l))mh2o / (rho_icergpzt(ig,l)*
367	& max(rice(ig,l),1.e-7)))
368
369	c get resistance growth
370	call growthrate(zt(ig,l),pplay(ig,l),
371	& real(ph2o/satu),rice(ig,l),res)
372
373	res_out(ig,l) = res
374
375	ccccccc implicit scheme of mass growth
376
377	dMice =
378	& (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))
379	& /(reszqsat(ig,l)/(csteNo*rice(ig,l)) + 1.)
380
381
382	! With the above scheme, dMice cannot be bigger than vapor,
383	! but can be bigger than all available ice.
384	dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
385	dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless...
386
387	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice
388	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice
389
390
391	countcells = countcells + 1
392
393	! latent heat release
394	lw=(2834.3-0.28*(zt(ig,l)-To)-
395	& 0.004(zt(ig,l)-To)(zt(ig,l)-To))*1.e+3
396	pdtcloud(ig,l)= dMice*lw/cpp/ptimestep
397
398
399
400	!=============================================================
401	! 5. Dust cores released, tendancies, latent heat, etc ...
402	!=============================================================
403
404	c If all the ice particles sublimate, all the condensation
405	c nuclei are released:
406	if (zq(ig,l,igcm_h2o_ice).le.1.e-8) then
407
408	c Water
409	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)
410	& + zq(ig,l,igcm_h2o_ice)
411	zq(ig,l,igcm_h2o_ice) = 0.
412	c Dust particles
413	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass)
414	& + zq(ig,l,igcm_ccn_mass)
415	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number)
416	& + zq(ig,l,igcm_ccn_number)
417	c CCNs
418	zq(ig,l,igcm_ccn_mass) = 0.
419	zq(ig,l,igcm_ccn_number) = 0.
420
421	endif
422
423	ENDIF !of if Nccn>1
424
425	ENDDO ! of ig loop
426	ENDDO ! of nlayer loop
427
428
429	! Get cloud tendencies
430	pdqcloud(1:ngrid,1:nlay,igcm_h2o_vap) =
431	& (zq(1:ngrid,1:nlay,igcm_h2o_vap) -
432	& zq0(1:ngrid,1:nlay,igcm_h2o_vap))/ptimestep
433	pdqcloud(1:ngrid,1:nlay,igcm_h2o_ice) =
434	& (zq(1:ngrid,1:nlay,igcm_h2o_ice) -
435	& zq0(1:ngrid,1:nlay,igcm_h2o_ice))/ptimestep
436	pdqcloud(1:ngrid,1:nlay,igcm_ccn_mass) =
437	& (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
438	& zq0(1:ngrid,1:nlay,igcm_ccn_mass))/ptimestep
439	pdqcloud(1:ngrid,1:nlay,igcm_ccn_number) =
440	& (zq(1:ngrid,1:nlay,igcm_ccn_number) -
441	& zq0(1:ngrid,1:nlay,igcm_ccn_number))/ptimestep
442
443	if (scavenging) then
444
445	pdqcloud(1:ngrid,1:nlay,igcm_dust_mass) =
446	& (zq(1:ngrid,1:nlay,igcm_dust_mass) -
447	& zq0(1:ngrid,1:nlay,igcm_dust_mass))/ptimestep
448	pdqcloud(1:ngrid,1:nlay,igcm_dust_number) =
449	& (zq(1:ngrid,1:nlay,igcm_dust_number) -
450	& zq0(1:ngrid,1:nlay,igcm_dust_number))/ptimestep
451
452	endif
453
454
455
456	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
457	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
458	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
459	IF (test_flag) then
460
461	! error2d(:) = 0.
462	DO l=1,nlay
463	DO ig=1,ngrid
464	! error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
465	satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l)
466	satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
467	ENDDO
468	ENDDO
469
470	print*, 'count is ',countcells, ' i.e. ',
471	& countcells100/(nlayngrid), '% for microphys computation'
472
473	! IF (ngrid.ne.1) THEN ! 3D
474	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
475	! & satu_out)
476	! call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
477	! & dM_out)
478	! call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
479	! & dN_out)
480	! call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2,
481	! & error2d)
482	! call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3,
483	! & zqsat)
484	! ENDIF
485
486	! IF (ngrid.eq.1) THEN ! 1D
487	! call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1,
488	! & error_out)
489	call WRITEdiagfi(ngrid,"resist","resistance","s/m2",1,
490	& res_out)
491	call WRITEdiagfi(ngrid,"satu_bf","satu before","kg/kg",1,
492	& satubf)
493	call WRITEdiagfi(ngrid,"satu_af","satu after","kg/kg",1,
494	& satuaf)
495	call WRITEdiagfi(ngrid,"vapbf","h2ovap before","kg/kg",1,
496	& zq0(1,:,igcm_h2o_vap))
497	call WRITEdiagfi(ngrid,"vapaf","h2ovap after","kg/kg",1,
498	& zq(1,:,igcm_h2o_vap))
499	call WRITEdiagfi(ngrid,"icebf","h2oice before","kg/kg",1,
500	& zq0(1,:,igcm_h2o_ice))
501	call WRITEdiagfi(ngrid,"iceaf","h2oice after","kg/kg",1,
502	& zq(1,:,igcm_h2o_ice))
503	call WRITEdiagfi(ngrid,"ccnbf","ccn before","/kg",1,
504	& zq0(1,:,igcm_ccn_number))
505	call WRITEdiagfi(ngrid,"ccnaf","ccn after","/kg",1,
506	& zq(1,:,igcm_ccn_number))
507	c call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
508	c & gr_out)
509	c call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
510	c & rate_out)
511	c call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
512	c & dM_out)
513	c call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
514	c & dN_out)
515	call WRITEdiagfi(ngrid,"zqsat","p vap sat","kg/kg",1,
516	& zqsat)
517	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
518	! & satu_out)
519	call WRITEdiagfi(ngrid,"rice","ice radius","m",1,
520	& rice)
521	! call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
522	! & rdust)
523	! call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
524	! & rsedcloud)
525	! call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
526	! & rhocloud)
527	! ENDIF
528
529	ENDIF ! endif test_flag
530	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
531	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
532	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
533
534	return
535	end
536
537
538
539	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
540	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
541	c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0
542	c It is an analytical solution to the ice radius growth equation,
543	c with the approximation of a constant 'reduced' cunningham correction factor
544	c (lambda in growthrate.F) taken at radius req instead of rice
545	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
546	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
547
548	c subroutine phi(rice,req,coeff1,coeff2,time)
549	c
550	c implicit none
551	c
552	c ! inputs
553	c real rice ! ice radius
554	c real req ! ice radius at equilibirum
555	c real coeff1 ! coeff for the log
556	c real coeff2 ! coeff for the arctan
557	c
558	c ! output
559	c real time
560	c
561	c !local
562	c real var
563	c
564	c ! 1.73205 is sqrt(3)
565	c
566	c var = max(
567	c & abs(rice-req) / sqrt(ricerice + ricereq + req*req),1e-30)
568	c
569	c time =
570	c & coeff1 *
571	c & log( var )
572	c & + coeff2 * 1.73205 *
573	c & atan( (2rice+req) / (1.73205req) )
574	c
575	c return
576	c end
577
578
579

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

Download in other formats:

Original Format