Context Navigation

callsedim_mod.F @ 2323

Last change on this file since 2323 was 2323, checked in by mvals, 5 years ago
Mars GCM: Follow-up of the last commit for the transport of the isotopic ratio: forgot the initialisation of zq0() in callsedim_mod.F ! MV
File size: 25.0 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,
	25	& iqfils,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
[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
	160	INTEGER,SAVE :: ico2_ice ! index of tracer containing CCN number
[38]	161
[1617]	162
[1005]	163	LOGICAL,SAVE :: firstcall=.true.
[38]	164
[1974]	165
	166
[38]	167	c ** un petit test de coherence
	168	c --------------------------
[1779]	169	! AS: firstcall OK absolute
[38]	170	IF (firstcall) THEN
[411]	171
[38]	172	c Doubleq: initialization
	173	IF (doubleq) THEN
	174	do ir=1,nr
	175	rd(ir)= rdmin(rdmax/rdmin)*(float(ir-1)/float(nr-1))
	176	end do
	177	rdi(1)=rdimin
	178	do ir=2,nr
	179	rdi(ir)= sqrt(rd(ir-1)*rd(ir))
	180	end do
	181	rdi(nr+1)=rdimax
	182
	183	do ir=1,nr
	184	do iint=1,ninter
	185	rr(iint,ir)=
	186	& rdi(ir)*
	187	& (rdi(ir+1)/rdi(ir))**(float(iint-1)/float(ninter-1))
	188	c write(,) rr(iint,ir)
	189	end do
	190	end do
	191
	192	! identify tracers corresponding to mass mixing ratio and
	193	! number mixing ratio
	194	idust_mass=0 ! dummy initialization
	195	idust_number=0 ! dummy initialization
	196
	197	do iq=1,nq
	198	if (noms(iq).eq."dust_mass") then
	199	idust_mass=iq
[1005]	200	write(,)"callsedim: idust_mass=",idust_mass
[38]	201	endif
	202	if (noms(iq).eq."dust_number") then
	203	idust_number=iq
[1005]	204	write(,)"callsedim: idust_number=",idust_number
[38]	205	endif
	206	enddo
	207
	208	! check that we did find the tracers
	209	if ((idust_mass.eq.0).or.(idust_number.eq.0)) then
	210	write(,) 'callsedim: error! could not identify'
	211	write(,) ' tracers for dust mass and number mixing'
	212	write(,) ' ratio and doubleq is activated!'
[2304]	213	call abort_physic(modname,"missing dust tracers",1)
[38]	214	endif
	215	ENDIF !of if (doubleq)
	216
[740]	217	IF (microphys) THEN
[358]	218	iccn_mass=0
	219	iccn_number=0
	220	do iq=1,nq
	221	if (noms(iq).eq."ccn_mass") then
	222	iccn_mass=iq
[1005]	223	write(,)"callsedim: iccn_mass=",iccn_mass
[358]	224	endif
	225	if (noms(iq).eq."ccn_number") then
	226	iccn_number=iq
[1005]	227	write(,)"callsedim: iccn_number=",iccn_number
[358]	228	endif
	229	enddo
	230	! check that we did find the tracers
	231	if ((iccn_mass.eq.0).or.(iccn_number.eq.0)) then
	232	write(,) 'callsedim: error! could not identify'
	233	write(,) ' tracers for ccn mass and number mixing'
[740]	234	write(,) ' ratio and microphys is activated!'
[2304]	235	call abort_physic(modname,"missing ccn tracers",1)
[358]	236	endif
[740]	237	ENDIF !of if (microphys)
[358]	238
[1720]	239	IF (co2clouds) THEN
[1617]	240	iccnco2_mass=0
	241	iccnco2_number=0
	242	ico2_ice=0
	243	do iq=1,nq
	244	if (noms(iq).eq."ccnco2_mass") then
	245	iccnco2_mass=iq
	246	write(,)"callsedim: iccnco2_mass=",iccnco2_mass
	247	endif
	248	if (noms(iq).eq."co2_ice") then
	249	ico2_ice=iq
	250	write(,)"callsedim: ico2_ice=",ico2_ice
	251	endif
	252	if (noms(iq).eq."ccnco2_number") then
	253	iccnco2_number=iq
	254	write(,)"callsedim: iccnco2_number=",iccnco2_number
	255	endif
	256	enddo
	257	! check that we did find the tracers
	258	if ((iccnco2_mass.eq.0).or.(iccnco2_number.eq.0)) then
	259	write(,) 'callsedim: error! could not identify'
	260	write(,) ' tracers for ccn co2 mass and number mixing'
[1720]	261	write(,) ' ratio and co2clouds are activated!'
[2304]	262	call abort_physic(modname,"missing co2 ccn tracers",1)
[1617]	263	endif
[1720]	264	ENDIF !of if (co2clouds)
[1617]	265
[1974]	266	IF (water) THEN
[626]	267	write(,) "correction for the shape of the ice particles ?"
	268	beta=0.75 ! default value
[2304]	269	call getin_p("ice_shape",beta)
[626]	270	write(,) " ice_shape = ",beta
	271
	272	write(,) "water_param nueff Sedimentation:", nuice_sed
	273	IF (activice) THEN
	274	write(,) "water_param nueff Radiative:", nuice_ref
	275	ENDIF
[1974]	276	ENDIF
	277
	278	IF (rdstorm) THEN ! identifying stormdust tracers for sedimentation
	279	istormdust_mass=0 ! dummy initialization
	280	istormdust_number=0 ! dummy initialization
	281
	282	do iq=1,nq
	283	if (noms(iq).eq."stormdust_mass") then
	284	istormdust_mass=iq
	285	write(,)"callsedim: istormdust_mass=",istormdust_mass
	286	endif
	287	if (noms(iq).eq."stormdust_number") then
	288	istormdust_number=iq
	289	write(,)"callsedim: istormdust_number=",
	290	& istormdust_number
	291	endif
	292	enddo
	293
	294	! check that we did find the tracers
	295	if ((istormdust_mass.eq.0).or.(istormdust_number.eq.0)) then
	296	write(,) 'callsedim: error! could not identify'
	297	write(,) ' tracers for stormdust mass and number mixing'
	298	write(,) ' ratio and rdstorm is activated!'
[2304]	299	call abort_physic(modname,"missing stormdust tracers",1)
[1974]	300	endif
	301	ENDIF !of if (rdstorm)
[2199]	302
	303	IF (slpwind) THEN ! identifying topdust tracers for sedimentation
	304	itopdust_mass=0 ! dummy initialization
	305	itopdust_number=0 ! dummy initialization
	306
	307	do iq=1,nq
	308	if (noms(iq).eq."topdust_mass") then
	309	itopdust_mass=iq
	310	write(,)"callsedim: itopdust_mass=",itopdust_mass
	311	endif
	312	if (noms(iq).eq."topdust_number") then
	313	itopdust_number=iq
	314	write(,)"callsedim: itopdust_number=",
	315	& itopdust_number
	316	endif
	317	enddo
	318
	319	! check that we did find the tracers
	320	if ((itopdust_mass.eq.0).or.(itopdust_number.eq.0)) then
	321	write(,) 'callsedim: error! could not identify'
	322	write(,) ' tracers for topdust mass and number mixing'
	323	write(,) ' ratio and slpwind is activated!'
[2304]	324	call abort_physic(modname,"missing topdust tracers",1)
[2199]	325	endif
	326	ENDIF !of if (slpwind)
	327
[38]	328	firstcall=.false.
	329	ENDIF ! of IF (firstcall)
	330
	331	c-----------------------------------------------------------------------
	332	c 1. Initialization
	333	c -----------------
	334
[1005]	335	! zqi(1:ngrid,1:nlay,1:nqmx) = 0.
	336	c Update the mass mixing ratio and temperature with the tendencies coming
[38]	337	c from other parameterizations:
	338	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1005]	339	zqi(1:ngrid,1:nlay,1:nq)=pq(1:ngrid,1:nlay,1:nq)
	340	& +pdqfi(1:ngrid,1:nlay,1:nq)*ptimestep
