Context Navigation

callsedim_mod.F @ 2429

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

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

Download in other formats:

Original Format