Context Navigation

callsedim.F @ 1818

Last change on this file since 1818 was 1779, checked in by aslmd, 7 years ago
LMDZ.MARS (purely comments) marked the absolute firstcalls not supposed to change with runtime (e.g. not domain-related). this is most of them. those firstcall can stay local and do not need to be linked with the caller's general firstcall.
File size: 16.7 KB

Rev	Line
[38]	1	SUBROUTINE callsedim(ngrid,nlay, ptimestep,
[1005]	2	& pplev,zlev, zlay, pt, pdt, rdust, rice,
[358]	3	& rsedcloud,rhocloud,
[411]	4	& pq, pdqfi, pdqsed,pdqs_sed,nq,
[740]	5	& tau,tauscaling)
[520]	6	! to use 'getin'
[1036]	7	USE ioipsl_getincom, only: getin
	8	USE updaterad, only: updaterdust,updaterice_micro,updaterice_typ
	9	USE tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number,
	10	& rho_dust, rho_q, radius, varian,
	11	& igcm_ccn_mass, igcm_ccn_number,
[1617]	12	& igcm_h2o_ice, nuice_sed, nuice_ref,
	13	& igcm_ccnco2_mass, igcm_ccnco2_number,
	14	& igcm_co2_ice
[1226]	15	USE comcstfi_h
[38]	16	IMPLICIT NONE
	17
	18	c=======================================================================
	19	c Sedimentation of the Martian aerosols
	20	c depending on their density and radius
	21	c
	22	c F.Forget 1999
	23	c
	24	c Modified by J.-B. Madeleine 2010: Now includes the doubleq
	25	c technique in order to have only one call to callsedim in
	26	c physiq.F.
	27	c
[1617]	28	c Modified by J. Audouard 09/16: Now includes the co2clouds case
	29	c If the co2 microphysics is on, then co2 theice & ccn tracers
	30	c are being sedimented in the microtimestep (co2cloud.F), not
	31	c in this routine.
	32	c
[38]	33	c=======================================================================
	34
	35	c-----------------------------------------------------------------------
	36	c declarations:
	37	c -------------
[520]	38
[38]	39	#include "callkeys.h"
	40
	41	c
	42	c arguments:
	43	c ----------
	44
[1005]	45	integer,intent(in) :: ngrid ! number of horizontal grid points
	46	integer,intent(in) :: nlay ! number of atmospheric layers
	47	real,intent(in) :: ptimestep ! physics time step (s)
	48	real,intent(in) :: pplev(ngrid,nlay+1) ! pressure at inter-layers (Pa)
	49	real,intent(in) :: zlev(ngrid,nlay+1) ! altitude at layer boundaries
	50	real,intent(in) :: zlay(ngrid,nlay) ! altitude at the middle of the layers
	51	real,intent(in) :: pt(ngrid,nlay) ! temperature at mid-layer (K)
	52	real,intent(in) :: pdt(ngrid,nlay) ! tendency on temperature, from
	53	! previous processes (K/s)
[38]	54	c Aerosol radius provided by the water ice microphysical scheme:
[1005]	55	real,intent(out) :: rdust(ngrid,nlay) ! Dust geometric mean radius (m)
	56	real,intent(out) :: rice(ngrid,nlay) ! H2O Ice geometric mean radius (m)
	57	c Sedimentation radius of water ice
[1047]	58	real,intent(in) :: rsedcloud(ngrid,nlay)
[1005]	59	c Cloud density (kg.m-3)
[1047]	60	real,intent(inout) :: rhocloud(ngrid,nlay)
[38]	61	c Traceurs :
[1005]	62	real,intent(in) :: pq(ngrid,nlay,nq) ! tracers (kg/kg)
	63	real,intent(in) :: pdqfi(ngrid,nlay,nq) ! tendency before sedimentation (kg/kg.s-1)
	64	real,intent(out) :: pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1)
	65	real,intent(out) :: pdqs_sed(ngrid,nq) ! flux at surface (kg.m-2.s-1)
	66	integer,intent(in) :: nq ! number of tracers
	67	real,intent(in) :: tau(ngrid,nlay) ! dust opacity
	68	real,intent(in) :: tauscaling(ngrid)
[38]	69
	70	c local:
	71	c ------
	72
	73	INTEGER l,ig, iq
	74
