Context Navigation

callsedim.F @ 1009

Last change on this file since 1009 was 1005, checked in by emillour, 11 years ago

Mars GCM:

Bug fix: when running with photochemistry, ccns did not sediment! Fixed initracer.F. Also added that callsedim/newsedim use updated temperatures.
Converted run0 and run_mcd scripts to bash.

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

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

Download in other formats:

Original Format