Context Navigation

callsedim_mod.F @ 3026

Last change on this file since 3026 was 2628, checked in by abierjon, 3 years ago

Mars GCM:
Big changes on mountain top dust flows for GCM6:

the scheme now activates only in grid meshes that contain a mountain among a hard-written list, instead of every meshes. This is done to prevent strong artificial reinjections of dust in places that don't present huge converging slopes enabling the concentration of dust (ex: Valles Marineris, Hellas). Topdust is now also detrained as soon as it leaves the column it originated from.
the list of the 19 allowed mountains is used by the subroutine topmons_setup in module topmons_mod, to compute a logical array contains_mons(ngrid). alpha_hmons and hsummit are also set up once and for all by this subroutine, which is called in physiq_mod's firstcall.
contains_mons, alpha_hmons and hsummit are now saved variables of the module surfdat_h, and are called as such and not as arguments in the subroutines using them.
the logical flag "slpwind" and the comments in the code have also been updated to the new terminology "mountain top dust flows", accordingly to ticket #71. The new flag read in callphys.def is "topflows".

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

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

Download in other formats:

Original Format