Context Navigation

callsedim_mod.F @ 2596

Last change on this file since 2596 was 2589, checked in by cmathe, 4 years ago
GCM Mars: add meteoritic particles as CN for CO2 microphysics
File size: 27.8 KB

Rev	Line
[1962]	1	MODULE callsedim_mod
	2
	3	IMPLICIT NONE
	4
	5	CONTAINS
[1974]	6
	7	SUBROUTINE callsedim(ngrid,nlay,ptimestep,
[2199]	8	& pplev,zlev,zlay,pt,pdt,
	9	& rdust,rstormdust,rtopdust,
	10	& rice,rsedcloud,rhocloud,
[1974]	11	& pq,pdqfi,pdqsed,pdqs_sed,nq,
[740]	12	& tau,tauscaling)
[1974]	13
[2304]	14	USE ioipsl_getin_p_mod, only: getin_p
[1036]	15	USE updaterad, only: updaterdust,updaterice_micro,updaterice_typ
	16	USE tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number,
	17	& rho_dust, rho_q, radius, varian,
	18	& igcm_ccn_mass, igcm_ccn_number,
[2312]	19	& igcm_h2o_ice, igcm_hdo_ice,
	20	& nuice_sed, nuice_ref,
[2589]	21	& igcm_co2_ice, igcm_stormdust_mass,
[2199]	22	& igcm_stormdust_number,igcm_topdust_mass,
[2322]	23	& igcm_topdust_number,
[2332]	24	& nqfils,qperemin,masseqmin ! MVals: variables isotopes
[1913]	25	USE newsedim_mod, ONLY: newsedim
	26	USE comcstfi_h, ONLY: g
[1983]	27	USE dimradmars_mod, only: naerkind
[2409]	28	USE dust_param_mod, ONLY: doubleq
[38]	29	IMPLICIT NONE
	30
	31	c=======================================================================
	32	c Sedimentation of the Martian aerosols
	33	c depending on their density and radius
	34	c
	35	c F.Forget 1999
	36	c
	37	c Modified by J.-B. Madeleine 2010: Now includes the doubleq
	38	c technique in order to have only one call to callsedim in
	39	c physiq.F.
	40	c
[1617]	41	c Modified by J. Audouard 09/16: Now includes the co2clouds case
	42	c If the co2 microphysics is on, then co2 theice & ccn tracers
	43	c are being sedimented in the microtimestep (co2cloud.F), not
	44	c in this routine.
	45	c
[38]	46	c=======================================================================
	47
	48	c-----------------------------------------------------------------------
	49	c declarations:
	50	c -------------
[520]	51
[1913]	52	include "callkeys.h"
[38]	53
	54	c
	55	c arguments:
	56	c ----------
	57
[1005]	58	integer,intent(in) :: ngrid ! number of horizontal grid points
	59	integer,intent(in) :: nlay ! number of atmospheric layers
	60	real,intent(in) :: ptimestep ! physics time step (s)
	61	real,intent(in) :: pplev(ngrid,nlay+1) ! pressure at inter-layers (Pa)
	62	real,intent(in) :: zlev(ngrid,nlay+1) ! altitude at layer boundaries
	63	real,intent(in) :: zlay(ngrid,nlay) ! altitude at the middle of the layers
	64	real,intent(in) :: pt(ngrid,nlay) ! temperature at mid-layer (K)
	65	real,intent(in) :: pdt(ngrid,nlay) ! tendency on temperature, from
	66	! previous processes (K/s)
[38]	67	c Aerosol radius provided by the water ice microphysical scheme:
[1005]	68	real,intent(out) :: rdust(ngrid,nlay) ! Dust geometric mean radius (m)
[1974]	69	real,intent(out) :: rstormdust(ngrid,nlay) ! Stormdust geometric mean radius (m)
[2199]	70	real,intent(out) :: rtopdust(ngrid,nlay) ! topdust geometric mean radius (m)
[1005]	71	real,intent(out) :: rice(ngrid,nlay) ! H2O Ice geometric mean radius (m)
	72	c Sedimentation radius of water ice
[1047]	73	real,intent(in) :: rsedcloud(ngrid,nlay)
[1005]	74	c Cloud density (kg.m-3)
[1047]	75	real,intent(inout) :: rhocloud(ngrid,nlay)
[38]	76	c Traceurs :
[1005]	77	real,intent(in) :: pq(ngrid,nlay,nq) ! tracers (kg/kg)
	78	real,intent(in) :: pdqfi(ngrid,nlay,nq) ! tendency before sedimentation (kg/kg.s-1)
	79	real,intent(out) :: pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1)
	80	real,intent(out) :: pdqs_sed(ngrid,nq) ! flux at surface (kg.m-2.s-1)
	81	integer,intent(in) :: nq ! number of tracers