[2323]	341	zq0(1:ngrid,1:nlay,1:nq)=pq(1:ngrid,1:nlay,1:nq) !MVals: keep the input value
	342	& +pdqfi(1:ngrid,1:nlay,1:nq)*ptimestep
[1005]	343	zt(1:ngrid,1:nlay)=pt(1:ngrid,1:nlay)
	344	& +pdt(1:ngrid,1:nlay)*ptimestep
[38]	345
	346	c Computing the different layer properties
	347	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	348	c Mass (kg.m-2), thickness(m), crossing time (s) etc.
	349
	350	do l=1,nlay
	351	do ig=1, ngrid
	352	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
	353	epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
	354	end do
	355	end do
	356
	357	c =================================================================
[358]	358	c Compute the geometric mean radius used for sedimentation
	359
	360	if (doubleq) then
	361	do l=1,nlay
	362	do ig=1, ngrid
[740]	363
	364	call updaterdust(zqi(ig,l,igcm_dust_mass),
	365	& zqi(ig,l,igcm_dust_number),r0dust(ig,l),
	366	& tauscaling(ig))
	367
[358]	368	end do
	369	end do
	370	endif
[2199]	371	c rocket dust storm
[1974]	372	if (rdstorm) then
	373	do l=1,nlay
	374	do ig=1, ngrid
	375
	376	call updaterdust(zqi(ig,l,igcm_stormdust_mass),
	377	& zqi(ig,l,igcm_stormdust_number),r0stormdust(ig,l),
	378	& tauscaling(ig))
	379
	380	end do
	381	end do
	382	endif
