Context Navigation

co2cloud.F @ 1713

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

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

Download in other formats:

Original Format