Context Navigation

co2cloud.F @ 1657

Last change on this file since 1657 was 1655, checked in by aslmd, 8 years ago
corrected bugs from rev 1651 that prevented the mesoscale model from compiling. getin does not take double precision arguments. removed meteoriticflux_mod and included the variables in tracer_mod.
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	& 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	& +dble(meteo_flux_mass)
192	pq(ig,meteo_lvl,igcm_dust_number)=
193	& pq(ig,meteo_lvl,igcm_dust_number)
194	& +dble(meteo_flux_number)
195	enddo
196
197	call WRITEDIAGFI(ngrid,"pzlay","altitude","km",3,
198	& pzlay)
199	beta=0.85
200	subpdq(1:ngrid,1:nlay,1:nq) = 0
201	subpdt(1:ngrid,1:nlay) = 0
202	subpdqcloudco2(1:ngrid,1:nlay,1:nq) = 0
203	subpdtcloudco2(1:ngrid,1:nlay) = 0
204
205
206	wq(:,:)=0
207	! default value if no ice
208	rhocloudco2(1:ngrid,1:nlay) = rho_dust
209	rhocloudco2t(1:ngrid,1:nlay) = rho_dust
210	epaisseur(1:ngrid,1:nlay)=0
211	masse(1:ngrid,1:nlay)=0
212
213	tempo_traceur_t(1:ngrid,1:nlay)=0
214	tempo_traceurs(1:ngrid,1:nlay,1:nq)=0
215	sav_trac(1:ngrid,1:nlay,1:nq)=0
216	pdqsed(1:ngrid,1:nlay,1:nq)=0
217
218	do l=1,nlay
219	do ig=1, ngrid
220	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
221	epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
222
223	enddo
224	enddo
225
226
227
228
229
230
231
232
233	c-------------------------------------------------------------------
234	c 1. Tendencies:
235	c------------------
236
237
238
239	c------------------------------------------------------------------
240	c Time subsampling for microphysics
241	c------------------------------------------------------------------
242	DO microstep=1,imicro
243	c------ Temperature tendency subpdt
244	! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt
245	! If imicro=1 subpdt is the same as pdt
246	DO l=1,nlay
247	DO ig=1,ngrid
248	c tempo_traceur_t(ig,l)=tempo_traceur_t(ig,l)
249	c & + subpdtcloudco2(ig,l)
250	!write(,) 'T micro= ', tempo_traceur_t(ig,l)
251	c tempo_traceurs(ig,l,:)=tempo_traceurs(ig,l,:)
252	c & +subpdqcloudco2(ig,l,:)
253
254	subpdt(ig,l) = subpdt(ig,l)
255	& + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
256
257	subpdq(ig,l,igcm_dust_mass) =
258	& subpdq(ig,l,igcm_dust_mass)
259	& + pdq(ig,l,igcm_dust_mass)
260
261	subpdq(ig,l,igcm_dust_number) =
262	& subpdq(ig,l,igcm_dust_number)
263	& + pdq(ig,l,igcm_dust_number)
264
265	subpdq(ig,l,igcm_ccnco2_mass) =
266	& subpdq(ig,l,igcm_ccnco2_mass)
267	& + pdq(ig,l,igcm_ccnco2_mass)
268
269	subpdq(ig,l,igcm_ccnco2_number) =
270	& subpdq(ig,l,igcm_ccnco2_number)
271	& + pdq(ig,l,igcm_ccnco2_number)
272
273	subpdq(ig,l,igcm_co2_ice) =
274	& subpdq(ig,l,igcm_co2_ice)
275	& + pdq(ig,l,igcm_co2_ice)
276	subpdq(ig,l,igcm_co2) =
277	& subpdq(ig,l,igcm_co2)
278	& + pdq(ig,l,igcm_co2)
279
280	tempo_traceur_t(ig,l)= pt(ig,l)+subpdt(ig,l)*microtimestep
281	tempo_traceurs(ig,l,:)= pq(ig,l,:)+subpdq(ig,l,:)
282	& *microtimestep
283	!Stepped entry for sedimentation
284	ENDDO
285	ENDDO
286
287	!RSEDCLOUD AND RICECO2 HERE
288
289	DO l=1, nlay
290	DO ig=1,ngrid
291	Niceco2=max(tempo_traceurs(ig,l,igcm_co2_ice),1.e-30)
292	Nccnco2=max(tempo_traceurs(ig,l,igcm_ccnco2_number),
293	& 1.e-30)
294	Qccnco2=max(tempo_traceurs(ig,l,igcm_ccnco2_mass),
295	& 1.e-30)
296	mdustJA= tempo_traceurs(ig,l,igcm_dust_mass)
297	ndustJA=tempo_traceurs(ig,l,igcm_dust_number)
298	if ((ndustJA .lt. tauscaling(ig)) .or. (mdustJA .lt.
299	& 1.e-30 *tauscaling(ig))) then
300	rdust(ig,l)=1.e-10
301	else
302	rdust(ig,l)=(3./4./pi/2500.mdustJA/ndustJA)*(1./3.)
303	rdust(ig,l)=min(rdust(ig,l),5.e-6)
304	rdust(ig,l)=max(rdust(ig,l),1.e-9)
305	end if
306	c rhocloudco2t(ig,l) = (Niceco2 *rho_ice_co2
307	c & + Qccnco2*rho_dust)
308	c & / (Niceco2 + Qccnco2)
309	riceco2(ig,l)= Niceco2*3.0/
310	& (4.0rho_ice_co2pi*Nccnco2)
311	& +rdust(ig,l)rdust(ig,l)rdust(ig,l)
312	riceco2(ig,l)=riceco2(ig,l)**(1.0/3.0)
313	c write(,) "in co2clouds, rice = ",riceco2(ig,l)
314	c write(,) "in co2clouds, rho = ",rhocloudco2t(ig,l)
315
316	call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2,
317	& tauscaling(ig),riceco2(ig,l),rhocloudco2t(ig,l))
318	c write(,) "in co2clouds, rice update = ",riceco2(ig,l)
319	c write(,) "in co2clouds, rho update = "
320	c & ,rhocloudco2t(ig,l)
321
322	rsedcloudco2(ig,l)=max(riceco2(ig,l)*
323	& (1.+nuiceco2_sed)(1.+nuiceco2_sed)(1.+nuiceco2_sed),
324	& rdust(ig,l))
325	rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),5.e-4)
326	c write(,) 'Rsedcloud = ',rsedcloudco2(ig,l)
327	!write(,) 'Rhocloudco2 = ',rhocloudco2t(ig,l)
328
329	ENDDO
330	ENDDO
331
332	! Gravitational sedimentation
333
334	! sedimentation computed from radius computed from q in module radii_mod
335	sav_trac(:,:,igcm_co2_ice)=tempo_traceurs(:,:,igcm_co2_ice)
336	sav_trac(:,:,igcm_ccnco2_mass)=
337	& tempo_traceurs(:,:,igcm_ccnco2_mass)
338	sav_trac(:,:,igcm_ccnco2_number)=
339	& tempo_traceurs(:,:,igcm_ccnco2_number)
340
341	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
342	& microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
343	& rsedcloudco2,rhocloudco2t,
344	& tempo_traceurs(:,:,igcm_co2_ice),wq,beta) ! 3 traceurs
345
346	! sedim at the surface of co2 ice
347	do ig=1,ngrid
348	pdqs_sedco2(ig)=pdqs_sedco2(ig)+ wq(ig,1)
349	end do
350
351	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
352	& microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
353	& rsedcloudco2,rhocloudco2t,
354	& tempo_traceurs(:,:,igcm_ccnco2_mass),wq,beta)
355
356	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
357	& microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
358	& rsedcloudco2,rhocloudco2t,
359	& tempo_traceurs(:,:,igcm_ccnco2_number),wq,beta)
360
361
362	DO l = 1, nlay
363	DO ig=1,ngrid
364	pdqsed(ig,l,igcm_ccnco2_mass)=
365	& (tempo_traceurs(ig,l,igcm_ccnco2_mass)-
366	& sav_trac(ig,l,igcm_ccnco2_mass))/microtimestep
367	pdqsed(ig,l,igcm_ccnco2_number)=
368	& (tempo_traceurs(ig,l,igcm_ccnco2_number)-
369	& sav_trac(ig,l,igcm_ccnco2_number))/microtimestep
370	pdqsed(ig,l,igcm_co2_ice)=
371	& (tempo_traceurs(ig,l,igcm_co2_ice)-
372	& sav_trac(ig,l,igcm_co2_ice))/microtimestep
373	ENDDO
374	ENDDO
375	!pdqsed est la tendance due a la sedimentation
376
377	DO l = 1, nlay
378	DO ig=1,ngrid
379	pdqsed(ig,l,igcm_ccnco2_mass)=
380	& (tempo_traceurs(ig,l,igcm_ccnco2_mass)-
381	& sav_trac(ig,l,igcm_ccnco2_mass))/microtimestep
382	pdqsed(ig,l,igcm_ccnco2_number)=
383	& (tempo_traceurs(ig,l,igcm_ccnco2_number)-
384	& sav_trac(ig,l,igcm_ccnco2_number))/microtimestep
385	pdqsed(ig,l,igcm_co2_ice)=
386	& (tempo_traceurs(ig,l,igcm_co2_ice)-
387	& sav_trac(ig,l,igcm_co2_ice))/microtimestep
388	ENDDO
389	ENDDO
390	!pdqsed est la tendance due a la sedimentation
391	DO l=1,nlay
392	DO ig=1,ngrid
393	subpdq(ig,l,igcm_ccnco2_mass) =
394	& subpdq(ig,l,igcm_ccnco2_mass)
395	& +pdqsed(ig,l,igcm_ccnco2_mass)
396
397	subpdq(ig,l,igcm_ccnco2_number) =
398	& subpdq(ig,l,igcm_ccnco2_number)
399	& +pdqsed(ig,l,igcm_ccnco2_number)
400
401	subpdq(ig,l,igcm_co2_ice) =
402	& subpdq(ig,l,igcm_co2_ice)
403	& +pdqsed(ig,l,igcm_co2_ice)
404	ENDDO
405	ENDDO
406	c-------------------------------------------------------------------
407	c 2. Main call to the different cloud schemes:
408	c------------------------------------------------
409	IF (microphysco2) THEN
410	CALL improvedCO2clouds(ngrid,nlay,microtimestep,
411	& pplay,pt,subpdt,
412	& pq,subpdq,subpdqcloudco2,subpdtcloudco2,
413	& nq,tauscaling)
414
415	ELSE
416
417	write(,) ' no simpleCO2clouds procedure: STOP' ! listo
418	STOP
419
420	c CALL simpleclouds(ngrid,nlay,microtimestep, ! for water-ice clouds
421	c & pplay,pzlay,pt,subpdt,
422	c & pq,subpdq,subpdqcloud,subpdtcloud,
423	c & nq,tau,riceco2)
424	ENDIF
425
426
427	c-------------------------------------------------------------------
428	c 3. Updating tendencies after cloud scheme:
429	c-----------------------------------------------
430
431	c IF (microphysco2) THEN
432	DO l=1,nlay
433	DO ig=1,ngrid
434	subpdq(ig,l,igcm_dust_mass) =
435	& subpdq(ig,l,igcm_dust_mass)
436	& + subpdqcloudco2(ig,l,igcm_dust_mass)
437
438	subpdq(ig,l,igcm_dust_number) =
439	& subpdq(ig,l,igcm_dust_number)
440	& + subpdqcloudco2(ig,l,igcm_dust_number)
441
442	subpdq(ig,l,igcm_ccnco2_mass) =
443	& subpdq(ig,l,igcm_ccnco2_mass)
444	& + subpdqcloudco2(ig,l,igcm_ccnco2_mass)
445	c & +pdqsed(ig,l,igcm_ccnco2_mass)
446
447	subpdq(ig,l,igcm_ccnco2_number) =
448	& subpdq(ig,l,igcm_ccnco2_number)
449	& + subpdqcloudco2(ig,l,igcm_ccnco2_number)
450	c & +pdqsed(ig,l,igcm_ccnco2_number)
451
452	subpdq(ig,l,igcm_co2_ice) =
453	& subpdq(ig,l,igcm_co2_ice)
454	& + subpdqcloudco2(ig,l,igcm_co2_ice)
455	c & +pdqsed(ig,l,igcm_co2_ice)
456
457	subpdq(ig,l,igcm_co2) =
458	& subpdq(ig,l,igcm_co2)
459	& + subpdqcloudco2(ig,l,igcm_co2)
460	ENDDO
461	ENDDO
462
463
464	!ici
465	! call WRITEdiagfi(ngrid,"co2cloud000","co2 traceur","kg/kg",1,
466	! & pq(1,:,igcm_co2_ice) + ptimestep*
467	! & ( subpdq(1,:,igcm_co2_ice)))
468
469
470	IF (activice) THEN
471	DO l=1,nlay
472	DO ig=1,ngrid
473	subpdt(ig,l) =
474	& subpdt(ig,l) + subpdtcloudco2(ig,l)
475	ENDDO
476	ENDDO
477	ENDIF
478
479
480	ENDDO ! of DO microstep=1,imicro
481
482	c-------------------------------------------------------------------
483	c 6. Compute final tendencies after time loop:
484	c------------------------------------------------
485	c CO2 flux at surface (kg.m-2.s-1)
486	do ig=1,ngrid
487	pdqs_sedco2(ig)=pdqs_sedco2(ig)/ptimestep
488	enddo
489
490	c------ Temperature tendency pdtcloud
491	DO l=1,nlay
492	DO ig=1,ngrid
493	pdtcloudco2(ig,l) =
494	& subpdt(ig,l)/imicro-pdt(ig,l)
495	ENDDO
496	ENDDO
497
498	c------ Tracers tendencies pdqcloud
499	DO l=1,nlay
500	DO ig=1,ngrid
501
502	pdqcloudco2(ig,l,igcm_co2_ice) =
503	& subpdq(ig,l,igcm_co2_ice)/imicro
504	& - pdq(ig,l,igcm_co2_ice)
505	pdqcloudco2(ig,l,igcm_co2) =
506	& subpdq(ig,l,igcm_co2)/imicro
507	& - pdq(ig,l,igcm_co2)
508	ENDDO
509	ENDDO
510
511
512	! call WRITEdiagfi(ngrid,"co2cloud00","co2 traceur","kg/kg",1,
513	! & pq(1,:,igcm_co2_ice) + ptimestep*
514	! & (pdq(1,:,igcm_co2_ice) + pdqcloudco2(1,:,igcm_co2_ice)))
515
516
517	IF(microphysco2) THEN
518	DO l=1,nlay
519	DO ig=1,ngrid
520	pdqcloudco2(ig,l,igcm_ccnco2_mass) =
521	& subpdq(ig,l,igcm_ccnco2_mass)/imicro
522	& - pdq(ig,l,igcm_ccnco2_mass)
523	pdqcloudco2(ig,l,igcm_ccnco2_number) =
524	& subpdq(ig,l,igcm_ccnco2_number)/imicro
525	& - pdq(ig,l,igcm_ccnco2_number)
526	ENDDO
527	ENDDO
528	ENDIF
529
530
531	IF(scavenging) THEN
532	DO l=1,nlay
533	DO ig=1,ngrid
534	pdqcloudco2(ig,l,igcm_dust_mass) =
535	& subpdq(ig,l,igcm_dust_mass)/real(imicro)
536	& - pdq(ig,l,igcm_dust_mass)
537	pdqcloudco2(ig,l,igcm_dust_number) =
538	& subpdq(ig,l,igcm_dust_number)/real(imicro)
539	& - pdq(ig,l,igcm_dust_number)
540	ENDDO
541	ENDDO
542	ENDIF
543
544	c ENDIF
545	c------- Due to stepped entry, other processes tendencies can add up to negative values
546	c------- Therefore, enforce positive values and conserve mass
547
548
549	IF(microphysco2) THEN
550	DO l=1,nlay
551	DO ig=1,ngrid
552	IF ((pq(ig,l,igcm_ccnco2_number) +
553	& ptimestep* (pdq(ig,l,igcm_ccnco2_number) +
554	& pdqcloudco2(ig,l,igcm_ccnco2_number))
555	& .lt. 0)
556	& .or. (pq(ig,l,igcm_ccnco2_mass) +
557	& ptimestep* (pdq(ig,l,igcm_ccnco2_mass) +
558	& pdqcloudco2(ig,l,igcm_ccnco2_mass))
559	& .lt. 0)) THEN
560
561	pdqcloudco2(ig,l,igcm_ccnco2_number) =
562	& - pq(ig,l,igcm_ccnco2_number)/ptimestep
563	& - pdq(ig,l,igcm_ccnco2_number) +0
564
565	pdqcloudco2(ig,l,igcm_dust_number) =
566	& -pdqcloudco2(ig,l,igcm_ccnco2_number)
567
568	pdqcloudco2(ig,l,igcm_ccnco2_mass) =
569	& - pq(ig,l,igcm_ccnco2_mass)/ptimestep
570	& - pdq(ig,l,igcm_ccnco2_mass)+0
571
572	pdqcloudco2(ig,l,igcm_dust_mass) =
573	& -pdqcloudco2(ig,l,igcm_ccnco2_mass)
574
575	ENDIF
576	ENDDO
577	ENDDO
578	ENDIF
579
580
581	IF(scavenging) THEN
582	DO l=1,nlay
583	DO ig=1,ngrid
584	IF ( (pq(ig,l,igcm_dust_number) +
585	& ptimestep* (pdq(ig,l,igcm_dust_number) +
586	& pdqcloudco2(ig,l,igcm_dust_number)) .lt. 0.)
587	& .or. (pq(ig,l,igcm_dust_mass)+
588	& ptimestep* (pdq(ig,l,igcm_dust_mass) +
589	& pdqcloudco2(ig,l,igcm_dust_mass))
590	& .lt. 0.)) then
591
592	pdqcloudco2(ig,l,igcm_dust_number) =
593	& - pq(ig,l,igcm_dust_number)/ptimestep
594	& - pdq(ig,l,igcm_dust_number)+0
595
596	pdqcloudco2(ig,l,igcm_ccnco2_number) =
597	& -pdqcloudco2(ig,l,igcm_dust_number)
598
599	pdqcloudco2(ig,l,igcm_dust_mass) =
600	& - pq(ig,l,igcm_dust_mass)/ptimestep
601	& - pdq(ig,l,igcm_dust_mass) +0
602
603	pdqcloudco2(ig,l,igcm_ccnco2_mass) =
604	& -pdqcloudco2(ig,l,igcm_dust_mass)
605	ENDIF
606	ENDDO
607	ENDDO
608	ENDIF !pq+ptime*(pdq+pdqc)=1 ! pdqc=1-pq/ptime-pdq
609
610
611	DO l=1,nlay
612	DO ig=1,ngrid
613	IF (pq(ig,l,igcm_co2_ice) + ptimestep*
614	& (pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice))
615	& .lt. 1.e-25) THEN
616	pdqcloudco2(ig,l,igcm_co2_ice) =
617	& - pq(ig,l,igcm_co2_ice)/ptimestep - pdq(ig,l,igcm_co2_ice)
618	& +1.e-25
619	pdqcloudco2(ig,l,igcm_co2) = -pdqcloudco2(ig,l,igcm_co2_ice)
620	ENDIF
621	ENDDO
622	ENDDO
623
624
625
626
627	c------Update the ice and dust particle size "riceco2" for output or photochemistry
628	c------Only rsedcloudco2 is used for the co2 (cloud) cycle
629
630	IF(scavenging) THEN
631	DO l=1, nlay
632	DO ig=1,ngrid
633
634	c call updaterdust(
635	c & pq(ig,l,igcm_dust_mass) + ! dust mass
636	c & (pdq(ig,l,igcm_dust_mass) + ! dust mass
637	c & pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep, ! dust mass
638	c & pq(ig,l,igcm_dust_number) + ! dust number
639	c & (pdq(ig,l,igcm_dust_number) + ! dust number
640	c & pdqcloudco2(ig,l,igcm_dust_number))*ptimestep, ! dust number
641	c & rdust(ig,l))
642	c write(,) "in co2clouds, rdust(ig,l)= ",rdust(ig,l)
643	mdustJA= pq(ig,l,igcm_dust_mass) +
644	& (pdq(ig,l,igcm_dust_mass) +
645	& pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep
646	ndustJA=pq(ig,l,igcm_dust_number) +
647	& (pdq(ig,l,igcm_dust_number) +
648	& pdqcloudco2(ig,l,igcm_dust_number))*ptimestep
649	if ((ndustJA .lt. tauscaling(ig)) .or. (mdustJA .lt.
650	& 1.e-30 *tauscaling(ig))) then
651	rdust(ig,l)=1.e-10
652	else
653	rdust(ig,l)=(3./4./pi/2500.mdustJA/ndustJA)*(1./3.)
654	rdust(ig,l)=min(rdust(ig,l),5.e-4)
655	rdust(ig,l)=max(rdust(ig,l),1.e-10)
656	endif
657	ENDDO
658	ENDDO
659	ENDIF
660
661
662	IF(microphysco2) THEN
663
664	DO l=1, nlay
665	DO ig=1,ngrid
666
667	c call updaterice_microco2(
668	c & pq(ig,l,igcm_co2_ice) + ! ice mass
669	c & (pdq(ig,l,igcm_co2_ice) + ! ice mass
670	c & pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep, ! ice mass
671	c & pq(ig,l,igcm_ccnco2_mass) + ! ccn mass
672	c & (pdq(ig,l,igcm_ccnco2_mass) + ! ccn mass
673	c & pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep, ! ccn mass
674	c & pq(ig,l,igcm_ccnco2_number) + ! ccn number
675	c & (pdq(ig,l,igcm_ccnco2_number) + ! ccn number
676	c & pdqcloudco2(ig,l,igcm_ccnco2_number))*ptimestep, ! ccn number
677	c & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
678	c write(,) "in co2clouds, riceco2(ig,l)= ",riceco2(ig,l)
679
680
681	Niceco2=pq(ig,l,igcm_co2_ice) +
682	& (pdq(ig,l,igcm_co2_ice) +
683	& pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep
684	Nccnco2=max((pq(ig,l,igcm_ccnco2_number) +
685	& (pdq(ig,l,igcm_ccnco2_number) +
686	& pdqcloudco2(ig,l,igcm_ccnco2_number))ptimestep)
687	& tauscaling(ig),1.e-30)
688	Qccnco2=max((pq(ig,l,igcm_ccnco2_mass) +
689	& (pdq(ig,l,igcm_ccnco2_mass) +
690	& pdqcloudco2(ig,l,igcm_ccnco2_mass))ptimestep)
691	& tauscaling(ig),1.e-30)
692	c rhocloudco2t(ig,l) = (Niceco2 rho_ice_co2 + Qccnco2rho_dust)
693	c & / (Niceco2 + Qccnco2)
694	c rhocloudco2(ig,l) = min(max(rhocloudco2t,rho_ice_co2),rho_dust)
695
696	c write(,) "test, nccnco2 =",nccnco22
697
698
699	riceco2(ig,l)= Niceco2*3.0/
700	& (4.0rho_ice_co2pi*Nccnco2)
701	& +rdust(ig,l)rdust(ig,l)rdust(ig,l)
702
703	riceco2(ig,l)=riceco2(ig,l)**(1.0/3.0)
704	c write(,) "In co2cloud, after loop, riceco2 =",riceco2(ig,l)
705	c write(,) "In co2cloud, after loop, rhoco2 ="
706	c & ,rhocloudco2t(ig,l)
707
708	call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2,
709	& tauscaling(ig),riceco2(ig,l),rhocloudco2t(ig,l))
710
711	c write(,) "In co2cloud, after loop and update, riceco2 ="
712	c & ,riceco2(ig,l)
713	c write(,) "In co2cloud, after loop and update, rhoco2 ="
714	c & ,rhocloudco2t(ig,l)
715
716	if ( Niceco2
717	& .le. 1.e-23 .or. riceco2(ig,l) .le. 1.e-10 .or.
718	& riceco2(ig,l) .ge. 4.999e-4) then ! .or. riceco2(ig,l) .gt. 1.e-4 ) then
719	riceco2(ig,l)=0.
720
721	!NO CLOUD : RESET TRACER AND CONSERVE MASS
722	pdqcloudco2(ig,l,igcm_co2)= pq(ig,l,igcm_co2_ice)
723	& /ptimestep+pdq(ig,l,igcm_co2_ice)
724
725	pdqcloudco2(ig,l,igcm_co2_ice)=-pq(ig,l,igcm_co2_ice)
726	& /ptimestep-pdq(ig,l,igcm_co2_ice)
727
728	pdqcloudco2(ig,l,igcm_ccnco2_mass)=
729	& -pq(ig,l,igcm_ccnco2_mass)
730	& /ptimestep-pdq(ig,l,igcm_ccnco2_mass)
731
732	pdqcloudco2(ig,l,igcm_ccnco2_number)=
733	& -pq(ig,l,igcm_ccnco2_number)
734	& /ptimestep-pdq(ig,l,igcm_ccnco2_number)
735
736	pdqcloudco2(ig,l,igcm_dust_number)=
737	& pq(ig,l,igcm_ccnco2_number)
738	& /ptimestep+pdq(ig,l,igcm_ccnco2_number)
739
740	pdqcloudco2(ig,l,igcm_dust_mass)=
741	& pq(ig,l,igcm_ccnco2_mass)
742	& /ptimestep+pdq(ig,l,igcm_ccnco2_mass)
743	c$$$
744
745	c$$$
746	endif
747
748	c write(,) "in co2clouds, riceco2(ig,l)v2= ",riceco2(ig,l)
749
750	ENDDO
751	ENDDO
752
753	ELSE ! no microphys ! not of concern for co2 clouds - listo
754
755	ENDIF ! of IF(microphysco2)
756
757
758	c TO CHEK for relevancy - listo
759
760	c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
761	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
762	c Then that should not affect the ice particle radius
763	do ig=1,ngrid
764	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,2)))then
765	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,3)))
766	& riceco2(ig,2)=riceco2(ig,3)
767	riceco2(ig,1)=riceco2(ig,2)
768	end if
769	end do
770
771
772	DO l=1,nlay
773	DO ig=1,ngrid
774	rsedcloudco2(ig,l)=max(riceco2(ig,l)*
775	& (1.+nuiceco2_sed)(1.+nuiceco2_sed)(1.+nuiceco2_sed),
776	& rdust(ig,l))
777	rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),1.e-4)
778	ENDDO
779	ENDDO
780
781	call co2sat(ngrid*nlay,pt,pplay,zqsatco2)
782	do ig=1,ngrid
783	do l=1,nlay
784	satuco2(ig,l) = pq(ig,l,igcm_co2)*
785	& (mmean(ig,l)/44.01)*pplay(ig,l)/zqsatco2(ig,l)
786
787	c write(,) "In CO2 pt,sat ",pt(ig,l),satuco2(ig,l)
788	enddo
789	enddo
790	call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",3,
791	& satuco2)
792	call WRITEdiagfi(ngrid,"riceco2","ice radius","m"
793	& ,3,riceco2)
794	! or output in diagfi.nc (for testphys1d)
795	c call WRITEDIAGFI(ngrid,'ps','Surface pressure','Pa',0,ps)
796	c call WRITEDIAGFI(ngrid,'temp','Temperature ',
797	c & 'K JA',1,pt)
798
799	call WRITEdiagfi(ngrid,"rsedcloudco2","rsed co2","m",3,
800	& rsedcloudco2)
801
802	! used for rad. transfer calculations
803	! nuice is constant because a lognormal distribution is prescribed
804	c nuice(1:ngrid,1:nlay)=nuice_ref
805
806
807
808	c=======================================================================
809
810	END
811

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

Download in other formats:

Original Format