[1983]	82	real,intent(in) :: tau(ngrid,naerkind) ! dust opacity
[1005]	83	real,intent(in) :: tauscaling(ngrid)
[38]	84
	85	c local:
	86	c ------
	87
	88	INTEGER l,ig, iq
	89
[1005]	90	real zqi(ngrid,nlay,nq) ! to locally store tracers
	91	real zt(ngrid,nlay) ! to locally store temperature
	92	real masse (ngrid,nlay) ! Layer mass (kg.m-2)
	93	real epaisseur (ngrid,nlay) ! Layer thickness (m)
[2322]	94	real wq(ngrid,nlay+1),w(ngrid,nlay) ! displaced tracer mass wq (kg.m-2), MVals: displaced "pere" tracer mass w (kg.m-2)
[1005]	95	real r0(ngrid,nlay) ! geometric mean radius used for
[358]	96	! sedimentation (m)
[1005]	97	real r0dust(ngrid,nlay) ! geometric mean radius used for
[358]	98	! dust (m)
[1974]	99	real r0stormdust(ngrid,nlay) ! Geometric mean radius used for stormdust (m)
[740]	100	! ! CCNs (m)
[2199]	101	real r0topdust(ngrid,nlay) ! Geometric mean radius used for topdust (m)
	102	! ! CCNs (m)
[1005]	103	real,save :: beta ! correction for the shape of the ice particles (cf. newsedim)
[411]	104	c for ice radius computation
	105	REAL Mo,No
	106	REAl ccntyp
[2304]	107	character(len=20),parameter :: modname="callsedim"
[2322]	108	c MVals: transport of the isotopic ratio
	109	REAL Ratio0(ngrid,nlay),Ratio(ngrid,nlay)
	110	REAL masseq(ngrid,nlay)
	111	REAL newmasse
	112	REAL zq0(ngrid,nlay,nq)
	113	INTEGER ifils,iq2
[38]	114
[358]	115
[38]	116	c Discrete size distributions (doubleq)
	117	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	118	c 1) Parameters used to represent the changes in fall
	119	c velocity as a function of particle size;
[1005]	120	integer ir
	121	integer,parameter :: nr=12 !(nr=7) ! number of bins
	122	real,save :: rd(nr)
	123	real qr(ngrid,nlay,nr)
	124	real,save :: rdi(nr+1) ! extreme and intermediate radii
	125	real Sq(ngrid,nlay)
	126	real,parameter :: rdmin=1.e-8
	127	real,parameter :: rdmax=30.e-6
	128	real,parameter :: rdimin=1.e-8 ! 1.e-7
	129	real,parameter :: rdimax=1.e-4
[38]	130
	131	c 2) Second size distribution for the log-normal integration
	132	c (the mass mixing ratio is computed for each radius)
	133
[1005]	134	integer iint
	135	integer,parameter :: ninter=4 ! number of points between each rdi radii
	136	real,save :: rr(ninter,nr)
[38]	137	integer radpower
[358]	138	real sigma0
[38]	139
	140	c 3) Other local variables used in doubleq
	141
[1005]	142	INTEGER,SAVE :: idust_mass ! index of tracer containing dust mass
	143	! mix. ratio
	144	INTEGER,SAVE :: idust_number ! index of tracer containing dust number
	145	! mix. ratio
	146	INTEGER,SAVE :: iccn_mass ! index of tracer containing CCN mass
	147	! mix. ratio
	148	INTEGER,SAVE :: iccn_number ! index of tracer containing CCN number
	149	! mix. ratio
[1974]	150	INTEGER,SAVE :: istormdust_mass ! index of tracer containing
	151	!stormdust mass mix. ratio
	152	INTEGER,SAVE :: istormdust_number ! index of tracer containing
[2199]	153	!stormdust number mix. ratio
	154	INTEGER,SAVE :: itopdust_mass ! index of tracer containing
	155	!topdust mass mix. ratio
	156	INTEGER,SAVE :: itopdust_number ! index of tracer containing
	157	!topdust number mix. ratio