[2199]	383	c entrainment by slope wind
	384	if (slpwind) then
	385	do l=1,nlay
	386	do ig=1, ngrid
	387
	388	call updaterdust(zqi(ig,l,igcm_topdust_mass),
	389	& zqi(ig,l,igcm_topdust_number),r0topdust(ig,l),
	390	& tauscaling(ig))
	391
	392	end do
	393	end do
	394	endif
[358]	395	c =================================================================
[38]	396	do iq=1,nq
[1617]	397	if(radius(iq).gt.1.e-9 .and.(iq.ne.ico2_ice) .and.
	398	& (iq .ne. iccnco2_mass) .and. (iq .ne.
[2322]	399	& iccnco2_number) .and. ! no sedim for gaz or CO2 clouds (done in microtimestep)
	400	& iq .ne. igcm_hdo_ice) then !MVals: hdo is transported by h2o
[38]	401	c -----------------------------------------------------------------
	402	c DOUBLEQ CASE
	403	c -----------------------------------------------------------------
[2199]	404	if ( doubleq.and.
[1974]	405	& ((iq.eq.idust_mass).or.(iq.eq.idust_number).or.
[2199]	406	& (iq.eq.istormdust_mass).or.(iq.eq.istormdust_number).or.
	407	& (iq.eq.itopdust_mass).or.(iq.eq.itopdust_number)) ) then
[740]	408
[38]	409	c Computing size distribution:
	410	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	411
[1974]	412	if ((iq.eq.idust_mass).or.(iq.eq.idust_number)) then
[358]	413	do l=1,nlay
	414	do ig=1, ngrid
	415	r0(ig,l)=r0dust(ig,l)
	416	end do
[38]	417	end do
[1974]	418	else if ((iq.eq.istormdust_mass).or.
	419	& (iq.eq.istormdust_number)) then
	420	do l=1,nlay
	421	do ig=1, ngrid
	422	r0(ig,l)=r0stormdust(ig,l)
	423	end do
	424	end do
