Context Navigation

improvedclouds.F @ 1922

Last change on this file since 1922 was 1909, checked in by mvals, 7 years ago

File size: 21.0 KB

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

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

Download in other formats:

Original Format