Context Navigation

co2cloud.F @ 1652

Last change on this file since 1652 was 1651, checked in by jaudouard, 8 years ago
Modification of co2 clouds microtimestep and addition of a meteoritic flux of dust
Property svn:executable set to ``*
File size: 28.5 KB

Line
1	SUBROUTINE co2cloud(ngrid,nlay,ptimestep,
2	& pplev,pplay,pdpsrf,pzlay,pt,pdt,
3	& pq,pdq,pdqcloudco2,pdtcloudco2,
4	& nq,tau,tauscaling,rdust,rice,riceco2,nuice,
5	& rsedcloudco2,rhocloudco2,zlev,pdqs_sedco2)
6	! to use 'getin'
7	use dimradmars_mod, only: naerkind
8	USE comcstfi_h
9	USE ioipsl_getincom
10	USE updaterad
11	use conc_mod, only: mmean
12	use tracer_mod, only: nqmx, igcm_co2, igcm_co2_ice,
13	& igcm_dust_mass, igcm_dust_number,
14	& igcm_ccnco2_mass, igcm_ccnco2_number,
15	& rho_dust, nuiceco2_sed, nuiceco2_ref,
16	& rho_ice_co2,r3n_q
17	USE meteoriticflux_mod, ONLY:meteo_flux_mass,meteo_flux_number
18	& ,meteo_alt
19
20	IMPLICIT NONE
21
22
23	c=======================================================================
24	c CO2 clouds formation
25	c
26	c There is a time loop specific to cloud formation
27	c due to timescales smaller than the GCM integration timestep.
28	c microphysics subroutine is improvedCO2clouds.F
29	c
30	c The co2 clouds tracers (co2_ice, ccn mass and concentration) are
31	c sedimented at each microtimestep. pdqs_sedco2 keeps track of the
32	c CO2 flux at the surface
33	c
34	c Authors: 09/2016 Joachim Audouard & Constantino Listowski
35	c Adaptation of the water ice clouds scheme (with specific microphysics)
36	c of Montmessin, Navarro & al.
37	c
38	c
39	c=======================================================================
40
41	c-----------------------------------------------------------------------
42	c declarations:
43	c -------------
44
45	!#include "dimensions.h"
46	!#include "dimphys.h"
47	#include "callkeys.h"
48	!#include "tracer.h"
49	!#include "comgeomfi.h"
50	!#include "dimradmars.h"
51	! naerkind is set in scatterers.h (built when compiling with makegcm -s #)
52	!#include"scatterers.h"
53	#include "microphys.h"
54
55
56	c Inputs:
57	c ------
58
59	INTEGER ngrid,nlay
60	INTEGER nq ! nombre de traceurs
61	REAL ptimestep ! pas de temps physique (s)
62	REAL pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa)
63	REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
64	REAL pdpsrf(ngrid) ! tendence surf pressure
65	REAL pzlay(ngrid,nlay) ! altitude at the middle of the layers
66	REAL pt(ngrid,nlay) ! temperature at the middle of the layers (K)
67	REAL pdt(ngrid,nlay) ! tendence temperature des autres param.
68	real,intent(in) :: zlev(ngrid,nlay+1) ! altitude at the boundaries of the layers
69
70	real pq(ngrid,nlay,nq) ! traceur (kg/kg)
71	real pdq(ngrid,nlay,nq) ! tendance avant condensation (kg/kg.s-1)
72
73	real rice(ngrid,nlay) ! Water Ice mass mean radius (m)
74	! used for nucleation of CO2 on ice-coated ccns
75
76	REAL tau(ngrid,naerkind) ! Column dust optical depth at each point
77	REAL tauscaling(ngrid) ! Convertion factor for dust amount
78	real rdust(ngrid,nlay) ! Dust geometric mean radius (m)
79
80	c Outputs:
81	c -------
82
83	real pdqcloudco2(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1)
84	REAL pdtcloudco2(ngrid,nlay) ! tendence temperature due
85	! a la chaleur latente
86
87	DOUBLE PRECISION riceco2(ngrid,nlay) ! Ice mass mean radius (m)
88	! (r_c in montmessin_2004)
89	REAL nuice(ngrid,nlay) ! Estimated effective variance
90	! of the size distribution
91	real rsedcloudco2(ngrid,nlay) ! Cloud sedimentation radius
92	real rhocloudco2(ngrid,nlay) ! Cloud density (kg.m-3)
93	real rhocloudco2t(ngrid,nlay) ! Cloud density (kg.m-3)
94	real pdqs_sedco2(ngrid) ! CO2 flux at the surface
95	c local:
96	c ------
97
98	! for ice radius computation
99	REAL Mo,No
100	REAl ccntyp
101
102	! for time loop
103	INTEGER microstep ! time subsampling step variable
104	INTEGER imicro ! time subsampling for coupled water microphysics & sedimentation
105	SAVE imicro
106	REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation
107	SAVE microtimestep
108
109	! tendency given by clouds (inside the micro loop)
110	REAL subpdqcloudco2(ngrid,nlay,nq) ! cf. pdqcloud
111	REAL subpdtcloudco2(ngrid,nlay) ! cf. pdtcloud
112
113	! global tendency (clouds+physics)
114	REAL subpdq(ngrid,nlay,nq) ! cf. pdqcloud
115	REAL subpdt(ngrid,nlay) ! cf. pdtcloud
116	real wq(ngrid,nlay+1) ! ! displaced tracer mass (kg.m-2) during microtimestep because sedim (?/m-2)
117
118	REAL satuco2(ngrid,nlay) ! co2 satu ratio for output
119	REAL zqsatco2(ngrid,nlay) ! saturation co2
120
121	INTEGER iq,ig,l
122	LOGICAL,SAVE :: firstcall=.true.
123	DOUBLE PRECISION Nccnco2, Niceco2,mdustJA,ndustJA
124	DOUBLE PRECISION Qccnco2
125	real :: beta
126
127	real epaisseur (ngrid,nlay) ! Layer thickness (m)
128	real masse (ngrid,nlay) ! Layer mass (kg.m-2)
129
130	double precision diff,diff0
131	integer meteo_lvl
132	real tempo_traceur_t(ngrid,nlay)
133	real tempo_traceurs(ngrid,nlay,nq)
134	real sav_trac(ngrid,nlay,nq)
135	real pdqsed(ngrid,nlay,nq)
136	c ** un petit test de coherence
137	c --------------------------
138
139	IF (firstcall) THEN
140
141	if (nq.gt.nqmx) then
142	write(,) 'stop in co2cloud (nq.gt.nqmx)!'
143	write(,) 'nq=',nq,' nqmx=',nqmx
144	stop
145	endif
146
147	write(,) "co2cloud: igcm_co2=",igcm_co2
148	write(,) " igcm_co2_ice=",igcm_co2_ice
149
150	write(,) "time subsampling for microphysic ?"
151	#ifdef MESOSCALE
152	imicro = 2
153	#else
154	imicro = 30
155	#endif
156	call getin("imicro",imicro)
157	write(,)"imicro = ",imicro
158
159	microtimestep = ptimestep/real(imicro)
160	write(,)"Physical timestep is",ptimestep
161	write(,)"CO2 Microphysics timestep is",microtimestep
162
163
164	write(,) "Warning ! Meteoritic flux of dust is turned on"
165	write(,) "Dust mass = ",meteo_flux_mass
166	write(,) "Dust number = ",meteo_flux_number
167	write(,) "Are added at the z-level (km)",meteo_alt
168	write(,) "Every timestep in co2cloud.F"
169
170	firstcall=.false.
171	ENDIF ! of IF (firstcall)
172
173	c-----Initialization
174
175	c add meteoritic flux of dust
176	!convert meteo_alt (in km) to z-level
177	!pzlay altitudes of the layers
178
179	do ig=1, ngrid
180	diff0=1000
181	meteo_lvl=1
182	do l=1,nlay
183	diff=pzlay(ig,l)/1000-meteo_alt
184	if (abs(diff) .lt. diff0) then
185	diff0=abs(diff)
186	meteo_lvl=l
187	endif
188	enddo
189	c write(,) "meteo_lvl=",meteo_lvl
190	pq(ig,meteo_lvl,igcm_dust_mass)=pq(ig,meteo_lvl,igcm_dust_mass)
191	& +meteo_flux_mass
192	pq(ig,meteo_lvl,igcm_dust_number)=
193	& pq(ig,meteo_lvl,igcm_dust_number)+meteo_flux_number
194	enddo
195
196	call WRITEDIAGFI(ngrid,"pzlay","altitude","km",3,
197	& pzlay)
198	beta=0.85
199	subpdq(1:ngrid,1:nlay,1:nq) = 0
200	subpdt(1:ngrid,1:nlay) = 0
201	subpdqcloudco2(1:ngrid,1:nlay,1:nq) = 0
202	subpdtcloudco2(1:ngrid,1:nlay) = 0
203
204
205	wq(:,:)=0
206	! default value if no ice
207	rhocloudco2(1:ngrid,1:nlay) = rho_dust
208	rhocloudco2t(1:ngrid,1:nlay) = rho_dust
209	epaisseur(1:ngrid,1:nlay)=0
210	masse(1:ngrid,1:nlay)=0
211
212	tempo_traceur_t(1:ngrid,1:nlay)=0
213	tempo_traceurs(1:ngrid,1:nlay,1:nq)=0
214	sav_trac(1:ngrid,1:nlay,1:nq)=0
215	pdqsed(1:ngrid,1:nlay,1:nq)=0
216
217	do l=1,nlay
218	do ig=1, ngrid
219	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
220	epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
221
222	enddo
223	enddo
224
225
226
227
228
229
230
231
232	c-------------------------------------------------------------------
233	c 1. Tendencies:
234	c------------------
235
236
237
238	c------------------------------------------------------------------
239	c Time subsampling for microphysics
240	c------------------------------------------------------------------
241	DO microstep=1,imicro
242	c------ Temperature tendency subpdt
243	! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt
244	! If imicro=1 subpdt is the same as pdt
245	DO l=1,nlay
246	DO ig=1,ngrid
247	c tempo_traceur_t(ig,l)=tempo_traceur_t(ig,l)
248	c & + subpdtcloudco2(ig,l)
249	!write(,) 'T micro= ', tempo_traceur_t(ig,l)
250	c tempo_traceurs(ig,l,:)=tempo_traceurs(ig,l,:)
251	c & +subpdqcloudco2(ig,l,:)
252
253	subpdt(ig,l) = subpdt(ig,l)
254	& + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
255
256	subpdq(ig,l,igcm_dust_mass) =
257	& subpdq(ig,l,igcm_dust_mass)
258	& + pdq(ig,l,igcm_dust_mass)
259
260	subpdq(ig,l,igcm_dust_number) =
261	& subpdq(ig,l,igcm_dust_number)
262	& + pdq(ig,l,igcm_dust_number)
263
264	subpdq(ig,l,igcm_ccnco2_mass) =
265	& subpdq(ig,l,igcm_ccnco2_mass)
266	& + pdq(ig,l,igcm_ccnco2_mass)
267
268	subpdq(ig,l,igcm_ccnco2_number) =
269	& subpdq(ig,l,igcm_ccnco2_number)
270	& + pdq(ig,l,igcm_ccnco2_number)
271
272	subpdq(ig,l,igcm_co2_ice) =
273	& subpdq(ig,l,igcm_co2_ice)
274	& + pdq(ig,l,igcm_co2_ice)
275	subpdq(ig,l,igcm_co2) =
276	& subpdq(ig,l,igcm_co2)
277	& + pdq(ig,l,igcm_co2)
278
279	tempo_traceur_t(ig,l)= pt(ig,l)+subpdt(ig,l)*microtimestep
280	tempo_traceurs(ig,l,:)= pq(ig,l,:)+subpdq(ig,l,:)
281	& *microtimestep
282	!Stepped entry for sedimentation
283	ENDDO
284	ENDDO
285
286	!RSEDCLOUD AND RICECO2 HERE
287
288	DO l=1, nlay
289	DO ig=1,ngrid
290	Niceco2=max(tempo_traceurs(ig,l,igcm_co2_ice),1.e-30)
291	Nccnco2=max(tempo_traceurs(ig,l,igcm_ccnco2_number),
292	& 1.e-30)
293	Qccnco2=max(tempo_traceurs(ig,l,igcm_ccnco2_mass),
294	& 1.e-30)
295	mdustJA= tempo_traceurs(ig,l,igcm_dust_mass)
296	ndustJA=tempo_traceurs(ig,l,igcm_dust_number)
297	if ((ndustJA .lt. tauscaling(ig)) .or. (mdustJA .lt.
298	& 1.e-30 *tauscaling(ig))) then
299	rdust(ig,l)=1.e-10
300	else
301	rdust(ig,l)=(3./4./pi/2500.mdustJA/ndustJA)*(1./3.)
302	rdust(ig,l)=min(rdust(ig,l),5.e-6)
303	rdust(ig,l)=max(rdust(ig,l),1.e-9)
304	end if
305	c rhocloudco2t(ig,l) = (Niceco2 *rho_ice_co2
306	c & + Qccnco2*rho_dust)
307	c & / (Niceco2 + Qccnco2)
308	riceco2(ig,l)= Niceco2*3.0/
309	& (4.0rho_ice_co2pi*Nccnco2)
310	& +rdust(ig,l)rdust(ig,l)rdust(ig,l)
311	riceco2(ig,l)=riceco2(ig,l)**(1.0/3.0)
312	c write(,) "in co2clouds, rice = ",riceco2(ig,l)
313	c write(,) "in co2clouds, rho = ",rhocloudco2t(ig,l)
314
315	call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2,
316	& tauscaling(ig),riceco2(ig,l),rhocloudco2t(ig,l))
317	c write(,) "in co2clouds, rice update = ",riceco2(ig,l)
318	c write(,) "in co2clouds, rho update = "
319	c & ,rhocloudco2t(ig,l)
320
321	rsedcloudco2(ig,l)=max(riceco2(ig,l)*
322	& (1.+nuiceco2_sed)(1.+nuiceco2_sed)(1.+nuiceco2_sed),
323	& rdust(ig,l))
324	rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),5.e-4)
325	c write(,) 'Rsedcloud = ',rsedcloudco2(ig,l)
326	!write(,) 'Rhocloudco2 = ',rhocloudco2t(ig,l)
327
328	ENDDO
329	ENDDO
330
331	! Gravitational sedimentation
332
333	! sedimentation computed from radius computed from q in module radii_mod
334	sav_trac(:,:,igcm_co2_ice)=tempo_traceurs(:,:,igcm_co2_ice)
335	sav_trac(:,:,igcm_ccnco2_mass)=
336	& tempo_traceurs(:,:,igcm_ccnco2_mass)
337	sav_trac(:,:,igcm_ccnco2_number)=
338	& tempo_traceurs(:,:,igcm_ccnco2_number)
339
340	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
341	& microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
342	& rsedcloudco2,rhocloudco2t,
343	& tempo_traceurs(:,:,igcm_co2_ice),wq,beta) ! 3 traceurs
344
345	! sedim at the surface of co2 ice
346	do ig=1,ngrid
347	pdqs_sedco2(ig)=pdqs_sedco2(ig)+ wq(ig,1)
348	end do
349
350	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
351	& microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
352	& rsedcloudco2,rhocloudco2t,
353	& tempo_traceurs(:,:,igcm_ccnco2_mass),wq,beta)
354
355	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
356	& microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
357	& rsedcloudco2,rhocloudco2t,
358	& tempo_traceurs(:,:,igcm_ccnco2_number),wq,beta)
359
360
361	DO l = 1, nlay
362	DO ig=1,ngrid
363	pdqsed(ig,l,igcm_ccnco2_mass)=
364	& (tempo_traceurs(ig,l,igcm_ccnco2_mass)-
365	& sav_trac(ig,l,igcm_ccnco2_mass))/microtimestep
366	pdqsed(ig,l,igcm_ccnco2_number)=
367	& (tempo_traceurs(ig,l,igcm_ccnco2_number)-
368	& sav_trac(ig,l,igcm_ccnco2_number))/microtimestep
369	pdqsed(ig,l,igcm_co2_ice)=
370	& (tempo_traceurs(ig,l,igcm_co2_ice)-
371	& sav_trac(ig,l,igcm_co2_ice))/microtimestep
372	ENDDO
373	ENDDO
374	!pdqsed est la tendance due a la sedimentation
375
376	DO l = 1, nlay
377	DO ig=1,ngrid
378	pdqsed(ig,l,igcm_ccnco2_mass)=
379	& (tempo_traceurs(ig,l,igcm_ccnco2_mass)-
380	& sav_trac(ig,l,igcm_ccnco2_mass))/microtimestep
381	pdqsed(ig,l,igcm_ccnco2_number)=
382	& (tempo_traceurs(ig,l,igcm_ccnco2_number)-
383	& sav_trac(ig,l,igcm_ccnco2_number))/microtimestep
384	pdqsed(ig,l,igcm_co2_ice)=
385	& (tempo_traceurs(ig,l,igcm_co2_ice)-
386	& sav_trac(ig,l,igcm_co2_ice))/microtimestep
387	ENDDO
388	ENDDO
389	!pdqsed est la tendance due a la sedimentation
390	DO l=1,nlay
391	DO ig=1,ngrid
392	subpdq(ig,l,igcm_ccnco2_mass) =
393	& subpdq(ig,l,igcm_ccnco2_mass)
394	& +pdqsed(ig,l,igcm_ccnco2_mass)
395
396	subpdq(ig,l,igcm_ccnco2_number) =
397	& subpdq(ig,l,igcm_ccnco2_number)
398	& +pdqsed(ig,l,igcm_ccnco2_number)
399
400	subpdq(ig,l,igcm_co2_ice) =
401	& subpdq(ig,l,igcm_co2_ice)
402	& +pdqsed(ig,l,igcm_co2_ice)
403	ENDDO
404	ENDDO
405	c-------------------------------------------------------------------
406	c 2. Main call to the different cloud schemes:
407	c------------------------------------------------
408	IF (microphysco2) THEN
409	CALL improvedCO2clouds(ngrid,nlay,microtimestep,
410	& pplay,pt,subpdt,
411	& pq,subpdq,subpdqcloudco2,subpdtcloudco2,
412	& nq,tauscaling)
413
414	ELSE
415
416	write(,) ' no simpleCO2clouds procedure: STOP' ! listo
417	STOP
418
419	c CALL simpleclouds(ngrid,nlay,microtimestep, ! for water-ice clouds
420	c & pplay,pzlay,pt,subpdt,
421	c & pq,subpdq,subpdqcloud,subpdtcloud,
422	c & nq,tau,riceco2)
423	ENDIF
424
425
426	c-------------------------------------------------------------------
427	c 3. Updating tendencies after cloud scheme:
428	c-----------------------------------------------
429
430	c IF (microphysco2) THEN
431	DO l=1,nlay
432	DO ig=1,ngrid
433	subpdq(ig,l,igcm_dust_mass) =
434	& subpdq(ig,l,igcm_dust_mass)
435	& + subpdqcloudco2(ig,l,igcm_dust_mass)
436
437	subpdq(ig,l,igcm_dust_number) =
438	& subpdq(ig,l,igcm_dust_number)
439	& + subpdqcloudco2(ig,l,igcm_dust_number)
440
441	subpdq(ig,l,igcm_ccnco2_mass) =
442	& subpdq(ig,l,igcm_ccnco2_mass)
443	& + subpdqcloudco2(ig,l,igcm_ccnco2_mass)
444	c & +pdqsed(ig,l,igcm_ccnco2_mass)
445
446	subpdq(ig,l,igcm_ccnco2_number) =
447	& subpdq(ig,l,igcm_ccnco2_number)
448	& + subpdqcloudco2(ig,l,igcm_ccnco2_number)
449	c & +pdqsed(ig,l,igcm_ccnco2_number)
450
451	subpdq(ig,l,igcm_co2_ice) =
452	& subpdq(ig,l,igcm_co2_ice)
453	& + subpdqcloudco2(ig,l,igcm_co2_ice)
454	c & +pdqsed(ig,l,igcm_co2_ice)
455
456	subpdq(ig,l,igcm_co2) =
457	& subpdq(ig,l,igcm_co2)
458	& + subpdqcloudco2(ig,l,igcm_co2)
459	ENDDO
460	ENDDO
461
462
463	!ici
464	! call WRITEdiagfi(ngrid,"co2cloud000","co2 traceur","kg/kg",1,
465	! & pq(1,:,igcm_co2_ice) + ptimestep*
466	! & ( subpdq(1,:,igcm_co2_ice)))
467
468
469	IF (activice) THEN
470	DO l=1,nlay
471	DO ig=1,ngrid
472	subpdt(ig,l) =
473	& subpdt(ig,l) + subpdtcloudco2(ig,l)
474	ENDDO
475	ENDDO
476	ENDIF
477
478
479	ENDDO ! of DO microstep=1,imicro
480
481	c-------------------------------------------------------------------
482	c 6. Compute final tendencies after time loop:
483	c------------------------------------------------
484	c CO2 flux at surface (kg.m-2.s-1)
485	do ig=1,ngrid
486	pdqs_sedco2(ig)=pdqs_sedco2(ig)/ptimestep
487	enddo
488
489	c------ Temperature tendency pdtcloud
490	DO l=1,nlay
491	DO ig=1,ngrid
492	pdtcloudco2(ig,l) =
493	& subpdt(ig,l)/imicro-pdt(ig,l)
494	ENDDO
495	ENDDO
496
497	c------ Tracers tendencies pdqcloud
498	DO l=1,nlay
499	DO ig=1,ngrid
500
501	pdqcloudco2(ig,l,igcm_co2_ice) =
502	& subpdq(ig,l,igcm_co2_ice)/imicro
503	& - pdq(ig,l,igcm_co2_ice)
504	pdqcloudco2(ig,l,igcm_co2) =
505	& subpdq(ig,l,igcm_co2)/imicro
506	& - pdq(ig,l,igcm_co2)
507	ENDDO
508	ENDDO
509
510
511	! call WRITEdiagfi(ngrid,"co2cloud00","co2 traceur","kg/kg",1,
512	! & pq(1,:,igcm_co2_ice) + ptimestep*
513	! & (pdq(1,:,igcm_co2_ice) + pdqcloudco2(1,:,igcm_co2_ice)))
514
515
516	IF(microphysco2) THEN
517	DO l=1,nlay
518	DO ig=1,ngrid
519	pdqcloudco2(ig,l,igcm_ccnco2_mass) =
520	& subpdq(ig,l,igcm_ccnco2_mass)/imicro
521	& - pdq(ig,l,igcm_ccnco2_mass)
522	pdqcloudco2(ig,l,igcm_ccnco2_number) =
523	& subpdq(ig,l,igcm_ccnco2_number)/imicro
524	& - pdq(ig,l,igcm_ccnco2_number)
525	ENDDO
526	ENDDO
527	ENDIF
528
529
530	IF(scavenging) THEN
531	DO l=1,nlay
532	DO ig=1,ngrid
533	pdqcloudco2(ig,l,igcm_dust_mass) =
534	& subpdq(ig,l,igcm_dust_mass)/real(imicro)
535	& - pdq(ig,l,igcm_dust_mass)
536	pdqcloudco2(ig,l,igcm_dust_number) =
537	& subpdq(ig,l,igcm_dust_number)/real(imicro)
538	& - pdq(ig,l,igcm_dust_number)
539	ENDDO
540	ENDDO
541	ENDIF
542
543	c ENDIF
544	c------- Due to stepped entry, other processes tendencies can add up to negative values
545	c------- Therefore, enforce positive values and conserve mass
546
547
548	IF(microphysco2) THEN
549	DO l=1,nlay
550	DO ig=1,ngrid
551	IF ((pq(ig,l,igcm_ccnco2_number) +
552	& ptimestep* (pdq(ig,l,igcm_ccnco2_number) +
553	& pdqcloudco2(ig,l,igcm_ccnco2_number))
554	& .lt. 0)
555	& .or. (pq(ig,l,igcm_ccnco2_mass) +
556	& ptimestep* (pdq(ig,l,igcm_ccnco2_mass) +
557	& pdqcloudco2(ig,l,igcm_ccnco2_mass))
558	& .lt. 0)) THEN
559
560	pdqcloudco2(ig,l,igcm_ccnco2_number) =
561	& - pq(ig,l,igcm_ccnco2_number)/ptimestep
562	& - pdq(ig,l,igcm_ccnco2_number) +0
563
564	pdqcloudco2(ig,l,igcm_dust_number) =
565	& -pdqcloudco2(ig,l,igcm_ccnco2_number)
566
567	pdqcloudco2(ig,l,igcm_ccnco2_mass) =
568	& - pq(ig,l,igcm_ccnco2_mass)/ptimestep
569	& - pdq(ig,l,igcm_ccnco2_mass)+0
570
571	pdqcloudco2(ig,l,igcm_dust_mass) =
572	& -pdqcloudco2(ig,l,igcm_ccnco2_mass)
573
574	ENDIF
575	ENDDO
576	ENDDO
577	ENDIF
578
579
580	IF(scavenging) THEN
581	DO l=1,nlay
582	DO ig=1,ngrid
583	IF ( (pq(ig,l,igcm_dust_number) +
584	& ptimestep* (pdq(ig,l,igcm_dust_number) +
585	& pdqcloudco2(ig,l,igcm_dust_number)) .lt. 0.)
586	& .or. (pq(ig,l,igcm_dust_mass)+
587	& ptimestep* (pdq(ig,l,igcm_dust_mass) +
588	& pdqcloudco2(ig,l,igcm_dust_mass))
589	& .lt. 0.)) then
590
591	pdqcloudco2(ig,l,igcm_dust_number) =
592	& - pq(ig,l,igcm_dust_number)/ptimestep
593	& - pdq(ig,l,igcm_dust_number)+0
594
595	pdqcloudco2(ig,l,igcm_ccnco2_number) =
596	& -pdqcloudco2(ig,l,igcm_dust_number)
597
598	pdqcloudco2(ig,l,igcm_dust_mass) =
599	& - pq(ig,l,igcm_dust_mass)/ptimestep
600	& - pdq(ig,l,igcm_dust_mass) +0
601
602	pdqcloudco2(ig,l,igcm_ccnco2_mass) =
603	& -pdqcloudco2(ig,l,igcm_dust_mass)
604	ENDIF
605	ENDDO
606	ENDDO
607	ENDIF !pq+ptime*(pdq+pdqc)=1 ! pdqc=1-pq/ptime-pdq
608
609
610	DO l=1,nlay
611	DO ig=1,ngrid
612	IF (pq(ig,l,igcm_co2_ice) + ptimestep*
613	& (pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice))
614	& .lt. 1.e-25) THEN
615	pdqcloudco2(ig,l,igcm_co2_ice) =
616	& - pq(ig,l,igcm_co2_ice)/ptimestep - pdq(ig,l,igcm_co2_ice)
617	& +1.e-25
618	pdqcloudco2(ig,l,igcm_co2) = -pdqcloudco2(ig,l,igcm_co2_ice)
619	ENDIF
620	ENDDO
621	ENDDO
622
623
624
625
626	c------Update the ice and dust particle size "riceco2" for output or photochemistry
627	c------Only rsedcloudco2 is used for the co2 (cloud) cycle
628
629	IF(scavenging) THEN
630	DO l=1, nlay
631	DO ig=1,ngrid
632
633	c call updaterdust(
634	c & pq(ig,l,igcm_dust_mass) + ! dust mass
635	c & (pdq(ig,l,igcm_dust_mass) + ! dust mass
636	c & pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep, ! dust mass
637	c & pq(ig,l,igcm_dust_number) + ! dust number
638	c & (pdq(ig,l,igcm_dust_number) + ! dust number
639	c & pdqcloudco2(ig,l,igcm_dust_number))*ptimestep, ! dust number
640	c & rdust(ig,l))
641	c write(,) "in co2clouds, rdust(ig,l)= ",rdust(ig,l)
642	mdustJA= pq(ig,l,igcm_dust_mass) +
643	& (pdq(ig,l,igcm_dust_mass) +
644	& pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep
645	ndustJA=pq(ig,l,igcm_dust_number) +
646	& (pdq(ig,l,igcm_dust_number) +
647	& pdqcloudco2(ig,l,igcm_dust_number))*ptimestep
648	if ((ndustJA .lt. tauscaling(ig)) .or. (mdustJA .lt.
649	& 1.e-30 *tauscaling(ig))) then
650	rdust(ig,l)=1.e-10
651	else
652	rdust(ig,l)=(3./4./pi/2500.mdustJA/ndustJA)*(1./3.)
653	rdust(ig,l)=min(rdust(ig,l),5.e-4)
654	rdust(ig,l)=max(rdust(ig,l),1.e-10)
655	endif
656	ENDDO
657	ENDDO
658	ENDIF
659
660
661	IF(microphysco2) THEN
662
663	DO l=1, nlay
664	DO ig=1,ngrid
665
666	c call updaterice_microco2(
667	c & pq(ig,l,igcm_co2_ice) + ! ice mass
668	c & (pdq(ig,l,igcm_co2_ice) + ! ice mass
669	c & pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep, ! ice mass
670	c & pq(ig,l,igcm_ccnco2_mass) + ! ccn mass
671	c & (pdq(ig,l,igcm_ccnco2_mass) + ! ccn mass
672	c & pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep, ! ccn mass
673	c & pq(ig,l,igcm_ccnco2_number) + ! ccn number
674	c & (pdq(ig,l,igcm_ccnco2_number) + ! ccn number
675	c & pdqcloudco2(ig,l,igcm_ccnco2_number))*ptimestep, ! ccn number
676	c & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
677	c write(,) "in co2clouds, riceco2(ig,l)= ",riceco2(ig,l)
678
679
680	Niceco2=pq(ig,l,igcm_co2_ice) +
681	& (pdq(ig,l,igcm_co2_ice) +
682	& pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep
683	Nccnco2=max((pq(ig,l,igcm_ccnco2_number) +
684	& (pdq(ig,l,igcm_ccnco2_number) +
685	& pdqcloudco2(ig,l,igcm_ccnco2_number))ptimestep)
686	& tauscaling(ig),1.e-30)
687	Qccnco2=max((pq(ig,l,igcm_ccnco2_mass) +
688	& (pdq(ig,l,igcm_ccnco2_mass) +
689	& pdqcloudco2(ig,l,igcm_ccnco2_mass))ptimestep)
690	& tauscaling(ig),1.e-30)
691	c rhocloudco2t(ig,l) = (Niceco2 rho_ice_co2 + Qccnco2rho_dust)
692	c & / (Niceco2 + Qccnco2)
693	c rhocloudco2(ig,l) = min(max(rhocloudco2t,rho_ice_co2),rho_dust)
694
695	c write(,) "test, nccnco2 =",nccnco22
696
697
698	riceco2(ig,l)= Niceco2*3.0/
699	& (4.0rho_ice_co2pi*Nccnco2)
700	& +rdust(ig,l)rdust(ig,l)rdust(ig,l)
701
702	riceco2(ig,l)=riceco2(ig,l)**(1.0/3.0)
703	c write(,) "In co2cloud, after loop, riceco2 =",riceco2(ig,l)
704	c write(,) "In co2cloud, after loop, rhoco2 ="
705	c & ,rhocloudco2t(ig,l)
706
707	call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2,
708	& tauscaling(ig),riceco2(ig,l),rhocloudco2t(ig,l))
709
710	c write(,) "In co2cloud, after loop and update, riceco2 ="
711	c & ,riceco2(ig,l)
712	c write(,) "In co2cloud, after loop and update, rhoco2 ="
713	c & ,rhocloudco2t(ig,l)
714
715	if ( Niceco2
716	& .le. 1.e-23 .or. riceco2(ig,l) .le. 1.e-10 .or.
717	& riceco2(ig,l) .ge. 4.999e-4) then ! .or. riceco2(ig,l) .gt. 1.e-4 ) then
718	riceco2(ig,l)=0.
719
720	!NO CLOUD : RESET TRACER AND CONSERVE MASS
721	pdqcloudco2(ig,l,igcm_co2)= pq(ig,l,igcm_co2_ice)
722	& /ptimestep+pdq(ig,l,igcm_co2_ice)
723
724	pdqcloudco2(ig,l,igcm_co2_ice)=-pq(ig,l,igcm_co2_ice)
725	& /ptimestep-pdq(ig,l,igcm_co2_ice)
726
727	pdqcloudco2(ig,l,igcm_ccnco2_mass)=
728	& -pq(ig,l,igcm_ccnco2_mass)
729	& /ptimestep-pdq(ig,l,igcm_ccnco2_mass)
730
731	pdqcloudco2(ig,l,igcm_ccnco2_number)=
732	& -pq(ig,l,igcm_ccnco2_number)
733	& /ptimestep-pdq(ig,l,igcm_ccnco2_number)
734
735	pdqcloudco2(ig,l,igcm_dust_number)=
736	& pq(ig,l,igcm_ccnco2_number)
737	& /ptimestep+pdq(ig,l,igcm_ccnco2_number)
738
739	pdqcloudco2(ig,l,igcm_dust_mass)=
740	& pq(ig,l,igcm_ccnco2_mass)
741	& /ptimestep+pdq(ig,l,igcm_ccnco2_mass)
742	c$$$
743
744	c$$$
745	endif
746
747	c write(,) "in co2clouds, riceco2(ig,l)v2= ",riceco2(ig,l)
748
749	ENDDO
750	ENDDO
751
752	ELSE ! no microphys ! not of concern for co2 clouds - listo
753
754	ENDIF ! of IF(microphysco2)
755
756
757	c TO CHEK for relevancy - listo
758
759	c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
760	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
761	c Then that should not affect the ice particle radius
762	do ig=1,ngrid
763	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,2)))then
764	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,3)))
765	& riceco2(ig,2)=riceco2(ig,3)
766	riceco2(ig,1)=riceco2(ig,2)
767	end if
768	end do
769
770
771	DO l=1,nlay
772	DO ig=1,ngrid
773	rsedcloudco2(ig,l)=max(riceco2(ig,l)*
774	& (1.+nuiceco2_sed)(1.+nuiceco2_sed)(1.+nuiceco2_sed),
775	& rdust(ig,l))
776	rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),1.e-4)
777	ENDDO
778	ENDDO
779
780	call co2sat(ngrid*nlay,pt,pplay,zqsatco2)
781	do ig=1,ngrid
782	do l=1,nlay
783	satuco2(ig,l) = pq(ig,l,igcm_co2)*
784	& (mmean(ig,l)/44.01)*pplay(ig,l)/zqsatco2(ig,l)
785
786	c write(,) "In CO2 pt,sat ",pt(ig,l),satuco2(ig,l)
787	enddo
788	enddo
789	call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",3,
790	& satuco2)
791	call WRITEdiagfi(ngrid,"riceco2","ice radius","m"
792	& ,3,riceco2)
793	! or output in diagfi.nc (for testphys1d)
794	c call WRITEDIAGFI(ngrid,'ps','Surface pressure','Pa',0,ps)
795	c call WRITEDIAGFI(ngrid,'temp','Temperature ',
796	c & 'K JA',1,pt)
797
798	call WRITEdiagfi(ngrid,"rsedcloudco2","rsed co2","m",3,
799	& rsedcloudco2)
800
801	! used for rad. transfer calculations
802	! nuice is constant because a lognormal distribution is prescribed
803	c nuice(1:ngrid,1:nlay)=nuice_ref
804
805
806
807	c=======================================================================
808
809	END
810

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

Download in other formats:

Original Format