[1617]	158	INTEGER,SAVE :: iccnco2_number ! index of tracer containing CCN number
	159	INTEGER,SAVE :: iccnco2_mass ! index of tracer containing CCN number
[2562]	160	INTEGER,SAVE :: iccnco2_h2o_number ! index of tracer containing CCN number
	161	INTEGER,SAVE :: iccnco2_h2o_mass_ice ! index of tracer containing CCN number
	162	INTEGER,SAVE :: iccnco2_h2o_mass_ccn ! index of tracer containing CCN number
[2589]	163	INTEGER,SAVE :: iccnco2_meteor_number ! index of tracer containing CCN number
	164	INTEGER,SAVE :: iccnco2_meteor_mass ! index of tracer containing CCN number
[1617]	165	INTEGER,SAVE :: ico2_ice ! index of tracer containing CCN number
[38]	166
[1617]	167
[1005]	168	LOGICAL,SAVE :: firstcall=.true.
[38]	169
[1974]	170
	171
[38]	172	c ** un petit test de coherence
	173	c --------------------------
[1779]	174	! AS: firstcall OK absolute
[38]	175	IF (firstcall) THEN
[411]	176
[38]	177	c Doubleq: initialization
	178	IF (doubleq) THEN
	179	do ir=1,nr
	180	rd(ir)= rdmin(rdmax/rdmin)*(float(ir-1)/float(nr-1))
	181	end do
	182	rdi(1)=rdimin
	183	do ir=2,nr
	184	rdi(ir)= sqrt(rd(ir-1)*rd(ir))
	185	end do
	186	rdi(nr+1)=rdimax
	187
	188	do ir=1,nr
	189	do iint=1,ninter
	190	rr(iint,ir)=
	191	& rdi(ir)*
	192	& (rdi(ir+1)/rdi(ir))**(float(iint-1)/float(ninter-1))
	193	c write(,) rr(iint,ir)
	194	end do
	195	end do
	196
	197	! identify tracers corresponding to mass mixing ratio and
	198	! number mixing ratio
	199	idust_mass=0 ! dummy initialization
	200	idust_number=0 ! dummy initialization
	201
	202	do iq=1,nq
	203	if (noms(iq).eq."dust_mass") then
	204	idust_mass=iq
[1005]	205	write(,)"callsedim: idust_mass=",idust_mass
[38]	206	endif
	207	if (noms(iq).eq."dust_number") then
	208	idust_number=iq
[1005]	209	write(,)"callsedim: idust_number=",idust_number
[38]	210	endif
	211	enddo
	212
	213	! check that we did find the tracers
	214	if ((idust_mass.eq.0).or.(idust_number.eq.0)) then
	215	write(,) 'callsedim: error! could not identify'
	216	write(,) ' tracers for dust mass and number mixing'
	217	write(,) ' ratio and doubleq is activated!'
[2304]	218	call abort_physic(modname,"missing dust tracers",1)
[38]	219	endif
	220	ENDIF !of if (doubleq)
	221
[740]	222	IF (microphys) THEN
[358]	223	iccn_mass=0
	224	iccn_number=0
	225	do iq=1,nq
	226	if (noms(iq).eq."ccn_mass") then
	227	iccn_mass=iq
[1005]	228	write(,)"callsedim: iccn_mass=",iccn_mass
[358]	229	endif
	230	if (noms(iq).eq."ccn_number") then
	231	iccn_number=iq
[1005]	232	write(,)"callsedim: iccn_number=",iccn_number
[358]	233	endif
	234	enddo
	235	! check that we did find the tracers
	236	if ((iccn_mass.eq.0).or.(iccn_number.eq.0)) then
	237	write(,) 'callsedim: error! could not identify'
	238	write(,) ' tracers for ccn mass and number mixing'
[740]	239	write(,) ' ratio and microphys is activated!'
[2304]	240	call abort_physic(modname,"missing ccn tracers",1)
[358]	241	endif
[740]	242	ENDIF !of if (microphys)
[358]	243
[1720]	244	IF (co2clouds) THEN
[1617]	245	iccnco2_mass=0
	246	iccnco2_number=0
[2562]	247	iccnco2_h2o_mass_ice=0
	248	iccnco2_h2o_mass_ccn=0
	249	iccnco2_h2o_number=0
[2589]	250	iccnco2_meteor_mass=0
	251	iccnco2_meteor_number=0