[1005]	75	real zqi(ngrid,nlay,nq) ! to locally store tracers
	76	real zt(ngrid,nlay) ! to locally store temperature
	77	real masse (ngrid,nlay) ! Layer mass (kg.m-2)
	78	real epaisseur (ngrid,nlay) ! Layer thickness (m)
	79	real wq(ngrid,nlay+1) ! displaced tracer mass (kg.m-2)
	80	real r0(ngrid,nlay) ! geometric mean radius used for
[358]	81	! sedimentation (m)
[1005]	82	real r0dust(ngrid,nlay) ! geometric mean radius used for
[358]	83	! dust (m)
[1005]	84	! real r0ccn(ngrid,nlay) ! geometric mean radius used for
[740]	85	! ! CCNs (m)
[1005]	86	real,save :: beta ! correction for the shape of the ice particles (cf. newsedim)
	87
[411]	88	c for ice radius computation
	89	REAL Mo,No
	90	REAl ccntyp
[38]	91
[358]	92
[411]	93
[38]	94	c Discrete size distributions (doubleq)
	95	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	96	c 1) Parameters used to represent the changes in fall
	97	c velocity as a function of particle size;
[1005]	98	integer ir
	99	integer,parameter :: nr=12 !(nr=7) ! number of bins
	100	real,save :: rd(nr)
	101	real qr(ngrid,nlay,nr)
	102	real,save :: rdi(nr+1) ! extreme and intermediate radii
	103	real Sq(ngrid,nlay)
	104	real,parameter :: rdmin=1.e-8
	105	real,parameter :: rdmax=30.e-6
	106	real,parameter :: rdimin=1.e-8 ! 1.e-7
	107	real,parameter :: rdimax=1.e-4
[38]	108
	109	c 2) Second size distribution for the log-normal integration
	110	c (the mass mixing ratio is computed for each radius)
	111
[1005]	112	integer iint
	113	integer,parameter :: ninter=4 ! number of points between each rdi radii
	114	real,save :: rr(ninter,nr)
[38]	115	integer radpower
[358]	116	real sigma0
[38]	117
	118	c 3) Other local variables used in doubleq
	119
[1005]	120	INTEGER,SAVE :: idust_mass ! index of tracer containing dust mass
	121	! mix. ratio
	122	INTEGER,SAVE :: idust_number ! index of tracer containing dust number
	123	! mix. ratio
	124	INTEGER,SAVE :: iccn_mass ! index of tracer containing CCN mass
	125	! mix. ratio
	126	INTEGER,SAVE :: iccn_number ! index of tracer containing CCN number
	127	! mix. ratio
[1617]	128	INTEGER,SAVE :: iccnco2_number ! index of tracer containing CCN number
	129	INTEGER,SAVE :: iccnco2_mass ! index of tracer containing CCN number
	130	INTEGER,SAVE :: ico2_ice ! index of tracer containing CCN number
[38]	131
[1617]	132
[1005]	133	LOGICAL,SAVE :: firstcall=.true.
[38]	134
	135	c ** un petit test de coherence
	136	c --------------------------
	137
[1779]	138	! AS: firstcall OK absolute
[38]	139	IF (firstcall) THEN
[411]	140
[38]	141	c Doubleq: initialization
	142	IF (doubleq) THEN
	143	do ir=1,nr
	144	rd(ir)= rdmin(rdmax/rdmin)*(float(ir-1)/float(nr-1))
	145	end do
	146	rdi(1)=rdimin
	147	do ir=2,nr
	148	rdi(ir)= sqrt(rd(ir-1)*rd(ir))
	149	end do
	150	rdi(nr+1)=rdimax
	151
	152	do ir=1,nr
	153	do iint=1,ninter
	154	rr(iint,ir)=
	155	& rdi(ir)*
	156	& (rdi(ir+1)/rdi(ir))**(float(iint-1)/float(ninter-1))
	157	c write(,) rr(iint,ir)
	158	end do
	159	end do
	160
	161	! identify tracers corresponding to mass mixing ratio and
	162	! number mixing ratio
	163	idust_mass=0 ! dummy initialization
	164	idust_number=0 ! dummy initialization
	165
	166	do iq=1,nq
	167	if (noms(iq).eq."dust_mass") then
	168	idust_mass=iq
