Context Navigation

callsedim.F @ 828

Last change on this file since 828 was 740, checked in by tnavarro, 13 years ago
module for ice and radius radius computation
File size: 14.4 KB

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

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

Download in other formats:

Original Format