[1617]	252	ico2_ice=0
	253	do iq=1,nq
	254	if (noms(iq).eq."ccnco2_mass") then
	255	iccnco2_mass=iq
	256	write(,)"callsedim: iccnco2_mass=",iccnco2_mass
	257	endif
	258	if (noms(iq).eq."co2_ice") then
	259	ico2_ice=iq
	260	write(,)"callsedim: ico2_ice=",ico2_ice
	261	endif
	262	if (noms(iq).eq."ccnco2_number") then
	263	iccnco2_number=iq
	264	write(,)"callsedim: iccnco2_number=",iccnco2_number
	265	endif
[2562]	266	if (noms(iq).eq."ccnco2_h2o_number") then
	267	iccnco2_h2o_number=iq
	268	write(,)"callsedim: iccnco2_h2o_number=",
	269	& iccnco2_h2o_number
	270	endif
	271	if (noms(iq).eq."ccnco2_h2o_mass_ice") then
	272	iccnco2_h2o_mass_ice=iq
	273	write(,)"callsedim: iccnco2_h2o_mass_ice=",
	274	& iccnco2_h2o_mass_ice
	275	endif
	276	if (noms(iq).eq."ccnco2_h2o_mass_ccn") then
	277	iccnco2_h2o_mass_ccn=iq
	278	write(,)"callsedim: iccnco2_h2o_mass_ccn=",
	279	& iccnco2_h2o_mass_ccn
	280	endif
[2589]	281	if (noms(iq).eq."ccnco2_meteor_number") then
	282	iccnco2_meteor_number=iq
	283	write(,)"callsedim: iccnco2_meteor_number=",
	284	& iccnco2_meteor_number
	285	endif
	286	if (noms(iq).eq."ccnco2_meteor_mass") then
	287	iccnco2_meteor_mass=iq
	288	write(,)"callsedim: iccnco2_meteor_mass=",
	289	& iccnco2_meteor_mass
	290	endif
[1617]	291	enddo
	292	! check that we did find the tracers
	293	if ((iccnco2_mass.eq.0).or.(iccnco2_number.eq.0)) then
	294	write(,) 'callsedim: error! could not identify'
	295	write(,) ' tracers for ccn co2 mass and number mixing'
[1720]	296	write(,) ' ratio and co2clouds are activated!'
[2304]	297	call abort_physic(modname,"missing co2 ccn tracers",1)
[1617]	298	endif
[2562]	299	if (co2useh2o) then
	300	if((iccnco2_h2o_mass_ccn.eq.0).or.
	301	& (iccnco2_h2o_mass_ice.eq.0).or.
	302	& (iccnco2_h2o_number.eq.0)) then
	303	write(,) 'callsedim: error! could not identify'
	304	write(,) ' tracers for ccn co2 mass and number mixing'
	305	write(,) ' ratio and co2clouds are activated!'
	306	call abort_physic(modname,"missing co2 ccn tracers",1)
	307	end if
	308	end if
[2589]	309	if (meteo_flux) then
	310	if((iccnco2_meteor_mass.eq.0).or.
	311	& (iccnco2_meteor_number.eq.0)) then
	312	write(,) 'callsedim: error! could not identify'
	313	write(,) ' tracers for ccn co2 mass and number mixing'
	314	write(,) ' ratio and co2clouds are activated!'
	315	call abort_physic(modname,"missing co2 ccn tracers",1)
	316	end if
	317	end if
	318
[1720]	319	ENDIF !of if (co2clouds)
[1617]	320
[1974]	321	IF (water) THEN
[626]	322	write(,) "correction for the shape of the ice particles ?"
	323	beta=0.75 ! default value
[2304]	324	call getin_p("ice_shape",beta)
[626]	325	write(,) " ice_shape = ",beta
	326
	327	write(,) "water_param nueff Sedimentation:", nuice_sed
	328	IF (activice) THEN
	329	write(,) "water_param nueff Radiative:", nuice_ref
	330	ENDIF
[1974]	331	ENDIF
	332
	333	IF (rdstorm) THEN ! identifying stormdust tracers for sedimentation
	334	istormdust_mass=0 ! dummy initialization
	335	istormdust_number=0 ! dummy initialization
	336
	337	do iq=1,nq
	338	if (noms(iq).eq."stormdust_mass") then
	339	istormdust_mass=iq
	340	write(,)"callsedim: istormdust_mass=",istormdust_mass
	341	endif
	342	if (noms(iq).eq."stormdust_number") then
	343	istormdust_number=iq
	344	write(,)"callsedim: istormdust_number=",
	345	& istormdust_number
	346	endif
	347	enddo
	348
	349	! check that we did find the tracers
	350	if ((istormdust_mass.eq.0).or.(istormdust_number.eq.0)) then
	351	write(,) 'callsedim: error! could not identify'
	352	write(,) ' tracers for stormdust mass and number mixing'
	353	write(,) ' ratio and rdstorm is activated!'
