Context Navigation

improvedclouds.F @ 1233

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

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

Download in other formats:

Original Format