Context Navigation

callsedim.F @ 1448

Last change on this file since 1448 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

Remove all reference/use of nlayermx and dimphys.h
Made use of automatic arrays whenever arrays are needed with dimension nlayer
Remove lots of obsolete reference to dimensions.h
Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

Deleted param.h and put contents into module param_v4_h
Adapted testphys1d, newstart, etc...
Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

Optimized computations in paramfoto_compact (twice less dlog10 calculations)
Checked consistency before/after modification in debug mode
Checked performance is not impacted (same as before)

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

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

Download in other formats:

Original Format