[2304]	354	call abort_physic(modname,"missing stormdust tracers",1)
[1974]	355	endif
	356	ENDIF !of if (rdstorm)
[2199]	357
	358	IF (slpwind) THEN ! identifying topdust tracers for sedimentation
	359	itopdust_mass=0 ! dummy initialization
	360	itopdust_number=0 ! dummy initialization
	361
	362	do iq=1,nq
	363	if (noms(iq).eq."topdust_mass") then
	364	itopdust_mass=iq
	365	write(,)"callsedim: itopdust_mass=",itopdust_mass
	366	endif
	367	if (noms(iq).eq."topdust_number") then
	368	itopdust_number=iq
	369	write(,)"callsedim: itopdust_number=",
	370	& itopdust_number
	371	endif
	372	enddo
	373
	374	! check that we did find the tracers
	375	if ((itopdust_mass.eq.0).or.(itopdust_number.eq.0)) then
	376	write(,) 'callsedim: error! could not identify'
	377	write(,) ' tracers for topdust mass and number mixing'
	378	write(,) ' ratio and slpwind is activated!'
[2304]	379	call abort_physic(modname,"missing topdust tracers",1)
[2199]	380	endif
	381	ENDIF !of if (slpwind)
	382
[38]	383	firstcall=.false.
	384	ENDIF ! of IF (firstcall)
	385
	386	c-----------------------------------------------------------------------
	387	c 1. Initialization
	388	c -----------------
	389
[1005]	390	! zqi(1:ngrid,1:nlay,1:nqmx) = 0.
	391	c Update the mass mixing ratio and temperature with the tendencies coming
[38]	392	c from other parameterizations:
	393	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1005]	394	zqi(1:ngrid,1:nlay,1:nq)=pq(1:ngrid,1:nlay,1:nq)
	395	& +pdqfi(1:ngrid,1:nlay,1:nq)*ptimestep
[2323]	396	zq0(1:ngrid,1:nlay,1:nq)=pq(1:ngrid,1:nlay,1:nq) !MVals: keep the input value
	397	& +pdqfi(1:ngrid,1:nlay,1:nq)*ptimestep
[1005]	398	zt(1:ngrid,1:nlay)=pt(1:ngrid,1:nlay)
	399	& +pdt(1:ngrid,1:nlay)*ptimestep
[38]	400
	401	c Computing the different layer properties
	402	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	403	c Mass (kg.m-2), thickness(m), crossing time (s) etc.
	404
	405	do l=1,nlay
	406	do ig=1, ngrid
	407	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
	408	epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
	409	end do
	410	end do
	411
	412	c =================================================================
[358]	413	c Compute the geometric mean radius used for sedimentation
	414
	415	if (doubleq) then
	416	do l=1,nlay
	417	do ig=1, ngrid
[740]	418
	419	call updaterdust(zqi(ig,l,igcm_dust_mass),
	420	& zqi(ig,l,igcm_dust_number),r0dust(ig,l),
	421	& tauscaling(ig))
	422
[358]	423	end do
	424	end do
	425	endif
[2199]	426	c rocket dust storm
[1974]	427	if (rdstorm) then
	428	do l=1,nlay
	429	do ig=1, ngrid
	430
	431	call updaterdust(zqi(ig,l,igcm_stormdust_mass),
	432	& zqi(ig,l,igcm_stormdust_number),r0stormdust(ig,l),
	433	& tauscaling(ig))
	434
	435	end do
	436	end do
	437	endif
[2199]	438	c entrainment by slope wind
	439	if (slpwind) then
	440	do l=1,nlay
	441	do ig=1, ngrid
	442
	443	call updaterdust(zqi(ig,l,igcm_topdust_mass),
	444	& zqi(ig,l,igcm_topdust_number),r0topdust(ig,l),
	445	& tauscaling(ig))
	446
	447	end do
	448	end do
	449	endif