[2199]	425	else if ((iq.eq.itopdust_mass).or.
	426	& (iq.eq.itopdust_number)) then
	427	do l=1,nlay
	428	do ig=1, ngrid
	429	r0(ig,l)=r0topdust(ig,l)
	430	end do
	431	end do
[1974]	432	endif
	433	sigma0 = varian
[38]	434
	435	c Computing mass mixing ratio for each particle size
	436	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2199]	437	IF ((iq.EQ.idust_mass).or.(iq.EQ.istormdust_mass)
	438	& .or.(iq.EQ.itopdust_mass)) then
[38]	439	radpower = 2
[520]	440	ELSE ! number
[38]	441	radpower = -1
	442	ENDIF
	443	Sq(1:ngrid,1:nlay) = 0.
	444	do ir=1,nr
	445	do l=1,nlay
	446	do ig=1,ngrid
	447	c ****************
	448	c Size distribution integration
	449	c (Trapezoid Integration Method)
	450	qr(ig,l,ir)=0.5(rr(2,ir)-rr(1,ir))
	451	& (rr(1,ir)*radpower)
[358]	452	& exp(-(log(rr(1,ir)/r0(ig,l)))*2/(2sigma0**2))
[38]	453	do iint=2,ninter-1
	454	qr(ig,l,ir)=qr(ig,l,ir) +
	455	& 0.5(rr(iint+1,ir)-rr(iint-1,ir))
	456	& (rr(iint,ir)*radpower)
	457	& exp(-(log(rr(iint,ir)/r0(ig,l)))**2/
[358]	458	& (2sigma0*2))
[38]	459	end do
	460	qr(ig,l,ir)=qr(ig,l,ir) +
	461	& 0.5(rr(ninter,ir)-rr(ninter-1,ir))
	462	& (rr(ninter,ir)*radpower)
	463	& exp(-(log(rr(ninter,ir)/r0(ig,l)))**2/
[358]	464	& (2sigma0*2))
[38]	465
	466	c **************** old method (not recommended!)
	467	c qr(ig,l,ir)=(rd(ir)*(5-3iq))*
[358]	468	c & exp( -(log(rd(ir)/r0(ig,l)))*2 / (2sigma0**2) )
[38]	469	c ******************************
	470
	471	Sq(ig,l)=Sq(ig,l)+qr(ig,l,ir)
	472	enddo
	473	enddo
	474	enddo
	475
	476	do ir=1,nr
	477	do l=1,nlay
	478	do ig=1,ngrid
	479	qr(ig,l,ir) = zqi(ig,l,iq)*qr(ig,l,ir)/Sq(ig,l)
	480	enddo
	481	enddo
	482	enddo
	483
	484	c Computing sedimentation for each tracer
	485	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	486
	487	zqi(1:ngrid,1:nlay,iq) = 0.
	488	pdqs_sed(1:ngrid,iq) = 0.
	489
	490	do ir=1,nr
[358]	491	call newsedim(ngrid,nlay,1,1,ptimestep,
[1913]	492	& pplev,masse,epaisseur,zt,rd(ir),(/rho_dust/),qr(1,1,ir),
[358]	493	& wq,0.5)
[38]	494
	495	c Tendencies
	496	c ~~~~~~~~~~
	497	do ig=1,ngrid
	498	pdqs_sed(ig,iq) = pdqs_sed(ig,iq)
	499	& + wq(ig,1)/ptimestep
	500	end do
	501	DO l = 1, nlay
	502	DO ig=1,ngrid
	503	zqi(ig,l,iq)=zqi(ig,l,iq)+qr(ig,l,ir)
	504	ENDDO
[1974]	505	ENDDO
[38]	506	enddo ! of do ir=1,nr
	507	c -----------------------------------------------------------------