[1005]	169	write(,)"callsedim: idust_mass=",idust_mass
[38]	170	endif
	171	if (noms(iq).eq."dust_number") then
	172	idust_number=iq
[1005]	173	write(,)"callsedim: idust_number=",idust_number
[38]	174	endif
	175	enddo
	176
	177	! check that we did find the tracers
	178	if ((idust_mass.eq.0).or.(idust_number.eq.0)) then
	179	write(,) 'callsedim: error! could not identify'
	180	write(,) ' tracers for dust mass and number mixing'
	181	write(,) ' ratio and doubleq is activated!'
	182	stop
	183	endif
	184	ENDIF !of if (doubleq)
	185
[740]	186	IF (microphys) THEN
[358]	187	iccn_mass=0
	188	iccn_number=0
	189	do iq=1,nq
	190	if (noms(iq).eq."ccn_mass") then
	191	iccn_mass=iq
[1005]	192	write(,)"callsedim: iccn_mass=",iccn_mass
[358]	193	endif
	194	if (noms(iq).eq."ccn_number") then
	195	iccn_number=iq
[1005]	196	write(,)"callsedim: iccn_number=",iccn_number
[358]	197	endif
	198	enddo
	199	! check that we did find the tracers
	200	if ((iccn_mass.eq.0).or.(iccn_number.eq.0)) then
	201	write(,) 'callsedim: error! could not identify'
	202	write(,) ' tracers for ccn mass and number mixing'
[740]	203	write(,) ' ratio and microphys is activated!'
[358]	204	stop
	205	endif
[740]	206	ENDIF !of if (microphys)
[358]	207
[1720]	208	IF (co2clouds) THEN
[1617]	209	iccnco2_mass=0
	210	iccnco2_number=0
	211	ico2_ice=0
	212	do iq=1,nq
	213	if (noms(iq).eq."ccnco2_mass") then
	214	iccnco2_mass=iq
	215	write(,)"callsedim: iccnco2_mass=",iccnco2_mass
	216	endif
	217	if (noms(iq).eq."co2_ice") then
	218	ico2_ice=iq
	219	write(,)"callsedim: ico2_ice=",ico2_ice
	220	endif
	221	if (noms(iq).eq."ccnco2_number") then
	222	iccnco2_number=iq
	223	write(,)"callsedim: iccnco2_number=",iccnco2_number
	224	endif
	225	enddo
	226	! check that we did find the tracers
	227	if ((iccnco2_mass.eq.0).or.(iccnco2_number.eq.0)) then
	228	write(,) 'callsedim: error! could not identify'
	229	write(,) ' tracers for ccn co2 mass and number mixing'
[1720]	230	write(,) ' ratio and co2clouds are activated!'
[1617]	231	stop
	232	endif
[1720]	233	ENDIF !of if (co2clouds)
[1617]	234
	235
[626]	236	IF (water) THEN
	237	write(,) "correction for the shape of the ice particles ?"
	238	beta=0.75 ! default value
	239	call getin("ice_shape",beta)
	240	write(,) " ice_shape = ",beta
	241
	242	write(,) "water_param nueff Sedimentation:", nuice_sed
	243	IF (activice) THEN
	244	write(,) "water_param nueff Radiative:", nuice_ref
	245	ENDIF
	246	ENDIF
[38]	247
	248	firstcall=.false.
	249	ENDIF ! of IF (firstcall)
	250
	251	c-----------------------------------------------------------------------
	252	c 1. Initialization
	253	c -----------------
	254
[1005]	255	! zqi(1:ngrid,1:nlay,1:nqmx) = 0.
	256	c Update the mass mixing ratio and temperature with the tendencies coming
[38]	257	c from other parameterizations:
	258	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1005]	259	zqi(1:ngrid,1:nlay,1:nq)=pq(1:ngrid,1:nlay,1:nq)
	260	& +pdqfi(1:ngrid,1:nlay,1:nq)*ptimestep
	261	zt(1:ngrid,1:nlay)=pt(1:ngrid,1:nlay)
	262	& +pdt(1:ngrid,1:nlay)*ptimestep
[38]	263
	264
	265	c Computing the different layer properties
	266	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	267	c Mass (kg.m-2), thickness(m), crossing time (s) etc.
	268
	269	do l=1,nlay
	270	do ig=1, ngrid
	271	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
	272	epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
	273	end do
	274	end do
	275
	276	c =================================================================