[358]	450	c =================================================================
[38]	451	do iq=1,nq
[1617]	452	if(radius(iq).gt.1.e-9 .and.(iq.ne.ico2_ice) .and.
	453	& (iq .ne. iccnco2_mass) .and. (iq .ne.
[2562]	454	& iccnco2_number) .and.
	455	& (iq.ne.iccnco2_h2o_number).and.
	456	& (iq.ne.iccnco2_h2o_mass_ice).and.
[2589]	457	& (iq.ne.iccnco2_h2o_mass_ccn).and.
	458	& (iq.ne.iccnco2_meteor_mass).and.
	459	& (iq.ne.iccnco2_meteor_number).and. ! no sedim for gaz or CO2 clouds (done in microtimestep)
[2322]	460	& iq .ne. igcm_hdo_ice) then !MVals: hdo is transported by h2o
[38]	461	c -----------------------------------------------------------------
	462	c DOUBLEQ CASE
	463	c -----------------------------------------------------------------
[2199]	464	if ( doubleq.and.
[1974]	465	& ((iq.eq.idust_mass).or.(iq.eq.idust_number).or.
[2199]	466	& (iq.eq.istormdust_mass).or.(iq.eq.istormdust_number).or.
	467	& (iq.eq.itopdust_mass).or.(iq.eq.itopdust_number)) ) then
[740]	468
[38]	469	c Computing size distribution:
	470	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	471
[1974]	472	if ((iq.eq.idust_mass).or.(iq.eq.idust_number)) then
[358]	473	do l=1,nlay
	474	do ig=1, ngrid
	475	r0(ig,l)=r0dust(ig,l)
	476	end do
[38]	477	end do
[1974]	478	else if ((iq.eq.istormdust_mass).or.
	479	& (iq.eq.istormdust_number)) then
	480	do l=1,nlay
	481	do ig=1, ngrid
	482	r0(ig,l)=r0stormdust(ig,l)
	483	end do
	484	end do
[2199]	485	else if ((iq.eq.itopdust_mass).or.
	486	& (iq.eq.itopdust_number)) then
	487	do l=1,nlay
	488	do ig=1, ngrid
	489	r0(ig,l)=r0topdust(ig,l)
	490	end do
	491	end do
[1974]	492	endif
	493	sigma0 = varian
[38]	494
	495	c Computing mass mixing ratio for each particle size
	496	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2199]	497	IF ((iq.EQ.idust_mass).or.(iq.EQ.istormdust_mass)
	498	& .or.(iq.EQ.itopdust_mass)) then
[38]	499	radpower = 2
[520]	500	ELSE ! number
[38]	501	radpower = -1
	502	ENDIF
	503	Sq(1:ngrid,1:nlay) = 0.
	504	do ir=1,nr
	505	do l=1,nlay
	506	do ig=1,ngrid
	507	c ****************
	508	c Size distribution integration
	509	c (Trapezoid Integration Method)
	510	qr(ig,l,ir)=0.5(rr(2,ir)-rr(1,ir))
	511	& (rr(1,ir)*radpower)
[358]	512	& exp(-(log(rr(1,ir)/r0(ig,l)))*2/(2sigma0**2))
[38]	513	do iint=2,ninter-1
	514	qr(ig,l,ir)=qr(ig,l,ir) +
	515	& 0.5(rr(iint+1,ir)-rr(iint-1,ir))
	516	& (rr(iint,ir)*radpower)
	517	& exp(-(log(rr(iint,ir)/r0(ig,l)))**2/
[358]	518	& (2sigma0*2))
[38]	519	end do
	520	qr(ig,l,ir)=qr(ig,l,ir) +
	521	& 0.5(rr(ninter,ir)-rr(ninter-1,ir))
	522	& (rr(ninter,ir)*radpower)
	523	& exp(-(log(rr(ninter,ir)/r0(ig,l)))**2/
[358]	524	& (2sigma0*2))
[38]	525
	526	c **************** old method (not recommended!)
	527	c qr(ig,l,ir)=(rd(ir)*(5-3iq))*
[358]	528	c & exp( -(log(rd(ir)/r0(ig,l)))*2 / (2sigma0**2) )
[38]	529	c ******************************
	530
	531	Sq(ig,l)=Sq(ig,l)+qr(ig,l,ir)
	532	enddo
	533	enddo
	534	enddo
	535
	536	do ir=1,nr
	537	do l=1,nlay
	538	do ig=1,ngrid
	539	qr(ig,l,ir) = zqi(ig,l,iq)*qr(ig,l,ir)/Sq(ig,l)
	540	enddo
	541	enddo
	542	enddo
	543
	544	c Computing sedimentation for each tracer
	545	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	546
	547	zqi(1:ngrid,1:nlay,iq) = 0.
	548	pdqs_sed(1:ngrid,iq) = 0.
	549
	550	do ir=1,nr