[626]	508	c WATER CYCLE CASE
[38]	509	c -----------------------------------------------------------------
[626]	510	else if ((iq .eq. iccn_mass) .or. (iq .eq. iccn_number)
[2322]	511	& .or. (iq .eq. igcm_h2o_ice)) then
[626]	512	if (microphys) then
[2322]	513	c water ice sedimentation
	514	c ~~~~~~~~~~
[626]	515	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
[1005]	516	& ptimestep,pplev,masse,epaisseur,zt,rsedcloud,rhocloud,
[626]	517	& zqi(1,1,iq),wq,beta)
[2322]	518	c Surface Tendencies
	519	c ~~~~~~~~~~
	520	do ig=1,ngrid
	521	pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
	522	end do
[626]	523	else
[2322]	524	c water ice sedimentation
	525	c ~~~~~~~~~~
[626]	526	call newsedim(ngrid,nlay,ngrid*nlay,1,
[1005]	527	& ptimestep,pplev,masse,epaisseur,zt,rsedcloud,rho_q(iq),
[626]	528	& zqi(1,1,iq),wq,beta)
[2322]	529	c Surface tendencies
	530	c ~~~~~~~~~~
	531	do ig=1,ngrid
	532	pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
	533	end do
	534	c Special case of isotopes
	535	c ~~~~~~~~~~
	536	!MVals: Loop over the sons ("fils")
	537	if (nqfils(iq).gt.0) then
	538	iq2=iqfils(nqfils(iq),iq) ! for now always nqfils(iq)=1 (special case of HDO only son of H2O)
	539	!MVals: input paramters in vlz_fi for hdo
	540	do l=1,nlay
	541	do ig=1,ngrid
	542	if (zq0(ig,l,iq).gt.qperemin) then
	543	Ratio0(ig,l)=zq0(ig,l,iq2)/zq0(ig,l,iq)
	544	else
	545	Ratio0(ig,l)=0.
	546	endif
	547	Ratio(ig,l)=Ratio0(ig,l)
	548	masseq(ig,l)=max(masse(ig,l)*zq0(ig,l,iq),masseqmin)
	549	w(ig,l)=wq(ig,l) !MVals: very important: hdo is transported by h2o (see vlsplt_p.F: correction bugg 15 mai 2015)
	550	enddo !ig=1,ngrid
	551	enddo !l=1,nlay
	552	!MVals: no need to enter newsedim as the transporting mass w has been already calculated
	553	call vlz_fi(ngrid,nlay,Ratio,2.,masseq,w,wq)
	554	zqi(:,nlay,iq2)=zqi(:,nlay,iq)*Ratio0(:,nlay)
	555	do l=1,nlay-1
	556	do ig=1,ngrid
	557	newmasse=max(masseq(ig,l)+w(ig,l+1)-w(ig,l),masseqmin)
	558	Ratio(ig,l)=(Ratio0(ig,l)*masseq(ig,l)
	559	& +wq(ig,l+1)-wq(ig,l))/newmasse
	560	zqi(ig,l,iq2)=zqi(ig,l,iq)*Ratio(ig,l)
	561	enddo
	562	enddo !l=1,nlay-1
	563	!MVals: hdo surface tendency
	564	do ig=1,ngrid
	565	if (w(ig,1).gt.masseqmin) then
	566	pdqs_sed(ig,iq2)=pdqs_sed(ig,iq)*(wq(ig,1)/w(ig,1))
	567	else
	568	pdqs_sed(ig,iq2)=pdqs_sed(ig,iq)*Ratio0(ig,1)
	569	endif
	570	end do
	571	endif !(nqfils(iq).gt.0)
[626]	572	endif ! of if (microphys)
[2322]	573
[38]	574	c -----------------------------------------------------------------
	575	c GENERAL CASE
	576	c -----------------------------------------------------------------
