Context Navigation

improvedclouds_mod.F @ 2156

Last change on this file since 2156 was 2151, checked in by aslmd, 5 years ago
changed old functions dexp dlog in generic exp and log
File size: 21.5 KB

Line
1	MODULE improvedclouds_mod
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	subroutine improvedclouds(ngrid,nlay,microtimestep,
8	& pplay,pteff,sum_subpdt,
9	& pqeff,sum_subpdq,subpdqcloud,subpdtcloud,
10	& nq,tauscaling)
11	! to use 'getin'
12	USE ioipsl_getincom
13	USE updaterad
14	USE watersat_mod, ONLY: watersat
15	use tracer_mod, only: rho_ice, nuice_sed, igcm_h2o_vap,
16	& igcm_h2o_ice, igcm_dust_mass,
17	& igcm_dust_number, igcm_ccn_mass,
18	& igcm_ccn_number
19	use conc_mod, only: mmean
20	USE comcstfi_h
21	implicit none
22
23
24	c------------------------------------------------------------------
25	c This routine is used to form clouds when a parcel of the GCM is
26	c saturated. It includes the ability to have supersaturation, a
27	c computation of the nucleation rates, growthrates and the
28	c scavenging of dust particles by clouds.
29	c It is worth noting that the amount of dust is computed using the
30	c dust optical depth computed in aeropacity.F. That's why
31	c the variable called "tauscaling" is used to convert
32	c pq(dust_mass) and pq(dust_number), which are relative
33	c quantities, to absolute and realistic quantities stored in zq.
34	c This has to be done to convert the inputs into absolute
35	c values, but also to convert the outputs back into relative
36	c values which are then used by the sedimentation and advection
37	c schemes.
38
39	c Authors: J.-B. Madeleine, based on the work by Franck Montmessin
40	c (October 2011)
41	c T. Navarro, debug,correction, new scheme (October-April 2011)
42	c A. Spiga, optimization (February 2012)
43	c------------------------------------------------------------------
44	#include "callkeys.h"
45	#include "microphys.h"
46	c------------------------------------------------------------------
47	c Inputs/outputs:
48
49	INTEGER, INTENT(IN) :: ngrid,nlay
50	INTEGER, INTENT(IN) :: nq ! nombre de traceurs
51	REAL, INTENT(IN) :: microtimestep ! pas de temps physique (s)
52	REAL, INTENT(IN) :: pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
53	REAL, INTENT(IN) :: pteff(ngrid,nlay) ! temperature at the middle of the
54	! layers (K)
55	REAL, INTENT(IN) :: sum_subpdt(ngrid,nlay)! tendance temperature des autres
56	! param.
57	REAL, INTENT(IN) :: pqeff(ngrid,nlay,nq) ! traceur (kg/kg)
58	REAL, INTENT(IN) :: sum_subpdq(ngrid,nlay,nq) ! tendance avant condensation
59	! (kg/kg.s-1)
60	REAL, INTENT(IN) :: tauscaling(ngrid) ! Convertion factor for qdust and Ndust
61
62	REAL, INTENT(OUT) :: subpdqcloud(ngrid,nlay,nq) ! tendance de la condensation
63	! H2O(kg/kg.s-1)
64	REAL, INTENT(OUT) :: subpdtcloud(ngrid,nlay) ! tendance temperature due
65	! a la chaleur latente
66
67	c------------------------------------------------------------------
68	c Local variables:
69
70	LOGICAL firstcall
71	DATA firstcall/.true./
72	SAVE firstcall
73
74	REAL*8 derf ! Error function
75	!external derf
76
77	INTEGER ig,l,i
78
79	REAL zq(ngrid,nlay,nq) ! local value of tracers
80	REAL zq0(ngrid,nlay,nq) ! local initial value of tracers
81	REAL zt(ngrid,nlay) ! local value of temperature
82	REAL zqsat(ngrid,nlay) ! saturation
83	REAL lw !Latent heat of sublimation (J.kg-1)
84	REAL cste
85	REAL dMice ! mass of condensed ice
86	! REAL sumcheck
87	REAL*8 ph2o ! Water vapor partial pressure (Pa)
88	REAL*8 satu ! Water vapor saturation ratio over ice
89	REAL*8 Mo,No
90	REAL*8 Rn, Rm, dev2, n_derf, m_derf
91	REAL*8 n_aer(nbin_cld) ! number conc. of particle/each size bin
92	REAL*8 m_aer(nbin_cld) ! mass mixing ratio of particle/each size bin
93
94	REAL*8 sig ! Water-ice/air surface tension (N.m)
95	EXTERNAL sig
96
97	REAL dN,dM
98	REAL rate(nbin_cld) ! nucleation rate
99	REAL seq
100
101	REAL rice(ngrid,nlay) ! Ice mass mean radius (m)
102	! (r_c in montmessin_2004)
103	REAL rhocloud(ngrid,nlay) ! Cloud density (kg.m-3)
104	REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m)
105
106	REAL res ! Resistance growth
107
108
109	c Parameters of the size discretization
110	c used by the microphysical scheme
111	DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m)
112	DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m)
113	DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.0001e-6
114	! Minimum boundary radius (m)
115	DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m)
116	DOUBLE PRECISION vrat_cld ! Volume ratio
117	DOUBLE PRECISION rb_cld(nbin_cld+1)! boundary values of each rad_cld bin (m)
118	SAVE rb_cld
119	DOUBLE PRECISION dr_cld(nbin_cld) ! width of each rad_cld bin (m)
120	DOUBLE PRECISION vol_cld(nbin_cld) ! particle volume for each bin (m3)
121
122
123	REAL sigma_ice ! Variance of the ice and CCN distributions
124	SAVE sigma_ice
125
126
127
128	c----------------------------------
129	c TESTS
130
131	INTEGER countcells
132
133	LOGICAL test_flag ! flag for test/debuging outputs
134	SAVE test_flag
135
136
137	REAL satubf(ngrid,nlay),satuaf(ngrid,nlay)
138	REAL res_out(ngrid,nlay)
139
140
141	c------------------------------------------------------------------
142
143	! AS: firstcall OK absolute
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 = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
160	vrat_cld = exp(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) = log(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*microtimestep
224
225	res_out(:,:) = 0
226	rice(:,:) = 1.e-8
227
228	c Initialize the tendencies
229	subpdqcloud(1:ngrid,1:nlay,1:nq)=0
230	subpdtcloud(1:ngrid,1:nlay)=0
231
232
233	zt(1:ngrid,1:nlay) =
234	& pteff(1:ngrid,1:nlay) +
235	& sum_subpdt(1:ngrid,1:nlay) * microtimestep
236
237	zq(1:ngrid,1:nlay,1:nq) =
238	& pqeff(1:ngrid,1:nlay,1:nq) +
239	& sum_subpdq(1:ngrid,1:nlay,1:nq) * microtimestep
240
241
242	WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
243	& zq(1:ngrid,1:nlay,1:nq) = 1.e-30
244
245	zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
246
247	!=============================================================
248	! 2. Compute saturation
249	!=============================================================
250
251
252	dev2 = 1. / ( sqrt(2.) * sigma_ice )
253
254	call watersat(ngrid*nlay,zt,pplay,zqsat)
255
256	countcells = 0
257
258	c Main loop over the GCM's grid
259	DO l=1,nlay
260	DO ig=1,ngrid
261
262	c Get the partial pressure of water vapor and its saturation ratio
263	ph2o = zq(ig,l,igcm_h2o_vap) * (mmean(ig,l)/18.) * pplay(ig,l)
264	satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
265
266	!=============================================================
267	! 3. Nucleation
268	!=============================================================
269
270	IF ( satu .ge. 1. ) THEN ! if there is condensation
271
272	call updaterccn(zq(ig,l,igcm_dust_mass),
273	& zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
274
275
276	c Expand the dust moments into a binned distribution
277	Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig) + 1.e-30
278	No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30
279	Rn = rdust(ig,l)
280	Rn = -log(Rn)
281	Rm = Rn - 3. * sigma_ice*sigma_ice
282	n_derf = derf( (rb_cld(1)+Rn) *dev2)
283	m_derf = derf( (rb_cld(1)+Rm) *dev2)
284	do i = 1, nbin_cld
285	n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
286	m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
287	n_derf = derf( (rb_cld(i+1)+Rn) *dev2)
288	m_derf = derf( (rb_cld(i+1)+Rm) *dev2)
289	n_aer(i) = n_aer(i) + 0.5 * No * n_derf
290	m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
291	enddo
292
293	! sumcheck = 0
294	! do i = 1, nbin_cld
295	! sumcheck = sumcheck + n_aer(i)
296	! enddo
297	! sumcheck = abs(sumcheck/No - 1)
298	! if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
299	! print*, "WARNING, No sumcheck PROBLEM"
300	! print*, "sumcheck, No",sumcheck, No
301	! print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
302	! print*, "Dust binned distribution", n_aer
303	! endif
304	!
305	! sumcheck = 0
306	! do i = 1, nbin_cld
307	! sumcheck = sumcheck + m_aer(i)
308	! enddo
309	! sumcheck = abs(sumcheck/Mo - 1)
310	! if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) then
311	! print*, "WARNING, Mo sumcheck PROBLEM"
312	! print*, "sumcheck, Mo",sumcheck, Mo
313	! print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
314	! print*, "Dust binned distribution", m_aer
315	! endif
316
317
318	c Get the rates of nucleation
319	call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
320
321	dN = 0.
322	dM = 0.
323	do i = 1, nbin_cld
324	n_aer(i) = n_aer(i)/( 1. + rate(i)*microtimestep)
325	m_aer(i) = m_aer(i)/( 1. + rate(i)*microtimestep)
326	dN = dN + n_aer(i) * rate(i) * microtimestep
327	dM = dM + m_aer(i) * rate(i) * microtimestep
328	enddo
329
330
331	c Update Dust particles
332	zq(ig,l,igcm_dust_mass) =
333	& zq(ig,l,igcm_dust_mass) - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
334	zq(ig,l,igcm_dust_number) =
335	& zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
336	c Update CCNs
337	zq(ig,l,igcm_ccn_mass) =
338	& zq(ig,l,igcm_ccn_mass) + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
339	zq(ig,l,igcm_ccn_number) =
340	& zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
341
342	ENDIF ! of is satu >1
343
344	!=============================================================
345	! 4. Ice growth: scheme for radius evolution
346	!=============================================================
347
348	c We trigger crystal growth if and only if there is at least one nuclei (N>1).
349	c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
350	c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
351
352	IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.) THEN ! we trigger crystal growth
353
354
355	call updaterice_micro(zq(ig,l,igcm_h2o_ice),
356	& zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
357	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
358
359	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
360
361	c saturation at equilibrium
362	c rice should not be too small, otherwise seq value is not valid
363	seq = exp(2.sig(zt(ig,l))mh2o / (rho_icergpzt(ig,l)*
364	& max(rice(ig,l),1.e-7)))
365
366	c get resistance growth
367	call growthrate(zt(ig,l),pplay(ig,l),
368	& real(ph2o/satu),rice(ig,l),res)
369
370	res_out(ig,l) = res
371
372	ccccccc implicit scheme of mass growth
373
374	dMice =
375	& (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))
376	& /(reszqsat(ig,l)/(csteNo*rice(ig,l)) + 1.)
377
378
379	! With the above scheme, dMice cannot be bigger than vapor,
380	! but can be bigger than all available ice.
381	dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
382	dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless...
383
384	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice
385	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice
386
387
388	countcells = countcells + 1
389
390	! latent heat release
391	lw=(2834.3-0.28*(zt(ig,l)-To)-
392	& 0.004(zt(ig,l)-To)(zt(ig,l)-To))*1.e+3
393	subpdtcloud(ig,l)= dMice*lw/cpp/microtimestep
394
395
396
397	!=============================================================
398	! 5. Dust cores released, tendancies, latent heat, etc ...
399	!=============================================================
400
401	c If all the ice particles sublimate, all the condensation
402	c nuclei are released:
403	if (zq(ig,l,igcm_h2o_ice).le.1.e-28) then
404
405	c Water
406	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)
407	& + zq(ig,l,igcm_h2o_ice)
408	zq(ig,l,igcm_h2o_ice) = 0.
409	c Dust particles
410	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass)
411	& + zq(ig,l,igcm_ccn_mass)
412	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number)
413	& + zq(ig,l,igcm_ccn_number)
414	c CCNs
415	zq(ig,l,igcm_ccn_mass) = 0.
416	zq(ig,l,igcm_ccn_number) = 0.
417
418	endif
419
420	ENDIF !of if Nccn>1
421
422	ENDDO ! of ig loop
423	ENDDO ! of nlayer loop
424
425
426	! Get cloud tendencies
427	subpdqcloud(1:ngrid,1:nlay,igcm_h2o_vap) =
428	& (zq(1:ngrid,1:nlay,igcm_h2o_vap) -
429	& zq0(1:ngrid,1:nlay,igcm_h2o_vap))/microtimestep
430	subpdqcloud(1:ngrid,1:nlay,igcm_h2o_ice) =
431	& (zq(1:ngrid,1:nlay,igcm_h2o_ice) -
432	& zq0(1:ngrid,1:nlay,igcm_h2o_ice))/microtimestep
433	subpdqcloud(1:ngrid,1:nlay,igcm_ccn_mass) =
434	& (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
435	& zq0(1:ngrid,1:nlay,igcm_ccn_mass))/microtimestep
436	subpdqcloud(1:ngrid,1:nlay,igcm_ccn_number) =
437	& (zq(1:ngrid,1:nlay,igcm_ccn_number) -
438	& zq0(1:ngrid,1:nlay,igcm_ccn_number))/microtimestep
439
440	if (scavenging) then
441
442	subpdqcloud(1:ngrid,1:nlay,igcm_dust_mass) =
443	& (zq(1:ngrid,1:nlay,igcm_dust_mass) -
444	& zq0(1:ngrid,1:nlay,igcm_dust_mass))/microtimestep
445	subpdqcloud(1:ngrid,1:nlay,igcm_dust_number) =
446	& (zq(1:ngrid,1:nlay,igcm_dust_number) -
447	& zq0(1:ngrid,1:nlay,igcm_dust_number))/microtimestep
448
449	endif
450
451
452
453	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
454	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
455	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
456	IF (test_flag) then
457
458	! error2d(:) = 0.
459	DO l=1,nlay
460	DO ig=1,ngrid
461	! error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
462	satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l)
463	satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
464	ENDDO
465	ENDDO
466
467	print*, 'count is ',countcells, ' i.e. ',
468	& countcells100/(nlayngrid), '% for microphys computation'
469
470	#ifndef MESOSCALE
471	! IF (ngrid.ne.1) THEN ! 3D
472	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
473	! & satu_out)
474	! call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
475	! & dM_out)
476	! call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
477	! & dN_out)
478	! call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2,
479	! & error2d)
480	! call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3,
481	! & zqsat)
482	! ENDIF
483
484	! IF (ngrid.eq.1) THEN ! 1D
485	! call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1,
486	! & error_out)
487	call WRITEdiagfi(ngrid,"resist","resistance","s/m2",1,
488	& res_out)
489	call WRITEdiagfi(ngrid,"satu_bf","satu before","kg/kg",1,
490	& satubf)
491	call WRITEdiagfi(ngrid,"satu_af","satu after","kg/kg",1,
492	& satuaf)
493	call WRITEdiagfi(ngrid,"vapbf","h2ovap before","kg/kg",1,
494	& zq0(1,1,igcm_h2o_vap))
495	call WRITEdiagfi(ngrid,"vapaf","h2ovap after","kg/kg",1,
496	& zq(1,1,igcm_h2o_vap))
497	call WRITEdiagfi(ngrid,"icebf","h2oice before","kg/kg",1,
498	& zq0(1,1,igcm_h2o_ice))
499	call WRITEdiagfi(ngrid,"iceaf","h2oice after","kg/kg",1,
500	& zq(1,1,igcm_h2o_ice))
501	call WRITEdiagfi(ngrid,"ccnbf","ccn before","/kg",1,
502	& zq0(1,1,igcm_ccn_number))
503	call WRITEdiagfi(ngrid,"ccnaf","ccn after","/kg",1,
504	& zq(1,1,igcm_ccn_number))
505	c call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
506	c & gr_out)
507	c call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
508	c & rate_out)
509	c call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
510	c & dM_out)
511	c call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
512	c & dN_out)
513	call WRITEdiagfi(ngrid,"zqsat","p vap sat","kg/kg",1,
514	& zqsat)
515	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
516	! & satu_out)
517	call WRITEdiagfi(ngrid,"rice","ice radius","m",1,
518	& rice)
519	! call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
520	! & rdust)
521	! call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
522	! & rsedcloud)
523	! call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
524	! & rhocloud)
525	! ENDIF
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
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
579	END SUBROUTINE improvedclouds
580
581	END MODULE improvedclouds_mod

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

Download in other formats:

Original Format