Context Navigation

improvedclouds_mod.F @ 2966

Last change on this file since 2966 was 2966, checked in by jnaar, 19 months ago
Architecture changes in watercloud_mod, improvedclouds_mod in preparation for adaptative timestep of clouds (not yet implemented). Changes prevent bit by bit correspondance with previous revision. "simpleclouds" routine broken.
File size: 22.6 KB

Line
1	MODULE improvedclouds_mod
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	subroutine improvedclouds(microtimestep,
8	& pplay,pt,pq,
9	& subpdqcloud,subpdtcloud,
10	& nq,tauscaling,mmean)
11	USE updaterad, ONLY: updaterice_micro, updaterccn
12	USE watersat_mod, ONLY: watersat
13	use tracer_mod, only: rho_ice, nuice_sed, igcm_h2o_vap,
14	& igcm_h2o_ice, igcm_dust_mass,
15	& igcm_dust_number, igcm_ccn_mass,
16	& igcm_ccn_number,
17	& igcm_hdo_vap,igcm_hdo_ice,
18	& qperemin
19	use comcstfi_h, only: pi, cpp
20	use write_output_mod, only: write_output
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 J. Naar, from global to local (no more ngrid, nlay) to allow
44	c different microphysical timesteps (May 2023)
45	c------------------------------------------------------------------
46	#include "callkeys.h"
47	#include "microphys.h"
48	c------------------------------------------------------------------
49	c Inputs/outputs:
50
51	INTEGER, INTENT(IN) :: nq ! nombre de traceurs
52	REAL, INTENT(IN) :: microtimestep ! pas de temps physique (s)
53	REAL, INTENT(IN) :: pplay ! pression au milieu des couches (Pa)
54	REAL, INTENT(IN) :: pt ! temperature at the middle of the
55	! layers (K)
56	! param.
57	REAL, INTENT(IN) :: pq(nq) ! traceur (kg/kg)
58	! (kg/kg.s-1)
59	REAL, INTENT(IN) :: tauscaling ! Convertion factor for qdust and Ndust
60	REAL, INTENT(IN) :: mmean ! Mean atmospheric mass
61
62	REAL, INTENT(OUT) :: subpdqcloud(nq) ! tendance de la condensation
63	! H2O(kg/kg.s-1)
64	REAL, INTENT(OUT) :: subpdtcloud ! 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	!$OMP THREADPRIVATE(firstcall)
75
76	REAL*8 derf ! Error function
77	!external derf
78
79	INTEGER ig,l,i
80
81	REAL zq(nq) ! local value of tracers
82	REAL zq0(nq) ! local initial value of tracers
83	REAL zt ! local value of temperature
84	REAL zqsat_tmp(1) ! saturation
85	REAL zqsat ! saturation
86	REAL lw !Latent heat of sublimation (J.kg-1)
87	REAL cste
88	REAL dMice ! mass of condensed ice
89	REAL dMice_hdo ! mass of condensed HDO ice
90	REAL alpha ! HDO equilibrium fractionation coefficient (Saturation=1)
91	REAL alpha_c ! HDO real fractionation coefficient
92	! REAL sumcheck
93	REAL*8 ph2o ! Water vapor partial pressure (Pa)
94	REAL*8 satu ! Water vapor saturation ratio over ice
95	REAL*8 Mo,No
96	REAL*8 Rn, Rm, dev2, n_derf, m_derf
97	REAL*8 n_aer(nbin_cld) ! number conc. of particle/each size bin
98	REAL*8 m_aer(nbin_cld) ! mass mixing ratio of particle/each size bin
99
100	REAL*8 sig ! Water-ice/air surface tension (N.m)
101	EXTERNAL sig
102
103	REAL dN,dM
104	REAL rate(nbin_cld) ! nucleation rate
105	REAL seq
106
107	REAL rice ! Ice mass mean radius (m)
108	! (r_c in montmessin_2004)
109	REAL rhocloud ! Cloud density (kg.m-3)
110	REAL rdust ! Dust geometric mean radius (m)
111
112	REAL res ! Resistance growth
113	REAL Dv,Dv_hdo ! Water/HDO vapor diffusion coefficient
114
115
116	c Parameters of the size discretization
117	c used by the microphysical scheme
118	DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m)
119	DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m)
120	DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.0001e-6
121	! Minimum boundary radius (m)
122	DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m)
123	DOUBLE PRECISION vrat_cld ! Volume ratio
124	DOUBLE PRECISION rb_cld(nbin_cld+1)! boundary values of each rad_cld bin (m)
125	SAVE rb_cld
126	DOUBLE PRECISION dr_cld(nbin_cld) ! width of each rad_cld bin (m)
127	DOUBLE PRECISION vol_cld(nbin_cld) ! particle volume for each bin (m3)
128
129	!$OMP THREADPRIVATE(rb_cld)
130
131
132	REAL sigma_ice ! Variance of the ice and CCN distributions
133	SAVE sigma_ice
134
135	!$OMP THREADPRIVATE(sigma_ice)
136
137
138
139	c----------------------------------
140	c TESTS
141
142	INTEGER countcells
143
144	LOGICAL test_flag ! flag for test/debuging outputs
145	SAVE test_flag
146
147	!$OMP THREADPRIVATE(test_flag)
148
149
150	REAL satubf,satuaf
151	REAL res_out
152
153
154	c------------------------------------------------------------------
155
156	! AS: firstcall OK absolute
157	IF (firstcall) THEN
158	!=============================================================
159	! 0. Definition of the size grid
160	!=============================================================
161	c rad_cld is the primary radius grid used for microphysics computation.
162	c The grid spacing is computed assuming a constant volume ratio
163	c between two consecutive bins; i.e. vrat_cld.
164	c vrat_cld is determined from the boundary values of the size grid:
165	c rmin_cld and rmax_cld.
166	c The rb_cld array contains the boundary values of each rad_cld bin.
167	c dr_cld is the width of each rad_cld bin.
168
169	c Volume ratio between two adjacent bins
170	! vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
171	! vrat_cld = exp(vrat_cld)
172	vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
173	vrat_cld = exp(vrat_cld)
174	write(,) "vrat_cld", vrat_cld
175
176	rb_cld(1) = rbmin_cld
177	rad_cld(1) = rmin_cld
178	vol_cld(1) = 4./3. * dble(pi) * rmin_cldrmin_cldrmin_cld
179	! vol_cld(1) = 4./3. * pi * rmin_cldrmin_cldrmin_cld
180
181	do i=1,nbin_cld-1
182	rad_cld(i+1) = rad_cld(i) * vrat_cld**(1./3.)
183	vol_cld(i+1) = vol_cld(i) * vrat_cld
184	enddo
185
186	do i=1,nbin_cld
187	rb_cld(i+1)= ( (2.vrat_cld) / (vrat_cld+1.) )(1./3.)
188	& rad_cld(i)
189	dr_cld(i) = rb_cld(i+1) - rb_cld(i)
190	enddo
191	rb_cld(nbin_cld+1) = rbmax_cld
192	dr_cld(nbin_cld) = rb_cld(nbin_cld+1) - rb_cld(nbin_cld)
193
194	print*, ' '
195	print*,'Microphysics: size bin information:'
196	print*,'i,rb_cld(i), rad_cld(i),dr_cld(i)'
197	print*,'-----------------------------------'
198	do i=1,nbin_cld
199	write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i),
200	& dr_cld(i)
201	enddo
202	write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1)
203	print*,'-----------------------------------'
204
205	do i=1,nbin_cld+1
206	! rb_cld(i) = log(rb_cld(i))
207	rb_cld(i) = log(rb_cld(i)) !! we save that so that it is not computed
208	!! at each timestep and gridpoint
209	enddo
210
211	c Contact parameter of water ice on dust ( m=cos(theta) )
212	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
213	! mteta = 0.95
214	write(,) 'water_param contact parameter:', mteta
215
216	c Volume of a water molecule (m3)
217	vo1 = mh2o / dble(rho_ice)
218	c Variance of the ice and CCN distributions
219	sigma_ice = sqrt(log(1.+nuice_sed))
220
221
222	write(,) 'Variance of ice & CCN distribs :', sigma_ice
223	write(,) 'nuice for sedimentation:', nuice_sed
224	write(,) 'Volume of a water molecule:', vo1
225
226
227	test_flag = .false.
228
229	firstcall=.false.
230	END IF
231
232
233	!=============================================================
234	! 1. Initialisation
235	!=============================================================
236	cste = 4pirho_ice*microtimestep
237
238	res_out = 0
239	rice = 1.e-8
240
241	c Initialize the tendencies
242	subpdqcloud(1:nq)=0
243	subpdtcloud=0
244
245	c temperature and tracers previously incremented here
246	c now done outside, in watercloud_mod.F
247	c zt = pteff + sum_subpdt * microtimestep
248	c zq(1:nq) = pqeff(1:nq) + sum_subpdq(1:nq) * microtimestep
249
250	zq(1:nq)=pq(1:nq)
251	zt=pt
252
253	WHERE( zq(1:nq) < 1.e-30 )
254	& zq(1:nq) = 1.e-30
255
256	zq0(1:nq) = zq(1:nq)
257
258	!=============================================================
259	! 2. Compute saturation
260	!=============================================================
261
262
263	dev2 = 1. / ( sqrt(2.) * sigma_ice )
264
265	call watersat(1,(/zt/),(/pplay/),zqsat_tmp)
266	zqsat=zqsat_tmp(1)
267	countcells = 0
268
269	c Get the partial pressure of water vapor and its saturation ratio
270	ph2o = zq(igcm_h2o_vap) * (mmean/18.) * pplay
271	satu = zq(igcm_h2o_vap) / zqsat
272
273	!=============================================================
274	! 3. Nucleation
275	!=============================================================
276
277	IF ( satu .ge. 1. ) THEN ! if there is condensation
278
279	call updaterccn(zq(igcm_dust_mass),
280	& zq(igcm_dust_number),rdust,tauscaling)
281
282
283	c Expand the dust moments into a binned distribution
284	Mo = zq(igcm_dust_mass)* tauscaling + 1.e-30
285	No = zq(igcm_dust_number)* tauscaling + 1.e-30
286	Rn = rdust
287	Rn = -log(Rn)
288	Rm = Rn - 3. * sigma_ice*sigma_ice
289	n_derf = derf( (rb_cld(1)+Rn) *dev2)
290	m_derf = derf( (rb_cld(1)+Rm) *dev2)
291	do i = 1, nbin_cld
292	n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
293	m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
294	n_derf = derf( (rb_cld(i+1)+Rn) *dev2)
295	m_derf = derf( (rb_cld(i+1)+Rm) *dev2)
296	n_aer(i) = n_aer(i) + 0.5 * No * n_derf
297	m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
298	enddo
299
300	! sumcheck = 0
301	! do i = 1, nbin_cld
302	! sumcheck = sumcheck + n_aer(i)
303	! enddo
304	! sumcheck = abs(sumcheck/No - 1)
305	! if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
306	! print*, "WARNING, No sumcheck PROBLEM"
307	! print*, "sumcheck, No",sumcheck, No
308	! print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
309	! print*, "Dust binned distribution", n_aer
310	! endif
311	!
312	! sumcheck = 0
313	! do i = 1, nbin_cld
314	! sumcheck = sumcheck + m_aer(i)
315	! enddo
316	! sumcheck = abs(sumcheck/Mo - 1)
317	! if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) then
318	! print*, "WARNING, Mo sumcheck PROBLEM"
319	! print*, "sumcheck, Mo",sumcheck, Mo
320	! print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
321	! print*, "Dust binned distribution", m_aer
322	! endif
323
324
325	c Get the rates of nucleation
326	call nuclea(ph2o,zt,satu,n_aer,rate)
327
328	dN = 0.
329	dM = 0.
330	do i = 1, nbin_cld
331	dN = dN + n_aer(i)(exp(-rate(i)microtimestep)-1.)
332	dM = dM + m_aer(i)(exp(-rate(i)microtimestep)-1.)
333	enddo
334
335
336	c Update Dust particles
337	zq(igcm_dust_mass) =
338	& zq(igcm_dust_mass) + dM/ tauscaling !max(tauscaling,1.e-10)
339	zq(igcm_dust_number) =
340	& zq(igcm_dust_number) + dN/ tauscaling !max(tauscaling,1.e-10)
341	c Update CCNs
342	zq(igcm_ccn_mass) =
343	& zq(igcm_ccn_mass) - dM/ tauscaling !max(tauscaling,1.e-10)
344	zq(igcm_ccn_number) =
345	& zq(igcm_ccn_number) - dN/ tauscaling !max(tauscaling,1.e-10)
346
347	ENDIF ! of is satu >1
348
349	!=============================================================
350	! 4. Ice growth: scheme for radius evolution
351	!=============================================================
352
353	c We trigger crystal growth if and only if there is at least one nuclei (N>1).
354	c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
355	c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
356
357	IF ( zq(igcm_ccn_number)*tauscaling.ge. 1.) THEN ! we trigger crystal growth
358
359
360	call updaterice_micro(zq(igcm_h2o_ice),
361	& zq(igcm_ccn_mass),zq(igcm_ccn_number),
362	& tauscaling,rice,rhocloud)
363
364	No = zq(igcm_ccn_number)* tauscaling + 1.e-30
365
366	c saturation at equilibrium
367	c rice should not be too small, otherwise seq value is not valid
368	seq = exp(2.sig(zt)mh2o / (rho_icergpzt*
369	& max(rice,1.e-7)))
370
371	c get resistance growth
372	call growthrate(zt,pplay,
373	& real(ph2o/satu),rice,res,Dv)
374
375	res_out = res
376
377	ccccccc implicit scheme of mass growth
378
379	dMice =
380	& (zq(igcm_h2o_vap)-seq*zqsat)
381	& /(reszqsat/(csteNo*rice) + 1.)
382
383
384	! With the above scheme, dMice cannot be bigger than vapor,
385	! but can be bigger than all available ice.
386	dMice = max(dMice,-zq(igcm_h2o_ice))
387	dMice = min(dMice,zq(igcm_h2o_vap)) ! this should be useless...
388
389	zq(igcm_h2o_ice) = zq(igcm_h2o_ice)+dMice
390	zq(igcm_h2o_vap) = zq(igcm_h2o_vap)-dMice
391
392
393	countcells = countcells + 1
394
395	! latent heat release
396	lw=(2834.3-0.28*(zt-To)-
397	& 0.004(zt-To)(zt-To))*1.e+3
398	subpdtcloud= dMice*lw/cpp/microtimestep
399
400	c Special case of the isotope of water HDO
401	if (hdo) then
402	!! condensation
403	if (dMice.gt.0.0) then
404	!! do we use fractionation?
405	if (hdofrac) then
406	!! Calculation of the HDO vapor coefficient
407	Dv_hdo = 1./3. * sqrt( 8kbzzt/(pi*mhdo/nav) )
408	& * kbz * zt /
409	& ( pi * pplay * (molco2+molhdo)*(molco2+molhdo)
410	& * sqrt(1.+mhdo/mco2) )
411	!! Calculation of the fractionnation coefficient at equilibrium
412	c alpha = exp(16288./zt**2.-9.34d-2) ! Merlivat and Nief et al. 1967
413	alpha = exp(13525./zt**2.-5.59d-2) ! Lamb et al. 2017
414	!! Calculation of the 'real' fractionnation coefficient
415	alpha_c = (alpha*satu)/
416	& ( (alpha(Dv/Dv_hdo)(satu-1.)) + 1.)
417	c alpha_c = alpha ! to test without the effect of cinetics
418	else
419	alpha_c = 1.d0
420	endif
421	if (zq0(igcm_h2o_vap).gt.qperemin) then
422	dMice_hdo=
423	& dMicealpha_c
424	& ( zq0(igcm_hdo_vap)
425	& /zq0(igcm_h2o_vap) )
426	else
427	dMice_hdo=0.
428	endif
429	!! sublimation
430	else
431	if (zq0(igcm_h2o_ice).gt.qperemin) then
432	dMice_hdo=
433	& dMice*
434	& ( zq0(igcm_hdo_ice)
435	& /zq0(igcm_h2o_ice) )
436	else
437	dMice_hdo=0.
438	endif
439	endif !if (dMice.gt.0.0)
440
441	dMice_hdo = max(dMice_hdo,-zq(igcm_hdo_ice))
442	dMice_hdo = min(dMice_hdo,zq(igcm_hdo_vap))
443
444	zq(igcm_hdo_ice) = zq(igcm_hdo_ice)+dMice_hdo
445	zq(igcm_hdo_vap) = zq(igcm_hdo_vap)-dMice_hdo
446
447	endif ! if (hdo)
448
449	!=============================================================
450	! 5. Dust cores released, tendancies, latent heat, etc ...
451	!=============================================================
452
453	c If all the ice particles sublimate, all the condensation
454	c nuclei are released:
455	if (zq(igcm_h2o_ice).le.1.e-28) then
456
457	c Water
458	zq(igcm_h2o_vap) = zq(igcm_h2o_vap)
459	& + zq(igcm_h2o_ice)
460	zq(igcm_h2o_ice) = 0.
461	if (hdo) then
462	zq(igcm_hdo_vap) = zq(igcm_hdo_vap)
463	& + zq(igcm_hdo_ice)
464	zq(igcm_hdo_ice) = 0.
465	endif
466	c Dust particles
467	zq(igcm_dust_mass) = zq(igcm_dust_mass)
468	& + zq(igcm_ccn_mass)
469	zq(igcm_dust_number) = zq(igcm_dust_number)
470	& + zq(igcm_ccn_number)
471	c CCNs
472	zq(igcm_ccn_mass) = 0.
473	zq(igcm_ccn_number) = 0.
474
475	endif
476
477	ENDIF !of if Nccn>1
478
479
480	! Get cloud tendencies
481	subpdqcloud(igcm_h2o_vap) =
482	& (zq(igcm_h2o_vap) -
483	& zq0(igcm_h2o_vap))/microtimestep
484	subpdqcloud(igcm_h2o_ice) =
485	& (zq(igcm_h2o_ice) -
486	& zq0(igcm_h2o_ice))/microtimestep
487	if (hdo) then
488	subpdqcloud(igcm_hdo_vap) =
489	& (zq(igcm_hdo_vap) -
490	& zq0(igcm_hdo_vap))/microtimestep
491	subpdqcloud(igcm_hdo_ice) =
492	& (zq(igcm_hdo_ice) -
493	& zq0(igcm_hdo_ice))/microtimestep
494	endif
495	subpdqcloud(igcm_ccn_mass) =
496	& (zq(igcm_ccn_mass) -
497	& zq0(igcm_ccn_mass))/microtimestep
498	subpdqcloud(igcm_ccn_number) =
499	& (zq(igcm_ccn_number) -
500	& zq0(igcm_ccn_number))/microtimestep
501
502	if (scavenging) then
503
504	subpdqcloud(igcm_dust_mass) =
505	& (zq(igcm_dust_mass) -
506	& zq0(igcm_dust_mass))/microtimestep
507	subpdqcloud(igcm_dust_number) =
508	& (zq(igcm_dust_number) -
509	& zq0(igcm_dust_number))/microtimestep
510
511	endif
512
513
514
515	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
516	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
517	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
518	IF (test_flag) then
519
520	! error2d(:) = 0.
521	! error2d(ig) = max(abs(error_out),error2d(ig))
522	satubf = zq0(igcm_h2o_vap)/zqsat
523	satuaf = zq(igcm_h2o_vap)/zqsat
524
525	c print*, 'count is ',countcells, ' i.e. ',
526	c & countcells100/(nlayngrid), '% for microphys computation'
527
528	#ifndef MESOSCALE
529	! IF (ngrid.ne.1) THEN ! 3D
530	! call write_output("satu","ratio saturation","",
531	! & satu_out)
532	! call write_output("dM","ccn variation","kg/kg",
533	! & dM_out)
534	! call write_output("dN","ccn variation","#",
535	! & dN_out)
536	! call write_output("error","dichotomy max error","%",
537	! & error2d)
538	! call write_output("zqsat","zqsat","kg",
539	! & zqsat)
540	! ENDIF
541
542	! IF (ngrid.eq.1) THEN ! 1D
543	! call write_output("error","incertitude sur glace","%",
544	! & error_out)
545	call write_output("resist","resistance","s/m2",
546	& res_out)
547	call write_output("satu_bf","satu before","kg/kg",
548	& satubf)
549	call write_output("satu_af","satu after","kg/kg",
550	& satuaf)
551	call write_output("vapbf","h2ovap before","kg/kg",
552	& zq0(igcm_h2o_vap))
553	call write_output("vapaf","h2ovap after","kg/kg",
554	& zq(igcm_h2o_vap))
555	call write_output("icebf","h2oice before","kg/kg",
556	& zq0(igcm_h2o_ice))
557	call write_output("iceaf","h2oice after","kg/kg",
558	& zq(igcm_h2o_ice))
559	call write_output("ccnbf","ccn before","/kg",
560	& zq0(igcm_ccn_number))
561	call write_output("ccnaf","ccn after","/kg",
562	& zq(igcm_ccn_number))
563	c call write_output("growthrate","growth rate","m^2/s",
564	c & gr_out)
565	c call write_output("nuclearate","nucleation rate","",
566	c & rate_out)
567	c call write_output("dM","ccn variation","kg",
568	c & dM_out)
569	c call write_output("dN","ccn variation","#",
570	c & dN_out)
571	call write_output("zqsat","p vap sat","kg/kg",
572	& zqsat)
573	! call write_output("satu","ratio saturation","",
574	! & satu_out(:,:))
575	call write_output("rice","ice radius","m",
576	& rice)
577	! call write_output("rdust_sca","rdust","m",
578	! & rdust)
579	! call write_output("rsedcloud","rsedcloud","m",
580	! & rsedcloud)
581	! call write_output("rhocloud","rhocloud","kg.m-3",
582	! & rhocloud)
583	! ENDIF
584	#endif
585
586	ENDIF ! endif test_flag
587	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
588	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
589	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
590
591	return
592
593
594
595
596	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
597	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
598	c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0
599	c It is an analytical solution to the ice radius growth equation,
600	c with the approximation of a constant 'reduced' cunningham correction factor
601	c (lambda in growthrate.F) taken at radius req instead of rice
602	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
603	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
604
605	c subroutine phi(rice,req,coeff1,coeff2,time)
606	c
607	c implicit none
608	c
609	c ! inputs
610	c real rice ! ice radius
611	c real req ! ice radius at equilibirum
612	c real coeff1 ! coeff for the log
613	c real coeff2 ! coeff for the arctan
614	c
615	c ! output
616	c real time
617	c
618	c !local
619	c real var
620	c
621	c ! 1.73205 is sqrt(3)
622	c
623	c var = max(
624	c & abs(rice-req) / sqrt(ricerice + ricereq + req*req),1e-30)
625	c
626	c time =
627	c & coeff1 *
628	c & log( var )
629	c & + coeff2 * 1.73205 *
630	c & atan( (2rice+req) / (1.73205req) )
631	c
632	c return
633	c end
634
635
636
637	END SUBROUTINE improvedclouds
638
639	END MODULE improvedclouds_mod

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

Download in other formats:

Original Format