[626]	577	else
[358]	578	call newsedim(ngrid,nlay,1,1,ptimestep,
[1005]	579	& pplev,masse,epaisseur,zt,radius(iq),rho_q(iq),
[38]	580	& zqi(1,1,iq),wq,1.0)
	581	c Tendencies
	582	c ~~~~~~~~~~
	583	do ig=1,ngrid
	584	pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
	585	end do
	586	endif ! of if doubleq and if water
	587	c -----------------------------------------------------------------
	588
[358]	589	c Compute the final tendency:
	590	c ---------------------------
[38]	591	DO l = 1, nlay
	592	DO ig=1,ngrid
	593	pdqsed(ig,l,iq)=(zqi(ig,l,iq)-
	594	$ (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
[2322]	595	!MVals: Special case of isotopes: for now only HDO
	596	if (nqfils(iq).gt.0) then
	597	iq2=iqfils(nqfils(iq),iq)
	598	pdqsed(ig,l,iq2)=(zqi(ig,l,iq2)-
	599	$ (pq(ig,l,iq2) + pdqfi(ig,l,iq2)*ptimestep))/ptimestep
	600	endif
[38]	601	ENDDO
	602	ENDDO
	603
	604	endif ! of if(radius(iq).gt.1.e-9)
	605	c =================================================================
	606	enddo ! of do iq=1,nq
[1005]	607
[358]	608	c Update the dust particle size "rdust"
	609	c -------------------------------------
[635]	610	if (doubleq) then
	611	DO l = 1, nlay
[358]	612	DO ig=1,ngrid
[740]	613
	614
	615	call updaterdust(zqi(ig,l,igcm_dust_mass),
	616	& zqi(ig,l,igcm_dust_number),rdust(ig,l),
	617	& tauscaling(ig))
	618
	619
[358]	620	ENDDO
[635]	621	ENDDO
	622	endif ! of if (doubleq)
[1974]	623
	624	if (rdstorm) then
	625	DO l = 1, nlay
	626	DO ig=1,ngrid
	627	call updaterdust(zqi(ig,l,igcm_stormdust_mass),
	628	& zqi(ig,l,igcm_stormdust_number),rstormdust(ig,l),
	629	& tauscaling(ig))
	630	ENDDO
	631	ENDDO
	632	endif ! of if (rdstorm)
[2199]	633
	634	if (slpwind) then
	635	DO l = 1, nlay
	636	DO ig=1,ngrid
	637	call updaterdust(zqi(ig,l,igcm_topdust_mass),
	638	& zqi(ig,l,igcm_topdust_number),rtopdust(ig,l),
	639	& tauscaling(ig))
	640	ENDDO
	641	ENDDO
	642	endif ! of if (slpwind)
[1974]	643
[411]	644	c Update the ice particle size "rice"
	645	c -------------------------------------
[635]	646	if (water) then
[740]	647	IF(microphys) THEN
	648
	649
[626]	650	DO l = 1, nlay
	651	DO ig=1,ngrid
[740]	652
	653	call updaterice_micro(zqi(ig,l,igcm_h2o_ice),
	654	& zqi(ig,l,igcm_ccn_mass),zqi(ig,l,igcm_ccn_number),
	655	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
[626]	656
	657	ENDDO
	658	ENDDO
[740]	659
[635]	660	ELSE
[740]	661
[626]	662	DO l = 1, nlay
	663	DO ig=1,ngrid
[740]	664
	665	call updaterice_typ(zqi(ig,l,igcm_h2o_ice),
	666	& tau(ig,1),zlay(ig,l),rice(ig,l))
	667
[626]	668	ENDDO
	669	ENDDO
[740]	670	ENDIF ! of IF(microphys)
[635]	671	endif ! of if (water)
[358]	672
[1962]	673	END SUBROUTINE callsedim
	674
	675	END MODULE callsedim_mod
[1974]	676

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

Download in other formats:

Original Format