[358]	277	c Compute the geometric mean radius used for sedimentation
	278
	279	if (doubleq) then
	280	do l=1,nlay
	281	do ig=1, ngrid
[740]	282
	283	call updaterdust(zqi(ig,l,igcm_dust_mass),
	284	& zqi(ig,l,igcm_dust_number),r0dust(ig,l),
	285	& tauscaling(ig))
	286
[358]	287	end do
	288	end do
	289	endif
	290
[740]	291
[358]	292	c =================================================================
[38]	293	do iq=1,nq
[1617]	294	if(radius(iq).gt.1.e-9 .and.(iq.ne.ico2_ice) .and.
	295	& (iq .ne. iccnco2_mass) .and. (iq .ne.
	296	& iccnco2_number)) then ! no sedim for gaz or CO2 clouds (done in microtimestep)
[38]	297
	298	c -----------------------------------------------------------------
	299	c DOUBLEQ CASE
	300	c -----------------------------------------------------------------
[358]	301	if ((doubleq.and.
[38]	302	& ((iq.eq.idust_mass).or.
[740]	303	& (iq.eq.idust_number)))) then
	304
[38]	305	c Computing size distribution:
	306	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	307
[520]	308	c if ((iq.eq.idust_mass).or.(iq.eq.idust_number)) then
[358]	309	do l=1,nlay
	310	do ig=1, ngrid
	311	r0(ig,l)=r0dust(ig,l)
	312	end do
[38]	313	end do
[358]	314	sigma0 = varian
[38]	315
	316	c Computing mass mixing ratio for each particle size
	317	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[358]	318	IF ((iq.EQ.idust_mass).or.(iq.EQ.iccn_mass)) then
[38]	319	radpower = 2
[520]	320	ELSE ! number
[38]	321	radpower = -1
	322	ENDIF
	323	Sq(1:ngrid,1:nlay) = 0.
	324	do ir=1,nr
	325	do l=1,nlay
	326	do ig=1,ngrid
	327	c ****************
	328	c Size distribution integration
	329	c (Trapezoid Integration Method)
	330	qr(ig,l,ir)=0.5(rr(2,ir)-rr(1,ir))
	331	& (rr(1,ir)*radpower)
[358]	332	& exp(-(log(rr(1,ir)/r0(ig,l)))*2/(2sigma0**2))
[38]	333	do iint=2,ninter-1
	334	qr(ig,l,ir)=qr(ig,l,ir) +
	335	& 0.5(rr(iint+1,ir)-rr(iint-1,ir))
	336	& (rr(iint,ir)*radpower)
	337	& exp(-(log(rr(iint,ir)/r0(ig,l)))**2/
[358]	338	& (2sigma0*2))
[38]	339	end do
	340	qr(ig,l,ir)=qr(ig,l,ir) +
	341	& 0.5(rr(ninter,ir)-rr(ninter-1,ir))
	342	& (rr(ninter,ir)*radpower)
	343	& exp(-(log(rr(ninter,ir)/r0(ig,l)))**2/
[358]	344	& (2sigma0*2))
[38]	345
	346	c **************** old method (not recommended!)
	347	c qr(ig,l,ir)=(rd(ir)*(5-3iq))*
[358]	348	c & exp( -(log(rd(ir)/r0(ig,l)))*2 / (2sigma0**2) )
[38]	349	c ******************************
	350
	351	Sq(ig,l)=Sq(ig,l)+qr(ig,l,ir)
	352	enddo
	353	enddo
	354	enddo
	355
	356	do ir=1,nr
	357	do l=1,nlay
	358	do ig=1,ngrid
	359	qr(ig,l,ir) = zqi(ig,l,iq)*qr(ig,l,ir)/Sq(ig,l)
	360	enddo
	361	enddo
	362	enddo
	363
	364	c Computing sedimentation for each tracer
	365	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	366
	367	zqi(1:ngrid,1:nlay,iq) = 0.
	368	pdqs_sed(1:ngrid,iq) = 0.
	369
	370	do ir=1,nr
