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

Rev	Line
[1617]	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,
[1655]	16	& rho_ice_co2,r3n_q,
	17	& meteo_flux_mass,meteo_flux_number,
	18	& meteo_alt
[1651]	19
[1617]	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
[1629]	87	DOUBLE PRECISION riceco2(ngrid,nlay) ! Ice mass mean radius (m)
[1617]	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
[1651]	130	double precision diff,diff0
	131	integer meteo_lvl
[1617]	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
[1651]	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
[1617]	170	firstcall=.false.
	171	ENDIF ! of IF (firstcall)
	172
	173	c-----Initialization
[1651]	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)
[1655]	191	& +dble(meteo_flux_mass)
[1651]	192	pq(ig,meteo_lvl,igcm_dust_number)=
[1655]	193	& pq(ig,meteo_lvl,igcm_dust_number)
	194	& +dble(meteo_flux_number)
[1651]	195	enddo
	196
	197	call WRITEDIAGFI(ngrid,"pzlay","altitude","km",3,
	198	& pzlay)
[1617]	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
[1649]	291	Niceco2=max(tempo_traceurs(ig,l,igcm_co2_ice),1.e-30)
[1617]	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
[1649]	306	c rhocloudco2t(ig,l) = (Niceco2 *rho_ice_co2
	307	c & + Qccnco2*rho_dust)
	308	c & / (Niceco2 + Qccnco2)
[1617]	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)
[1649]	313	c write(,) "in co2clouds, rice = ",riceco2(ig,l)
	314	c write(,) "in co2clouds, rho = ",rhocloudco2t(ig,l)
[1617]	315
[1649]	316	call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2,
[1617]	317	& tauscaling(ig),riceco2(ig,l),rhocloudco2t(ig,l))
[1649]	318	c write(,) "in co2clouds, rice update = ",riceco2(ig,l)
	319	c write(,) "in co2clouds, rho update = "
	320	c & ,rhocloudco2t(ig,l)
[1617]	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)
[1649]	326	c write(,) 'Rsedcloud = ',rsedcloudco2(ig,l)
[1617]	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)
[1649]	692	c rhocloudco2t(ig,l) = (Niceco2 rho_ice_co2 + Qccnco2rho_dust)
	693	c & / (Niceco2 + Qccnco2)
[1617]	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)
[1649]	704	c write(,) "In co2cloud, after loop, riceco2 =",riceco2(ig,l)
	705	c write(,) "In co2cloud, after loop, rhoco2 ="
	706	c & ,rhocloudco2t(ig,l)
[1617]	707
	708	call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2,
	709	& tauscaling(ig),riceco2(ig,l),rhocloudco2t(ig,l))
	710
[1649]	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)
[1617]	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
[1649]	787	c write(,) "In CO2 pt,sat ",pt(ig,l),satuco2(ig,l)
[1617]	788	enddo
	789	enddo
[1649]	790	call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",3,
[1629]	791	& satuco2)
	792	call WRITEdiagfi(ngrid,"riceco2","ice radius","m"
[1649]	793	& ,3,riceco2)
[1617]	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
[1649]	799	call WRITEdiagfi(ngrid,"rsedcloudco2","rsed co2","m",3,
[1617]	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