[358]	551	call newsedim(ngrid,nlay,1,1,ptimestep,
[1913]	552	& pplev,masse,epaisseur,zt,rd(ir),(/rho_dust/),qr(1,1,ir),
[358]	553	& wq,0.5)
[38]	554
	555	c Tendencies
	556	c ~~~~~~~~~~
	557	do ig=1,ngrid
	558	pdqs_sed(ig,iq) = pdqs_sed(ig,iq)
	559	& + wq(ig,1)/ptimestep
	560	end do
	561	DO l = 1, nlay
	562	DO ig=1,ngrid
	563	zqi(ig,l,iq)=zqi(ig,l,iq)+qr(ig,l,ir)
	564	ENDDO
[1974]	565	ENDDO
[38]	566	enddo ! of do ir=1,nr
	567	c -----------------------------------------------------------------
[626]	568	c WATER CYCLE CASE
[38]	569	c -----------------------------------------------------------------
[626]	570	else if ((iq .eq. iccn_mass) .or. (iq .eq. iccn_number)
[2322]	571	& .or. (iq .eq. igcm_h2o_ice)) then
[626]	572	if (microphys) then
[2322]	573	c water ice sedimentation
	574	c ~~~~~~~~~~
[626]	575	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
[1005]	576	& ptimestep,pplev,masse,epaisseur,zt,rsedcloud,rhocloud,
[626]	577	& zqi(1,1,iq),wq,beta)
	578	else
[2322]	579	c water ice sedimentation
	580	c ~~~~~~~~~~
[626]	581	call newsedim(ngrid,nlay,ngrid*nlay,1,
[1005]	582	& ptimestep,pplev,masse,epaisseur,zt,rsedcloud,rho_q(iq),
[626]	583	& zqi(1,1,iq),wq,beta)
[2407]	584	endif ! of if (microphys)
	585	c Surface tendencies
	586	c ~~~~~~~~~~
	587	do ig=1,ngrid
	588	pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
	589	end do
	590	c Special case of isotopes
	591	c ~~~~~~~~~~
	592	!MVals: for now only water can have an isotope, a son ("fils"), and it has to be hdo
	593	if (nqfils(iq).gt.0) then
	594	if (iq.eq.igcm_h2o_ice) then
	595	iq2=igcm_hdo_ice
	596	else
	597	call abort_physic("callsedim_mod","invalid isotope",1)
	598	endif
	599	!MVals: input parameters in vlz_fi for hdo
	600	do l=1,nlay
	601	do ig=1,ngrid
	602	if (zq0(ig,l,iq).gt.qperemin) then
	603	Ratio0(ig,l)=zq0(ig,l,iq2)/zq0(ig,l,iq)
	604	else
	605	Ratio0(ig,l)=0.
	606	endif
	607	Ratio(ig,l)=Ratio0(ig,l)
	608	masseq(ig,l)=max(masse(ig,l)*zq0(ig,l,iq),masseqmin)
	609	w(ig,l)=wq(ig,l) !MVals: very important: hdo is transported by h2o (see vlsplt_p.F: correction bugg 15 mai 2015)
	610	enddo !ig=1,ngrid
	611	enddo !l=1,nlay
	612	!MVals: no need to enter newsedim as the transporting mass w has been already calculated
	613	call vlz_fi(ngrid,nlay,Ratio,2.,masseq,w,wq)
	614	zqi(:,nlay,iq2)=zqi(:,nlay,iq)*Ratio0(:,nlay)
	615	do l=1,nlay-1
	616	do ig=1,ngrid
	617	newmasse=max(masseq(ig,l)+w(ig,l+1)-w(ig,l),masseqmin)
	618	Ratio(ig,l)=(Ratio0(ig,l)*masseq(ig,l)
	619	& +wq(ig,l+1)-wq(ig,l))/newmasse
	620	zqi(ig,l,iq2)=zqi(ig,l,iq)*Ratio(ig,l)
	621	enddo
	622	enddo !l=1,nlay-1
	623	!MVals: hdo surface tendency
	624	do ig=1,ngrid
	625	if (w(ig,1).gt.masseqmin) then
	626	pdqs_sed(ig,iq2)=pdqs_sed(ig,iq)*(wq(ig,1)/w(ig,1))
	627	else
	628	pdqs_sed(ig,iq2)=pdqs_sed(ig,iq)*Ratio0(ig,1)
	629	endif
