Context Navigation

improvedclouds_mod.F @ 1963

Last change on this file since 1963 was 1963, checked in by mvals, 6 years ago

Mars GCM:

improvedclouds.F changed to module improvedclouds_mod.F

File size: 21.2 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 tracer_mod, only: rho_ice, nuice_sed, igcm_h2o_vap,
15	& igcm_h2o_ice, igcm_dust_mass,
16	& igcm_dust_number, igcm_ccn_mass,
17	& igcm_ccn_number
18	use conc_mod, only: mmean
19	USE comcstfi_h
20	implicit none
21
22
23	c------------------------------------------------------------------
24	c This routine is used to form clouds when a parcel of the GCM is
25	c saturated. It includes the ability to have supersaturation, a
26	c computation of the nucleation rates, growthrates and the
27	c scavenging of dust particles by clouds.
28	c It is worth noting that the amount of dust is computed using the
29	c dust optical depth computed in aeropacity.F. That's why
30	c the variable called "tauscaling" is used to convert
31	c pq(dust_mass) and pq(dust_number), which are relative
32	c quantities, to absolute and realistic quantities stored in zq.
33	c This has to be done to convert the inputs into absolute
34	c values, but also to convert the outputs back into relative
35	c values which are then used by the sedimentation and advection
36	c schemes.
37
38	c Authors: J.-B. Madeleine, based on the work by Franck Montmessin
39	c (October 2011)
40	c T. Navarro, debug,correction, new scheme (October-April 2011)
41	c A. Spiga, optimization (February 2012)
42	c------------------------------------------------------------------
43	#include "callkeys.h"
44	#include "microphys.h"
45	c------------------------------------------------------------------
46	c Inputs:
47
48	INTEGER ngrid,nlay
49	integer nq ! nombre de traceurs
50	REAL microtimestep ! pas de temps physique (s)
51	REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
52
53	REAL pteff(ngrid,nlay) ! temperature at the middle of the
54	! layers (K)
55	REAL sum_subpdt(ngrid,nlay) ! tendance temperature des autres
56	! param.
57	REAL pqeff(ngrid,nlay,nq) ! traceur (kg/kg)
58	REAL sum_subpdq(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 subpdqcloud(ngrid,nlay,nq) ! tendance de la condensation
64	! H2O(kg/kg.s-1)
65	REAL subpdtcloud(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	! AS: firstcall OK absolute
145	IF (firstcall) THEN
146	!=============================================================
147	! 0. Definition of the size grid
148	!=============================================================
149	c rad_cld is the primary radius grid used for microphysics computation.
150	c The grid spacing is computed assuming a constant volume ratio
151	c between two consecutive bins; i.e. vrat_cld.
152	c vrat_cld is determined from the boundary values of the size grid:
153	c rmin_cld and rmax_cld.
154	c The rb_cld array contains the boundary values of each rad_cld bin.
155	c dr_cld is the width of each rad_cld bin.
156
157	c Volume ratio between two adjacent bins
158	! vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
159	! vrat_cld = exp(vrat_cld)
160	vrat_cld = dlog(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
161	vrat_cld = dexp(vrat_cld)
162	write(,) "vrat_cld", vrat_cld
163
164	rb_cld(1) = rbmin_cld
165	rad_cld(1) = rmin_cld
166	vol_cld(1) = 4./3. * dble(pi) * rmin_cldrmin_cldrmin_cld
167	! vol_cld(1) = 4./3. * pi * rmin_cldrmin_cldrmin_cld
168
169	do i=1,nbin_cld-1
170	rad_cld(i+1) = rad_cld(i) * vrat_cld**(1./3.)
171	vol_cld(i+1) = vol_cld(i) * vrat_cld
172	enddo
173
174	do i=1,nbin_cld
175	rb_cld(i+1)= ( (2.vrat_cld) / (vrat_cld+1.) )(1./3.)
176	& rad_cld(i)
177	dr_cld(i) = rb_cld(i+1) - rb_cld(i)
178	enddo
179	rb_cld(nbin_cld+1) = rbmax_cld
180	dr_cld(nbin_cld) = rb_cld(nbin_cld+1) - rb_cld(nbin_cld)
181
182	print*, ' '
183	print*,'Microphysics: size bin information:'
184	print*,'i,rb_cld(i), rad_cld(i),dr_cld(i)'
185	print*,'-----------------------------------'
186	do i=1,nbin_cld
187	write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i),
188	& dr_cld(i)
189	enddo
190	write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1)
191	print*,'-----------------------------------'
192
193	do i=1,nbin_cld+1
194	! rb_cld(i) = log(rb_cld(i))
195	rb_cld(i) = dlog(rb_cld(i)) !! we save that so that it is not computed
196	!! at each timestep and gridpoint
197	enddo
198
199	c Contact parameter of water ice on dust ( m=cos(theta) )
200	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
201	! mteta = 0.95
202	write(,) 'water_param contact parameter:', mteta
203
204	c Volume of a water molecule (m3)
205	vo1 = mh2o / dble(rho_ice)
206	c Variance of the ice and CCN distributions
207	sigma_ice = sqrt(log(1.+nuice_sed))
208
209
210	write(,) 'Variance of ice & CCN distribs :', sigma_ice
211	write(,) 'nuice for sedimentation:', nuice_sed
212	write(,) 'Volume of a water molecule:', vo1
213
214
215	test_flag = .false.
216
217	firstcall=.false.
218	END IF
219
220
221	!=============================================================
222	! 1. Initialisation
223	!=============================================================
224	cste = 4pirho_ice*microtimestep
225
226	res_out(:,:) = 0
227	rice(:,:) = 1.e-8
228
229	c Initialize the tendencies
230	subpdqcloud(1:ngrid,1:nlay,1:nq)=0
231	subpdtcloud(1:ngrid,1:nlay)=0
232
233
234	zt(1:ngrid,1:nlay) =
235	& pteff(1:ngrid,1:nlay) +
236	& sum_subpdt(1:ngrid,1:nlay) * microtimestep
237
238	zq(1:ngrid,1:nlay,1:nq) =
239	& pqeff(1:ngrid,1:nlay,1:nq) +
240	& sum_subpdq(1:ngrid,1:nlay,1:nq) * microtimestep
241
242
243	WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
244	& zq(1:ngrid,1:nlay,1:nq) = 1.e-30
245
246	zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
247
248	!=============================================================
249	! 2. Compute saturation
250	!=============================================================
251
252
253	dev2 = 1. / ( sqrt(2.) * sigma_ice )
254
255	call watersat(ngrid*nlay,zt,pplay,zqsat)
256
257	countcells = 0
258
259	c Main loop over the GCM's grid
260	DO l=1,nlay
261	DO ig=1,ngrid
262
263	c Get the partial pressure of water vapor and its saturation ratio
264	ph2o = zq(ig,l,igcm_h2o_vap) * (mmean(ig,l)/18.) * pplay(ig,l)
265	satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
266
267	!=============================================================
268	! 3. Nucleation
269	!=============================================================
270
271	IF ( satu .ge. 1. ) THEN ! if there is condensation
272
273	call updaterccn(zq(ig,l,igcm_dust_mass),
274	& zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
275
276
277	c Expand the dust moments into a binned distribution
278	Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig) + 1.e-30
279	No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30
280	Rn = rdust(ig,l)
281	Rn = -dlog(Rn)
282	Rm = Rn - 3. * sigma_ice*sigma_ice
283	n_derf = derf( (rb_cld(1)+Rn) *dev2)
284	m_derf = derf( (rb_cld(1)+Rm) *dev2)
285	do i = 1, nbin_cld
286	n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
287	m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
288	n_derf = derf( (rb_cld(i+1)+Rn) *dev2)
289	m_derf = derf( (rb_cld(i+1)+Rm) *dev2)
290	n_aer(i) = n_aer(i) + 0.5 * No * n_derf
291	m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
292	enddo
293
294	! sumcheck = 0
295	! do i = 1, nbin_cld
296	! sumcheck = sumcheck + n_aer(i)
297	! enddo
298	! sumcheck = abs(sumcheck/No - 1)
299	! if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
300	! print*, "WARNING, No sumcheck PROBLEM"
301	! print*, "sumcheck, No",sumcheck, No
302	! print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
303	! print*, "Dust binned distribution", n_aer
304	! endif
305	!
306	! sumcheck = 0
307	! do i = 1, nbin_cld
308	! sumcheck = sumcheck + m_aer(i)
309	! enddo
310	! sumcheck = abs(sumcheck/Mo - 1)
311	! if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) then
312	! print*, "WARNING, Mo sumcheck PROBLEM"
313	! print*, "sumcheck, Mo",sumcheck, Mo
314	! print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
315	! print*, "Dust binned distribution", m_aer
316	! endif
317
318
319	c Get the rates of nucleation
320	call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
321
322	dN = 0.
323	dM = 0.
324	do i = 1, nbin_cld
325	n_aer(i) = n_aer(i)/( 1. + rate(i)*microtimestep)
326	m_aer(i) = m_aer(i)/( 1. + rate(i)*microtimestep)
327	dN = dN + n_aer(i) * rate(i) * microtimestep
328	dM = dM + m_aer(i) * rate(i) * microtimestep
329	enddo
330
331
332	c Update Dust particles
333	zq(ig,l,igcm_dust_mass) =
334	& zq(ig,l,igcm_dust_mass) - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
335	zq(ig,l,igcm_dust_number) =
336	& zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
337	c Update CCNs
338	zq(ig,l,igcm_ccn_mass) =
339	& zq(ig,l,igcm_ccn_mass) + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
340	zq(ig,l,igcm_ccn_number) =
341	& zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
342
343	ENDIF ! of is satu >1
344
345	!=============================================================
346	! 4. Ice growth: scheme for radius evolution
347	!=============================================================
348
349	c We trigger crystal growth if and only if there is at least one nuclei (N>1).
350	c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
351	c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
352
353	IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.) THEN ! we trigger crystal growth
354
355
356	call updaterice_micro(zq(ig,l,igcm_h2o_ice),
357	& zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
358	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
359
360	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
361
362	c saturation at equilibrium
363	c rice should not be too small, otherwise seq value is not valid
364	seq = exp(2.sig(zt(ig,l))mh2o / (rho_icergpzt(ig,l)*
365	& max(rice(ig,l),1.e-7)))
366
367	c get resistance growth
368	call growthrate(zt(ig,l),pplay(ig,l),
369	& real(ph2o/satu),rice(ig,l),res)
370
371	res_out(ig,l) = res
372
373	ccccccc implicit scheme of mass growth
374
375	dMice =
376	& (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))
377	& /(reszqsat(ig,l)/(csteNo*rice(ig,l)) + 1.)
378
379
380	! With the above scheme, dMice cannot be bigger than vapor,
381	! but can be bigger than all available ice.
382	dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
383	dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless...
384
385	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice
386	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice
387
388
389	countcells = countcells + 1
390
391	! latent heat release
392	lw=(2834.3-0.28*(zt(ig,l)-To)-
393	& 0.004(zt(ig,l)-To)(zt(ig,l)-To))*1.e+3
394	subpdtcloud(ig,l)= dMice*lw/cpp/microtimestep
395
396
397
398	!=============================================================
399	! 5. Dust cores released, tendancies, latent heat, etc ...
400	!=============================================================
401
402	c If all the ice particles sublimate, all the condensation
403	c nuclei are released:
404	if (zq(ig,l,igcm_h2o_ice).le.1.e-28) then
405
406	c Water
407	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)
408	& + zq(ig,l,igcm_h2o_ice)
409	zq(ig,l,igcm_h2o_ice) = 0.
410	c Dust particles
411	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass)
412	& + zq(ig,l,igcm_ccn_mass)
413	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number)
414	& + zq(ig,l,igcm_ccn_number)
415	c CCNs
416	zq(ig,l,igcm_ccn_mass) = 0.
417	zq(ig,l,igcm_ccn_number) = 0.
418
419	endif
420
421	ENDIF !of if Nccn>1
422
423	ENDDO ! of ig loop
424	ENDDO ! of nlayer loop
425
426
427	! Get cloud tendencies
428	subpdqcloud(1:ngrid,1:nlay,igcm_h2o_vap) =
429	& (zq(1:ngrid,1:nlay,igcm_h2o_vap) -
430	& zq0(1:ngrid,1:nlay,igcm_h2o_vap))/microtimestep
431	subpdqcloud(1:ngrid,1:nlay,igcm_h2o_ice) =
432	& (zq(1:ngrid,1:nlay,igcm_h2o_ice) -
433	& zq0(1:ngrid,1:nlay,igcm_h2o_ice))/microtimestep
434	subpdqcloud(1:ngrid,1:nlay,igcm_ccn_mass) =
435	& (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
436	& zq0(1:ngrid,1:nlay,igcm_ccn_mass))/microtimestep
437	subpdqcloud(1:ngrid,1:nlay,igcm_ccn_number) =
438	& (zq(1:ngrid,1:nlay,igcm_ccn_number) -
439	& zq0(1:ngrid,1:nlay,igcm_ccn_number))/microtimestep
440
441	if (scavenging) then
442
443	subpdqcloud(1:ngrid,1:nlay,igcm_dust_mass) =
444	& (zq(1:ngrid,1:nlay,igcm_dust_mass) -
445	& zq0(1:ngrid,1:nlay,igcm_dust_mass))/microtimestep
446	subpdqcloud(1:ngrid,1:nlay,igcm_dust_number) =
447	& (zq(1:ngrid,1:nlay,igcm_dust_number) -
448	& zq0(1:ngrid,1:nlay,igcm_dust_number))/microtimestep
449
450	endif
451
452
453
454	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
455	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
456	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
457	IF (test_flag) then
458
459	! error2d(:) = 0.
460	DO l=1,nlay
461	DO ig=1,ngrid
462	! error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
463	satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l)
464	satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
465	ENDDO
466	ENDDO
467
468	print*, 'count is ',countcells, ' i.e. ',
469	& countcells100/(nlayngrid), '% for microphys computation'
470
471	#ifndef MESOSCALE
472	! IF (ngrid.ne.1) THEN ! 3D
473	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
474	! & satu_out)
475	! call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
476	! & dM_out)
477	! call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
478	! & dN_out)
479	! call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2,
480	! & error2d)
481	! call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3,
482	! & zqsat)
483	! ENDIF
484
485	! IF (ngrid.eq.1) THEN ! 1D
486	! call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1,
487	! & error_out)
488	call WRITEdiagfi(ngrid,"resist","resistance","s/m2",1,
489	& res_out)
490	call WRITEdiagfi(ngrid,"satu_bf","satu before","kg/kg",1,
491	& satubf)
492	call WRITEdiagfi(ngrid,"satu_af","satu after","kg/kg",1,
493	& satuaf)
494	call WRITEdiagfi(ngrid,"vapbf","h2ovap before","kg/kg",1,
495	& zq0(1,1,igcm_h2o_vap))
496	call WRITEdiagfi(ngrid,"vapaf","h2ovap after","kg/kg",1,
497	& zq(1,1,igcm_h2o_vap))
498	call WRITEdiagfi(ngrid,"icebf","h2oice before","kg/kg",1,
499	& zq0(1,1,igcm_h2o_ice))
500	call WRITEdiagfi(ngrid,"iceaf","h2oice after","kg/kg",1,
501	& zq(1,1,igcm_h2o_ice))
502	call WRITEdiagfi(ngrid,"ccnbf","ccn before","/kg",1,
503	& zq0(1,1,igcm_ccn_number))
504	call WRITEdiagfi(ngrid,"ccnaf","ccn after","/kg",1,
505	& zq(1,1,igcm_ccn_number))
506	c call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
507	c & gr_out)
508	c call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
509	c & rate_out)
510	c call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
511	c & dM_out)
512	c call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
513	c & dN_out)
514	call WRITEdiagfi(ngrid,"zqsat","p vap sat","kg/kg",1,
515	& zqsat)
516	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
517	! & satu_out)
518	call WRITEdiagfi(ngrid,"rice","ice radius","m",1,
519	& rice)
520	! call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
521	! & rdust)
522	! call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
523	! & rsedcloud)
524	! call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
525	! & rhocloud)
526	! ENDIF
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
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
580	END SUBROUTINE improvedclouds
581
582	END MODULE improvedclouds_mod

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

Download in other formats:

Original Format