[740]	371
[358]	372	call newsedim(ngrid,nlay,1,1,ptimestep,
[1005]	373	& pplev,masse,epaisseur,zt,rd(ir),rho_dust,qr(1,1,ir),
[358]	374	& wq,0.5)
[38]	375
	376	c Tendencies
	377	c ~~~~~~~~~~
	378	do ig=1,ngrid
	379	pdqs_sed(ig,iq) = pdqs_sed(ig,iq)
	380	& + wq(ig,1)/ptimestep
	381	end do
	382	DO l = 1, nlay
	383	DO ig=1,ngrid
	384	zqi(ig,l,iq)=zqi(ig,l,iq)+qr(ig,l,ir)
	385	ENDDO
	386	ENDDO
	387	enddo ! of do ir=1,nr
	388	c -----------------------------------------------------------------
[626]	389	c WATER CYCLE CASE
[38]	390	c -----------------------------------------------------------------
[626]	391	else if ((iq .eq. iccn_mass) .or. (iq .eq. iccn_number)
	392	& .or. (iq .eq. igcm_h2o_ice)) then
	393	if (microphys) then
	394	! water ice sedimentation
	395	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
[1005]	396	& ptimestep,pplev,masse,epaisseur,zt,rsedcloud,rhocloud,
[626]	397	& zqi(1,1,iq),wq,beta)
	398	else
	399	! water ice sedimentation
	400	call newsedim(ngrid,nlay,ngrid*nlay,1,
[1005]	401	& ptimestep,pplev,masse,epaisseur,zt,rsedcloud,rho_q(iq),
[626]	402	& zqi(1,1,iq),wq,beta)
	403	endif ! of if (microphys)
[38]	404	c Tendencies
	405	c ~~~~~~~~~~
[626]	406	do ig=1,ngrid
	407	pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
	408	end do
[38]	409	c -----------------------------------------------------------------
	410	c GENERAL CASE
	411	c -----------------------------------------------------------------
[626]	412	else
[358]	413	call newsedim(ngrid,nlay,1,1,ptimestep,
[1005]	414	& pplev,masse,epaisseur,zt,radius(iq),rho_q(iq),
[38]	415	& zqi(1,1,iq),wq,1.0)
	416	c Tendencies
	417	c ~~~~~~~~~~
	418	do ig=1,ngrid
	419	pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
	420	end do
	421	endif ! of if doubleq and if water
	422	c -----------------------------------------------------------------
	423
[358]	424	c Compute the final tendency:
	425	c ---------------------------
[38]	426	DO l = 1, nlay
	427	DO ig=1,ngrid
	428	pdqsed(ig,l,iq)=(zqi(ig,l,iq)-
	429	$ (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
	430	ENDDO
	431	ENDDO
	432
	433	endif ! of if(radius(iq).gt.1.e-9)
	434	c =================================================================
	435	enddo ! of do iq=1,nq
[1005]	436
[358]	437	c Update the dust particle size "rdust"
	438	c -------------------------------------
[635]	439	if (doubleq) then
	440	DO l = 1, nlay
[358]	441	DO ig=1,ngrid
[740]	442
	443
	444	call updaterdust(zqi(ig,l,igcm_dust_mass),
	445	& zqi(ig,l,igcm_dust_number),rdust(ig,l),
	446	& tauscaling(ig))
	447
	448
[358]	449	ENDDO
[635]	450	ENDDO
	451	endif ! of if (doubleq)
[411]	452
	453	c Update the ice particle size "rice"
	454	c -------------------------------------
[635]	455	if (water) then
[740]	456	IF(microphys) THEN
	457
	458
[626]	459	DO l = 1, nlay
	460	DO ig=1,ngrid
[740]	461
	462	call updaterice_micro(zqi(ig,l,igcm_h2o_ice),
	463	& zqi(ig,l,igcm_ccn_mass),zqi(ig,l,igcm_ccn_number),
	464	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
[626]	465
	466	ENDDO
	467	ENDDO
[740]	468
[635]	469	ELSE
[740]	470
[626]	471	DO l = 1, nlay
	472	DO ig=1,ngrid
[740]	473
	474	call updaterice_typ(zqi(ig,l,igcm_h2o_ice),
	475	& tau(ig,1),zlay(ig,l),rice(ig,l))
	476
[626]	477	ENDDO
	478	ENDDO
[740]	479	ENDIF ! of IF(microphys)
[635]	480	endif ! of if (water)
[358]	481
[38]	482	END
	483

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

Download in other formats:

Original Format