[2322]	630	end do
[2407]	631	endif !(nqfils(iq).gt.0)
[38]	632	c -----------------------------------------------------------------
	633	c GENERAL CASE
	634	c -----------------------------------------------------------------
[626]	635	else
[358]	636	call newsedim(ngrid,nlay,1,1,ptimestep,
[1005]	637	& pplev,masse,epaisseur,zt,radius(iq),rho_q(iq),
[38]	638	& zqi(1,1,iq),wq,1.0)
	639	c Tendencies
	640	c ~~~~~~~~~~
	641	do ig=1,ngrid
	642	pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
	643	end do
	644	endif ! of if doubleq and if water
	645	c -----------------------------------------------------------------
	646
[358]	647	c Compute the final tendency:
	648	c ---------------------------
[38]	649	DO l = 1, nlay
	650	DO ig=1,ngrid
	651	pdqsed(ig,l,iq)=(zqi(ig,l,iq)-
	652	$ (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
[2322]	653	!MVals: Special case of isotopes: for now only HDO
	654	if (nqfils(iq).gt.0) then
[2332]	655	if (iq.eq.igcm_h2o_ice) then
	656	iq2=igcm_hdo_ice
	657	else
	658	call abort_physic("callsedim_mod","invalid isotope",1)
	659	endif
[2322]	660	pdqsed(ig,l,iq2)=(zqi(ig,l,iq2)-
	661	$ (pq(ig,l,iq2) + pdqfi(ig,l,iq2)*ptimestep))/ptimestep
	662	endif
[38]	663	ENDDO
	664	ENDDO
	665
	666	endif ! of if(radius(iq).gt.1.e-9)
	667	c =================================================================
	668	enddo ! of do iq=1,nq
[1005]	669
[358]	670	c Update the dust particle size "rdust"
	671	c -------------------------------------
[635]	672	if (doubleq) then
	673	DO l = 1, nlay
[358]	674	DO ig=1,ngrid
[740]	675
	676
	677	call updaterdust(zqi(ig,l,igcm_dust_mass),
	678	& zqi(ig,l,igcm_dust_number),rdust(ig,l),
	679	& tauscaling(ig))
	680
	681
[358]	682	ENDDO
[635]	683	ENDDO
	684	endif ! of if (doubleq)
[1974]	685
	686	if (rdstorm) then
	687	DO l = 1, nlay
	688	DO ig=1,ngrid
	689	call updaterdust(zqi(ig,l,igcm_stormdust_mass),
	690	& zqi(ig,l,igcm_stormdust_number),rstormdust(ig,l),
	691	& tauscaling(ig))
	692	ENDDO
	693	ENDDO
	694	endif ! of if (rdstorm)
[2199]	695
	696	if (slpwind) then
	697	DO l = 1, nlay
	698	DO ig=1,ngrid
	699	call updaterdust(zqi(ig,l,igcm_topdust_mass),
	700	& zqi(ig,l,igcm_topdust_number),rtopdust(ig,l),
	701	& tauscaling(ig))
	702	ENDDO
	703	ENDDO
	704	endif ! of if (slpwind)
[1974]	705
[411]	706	c Update the ice particle size "rice"
	707	c -------------------------------------
[635]	708	if (water) then
[740]	709	IF(microphys) THEN
	710
	711
[626]	712	DO l = 1, nlay
	713	DO ig=1,ngrid
[740]	714
	715	call updaterice_micro(zqi(ig,l,igcm_h2o_ice),
	716	& zqi(ig,l,igcm_ccn_mass),zqi(ig,l,igcm_ccn_number),
	717	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
[626]	718
	719	ENDDO
	720	ENDDO
[740]	721
[635]	722	ELSE
[740]	723
[626]	724	DO l = 1, nlay
	725	DO ig=1,ngrid
[740]	726
	727	call updaterice_typ(zqi(ig,l,igcm_h2o_ice),
	728	& tau(ig,1),zlay(ig,l),rice(ig,l))
	729
[626]	730	ENDDO
	731	ENDDO
[740]	732	ENDIF ! of IF(microphys)
[635]	733	endif ! of if (water)
[1962]	734	END SUBROUTINE callsedim
	735
	736	END MODULE callsedim_mod
[1974]	737

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

Download in other formats:

Original Format