Context Navigation

aeropacity.F @ 73

Last change on this file since 73 was 38, checked in by emillour, 15 years ago
Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs. EM
File size: 19.3 KB

Rev	Line
[38]	1	SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
	2	& pq,ccn,tauref,tau,aerosol,reffrad,nueffrad,
	3	& QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d)
	4
	5	! to use 'getin'
	6	USE ioipsl_getincom
	7	IMPLICIT NONE
	8	c=======================================================================
	9	c subject:
	10	c --------
	11	c Computing aerosol optical depth in each gridbox.
	12	c
	13	c author: F.Forget
	14	c ------
	15	c update F. Montmessin (water ice scheme)
	16	c and S. Lebonnois (12/06/2003) compatibility dust/ice/chemistry
	17	c update J.-B. Madeleine 2008-2009:
	18	c - added 3D scattering by aerosols;
	19	c - dustopacity transferred from physiq.F to callradite.F,
	20	c and renamed into aeropacity.F;
	21	c
	22	c input:
	23	c -----
	24	c ngrid Number of gridpoint of horizontal grid
	25	c nlayer Number of layer
	26	c nq Number of tracer
	27	c zday Date (time since Ls=0, in martian days)
	28	c ls Solar longitude (Ls) , radian
	29	c pplay,pplev pressure (Pa) in the middle and boundary of each layer
	30	c pq Dust mixing ratio (used if tracer =T and active=T).
	31	c reffrad(ngrid,nlayer,naerkind) Aerosol effective radius
	32	c QREFvis3d(ngridmx,nlayermx,naerkind) \ 3d extinction coefficients
	33	c QREFir3d(ngridmx,nlayermx,naerkind) / at reference wavelengths;
	34	c omegaREFvis3d(ngridmx,nlayermx,naerkind) \ 3d single scat. albedo
	35	c omegaREFir3d(ngridmx,nlayermx,naerkind) / at reference wavelengths;
	36	c
	37	c output:
	38	c -------
	39	c tauref Prescribed mean column optical depth at 700 Pa
	40	c tau Column total visible dust optical depth at each point
	41	c aerosol aerosol(ig,l,1) is the dust optical
	42	c depth in layer l, grid point ig
	43
	44	c
	45	c=======================================================================
	46	#include "dimensions.h"
	47	#include "dimphys.h"
	48	#include "callkeys.h"
	49	#include "comcstfi.h"
	50	#include "comgeomfi.h"
	51	#include "dimradmars.h"
	52	#include "yomaer.h"
	53	#include "tracer.h"
	54	#include "planete.h"
	55	#include "aerkind.h"
	56
	57	c-----------------------------------------------------------------------
	58	c
	59	c Declarations :
	60	c --------------
	61	c
	62	c Input/Output
	63	c ------------
	64	INTEGER ngrid,nlayer,nq
	65
	66	REAL ls,zday,expfactor
	67	REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
	68	REAL pq(ngrid,nlayer,nq)
	69	REAL tauref(ngrid), tau(ngrid,naerkind)
	70	REAL aerosol(ngrid,nlayer,naerkind)
	71	REAL dsodust(ngridmx,nlayermx)
	72	REAL reffrad(ngrid,nlayer,naerkind)
	73	REAL nueffrad(ngrid,nlayer,naerkind)
	74	REAL QREFvis3d(ngridmx,nlayermx,naerkind)
	75	REAL QREFir3d(ngridmx,nlayermx,naerkind)
	76	REAL omegaREFvis3d(ngridmx,nlayermx,naerkind)
	77	REAL omegaREFir3d(ngridmx,nlayermx,naerkind)
	78	c
	79	c Local variables :
	80	c -----------------
	81	INTEGER l,ig,iq,i,j
	82	INTEGER iaer ! Aerosol index
	83	real topdust(ngridmx)
	84	real zlsconst, zp
	85	real taueq,tauS,tauN
	86	c Mean Qext(vis)/Qext(ir) profile
	87	real msolsir(nlayermx,naerkind)
	88	c Mean Qext(ir)/Qabs(ir) profile
	89	real mqextsqabs(nlayermx,naerkind)
	90	c Variables used when multiple particle sizes are used
	91	c for dust or water ice particles in the radiative transfer
	92	c (see callradite.F for more information).
	93	REAL taudusttmp(ngridmx)! Temporary dust opacity
	94	! used before scaling
	95	REAL taudustvis(ngridmx) ! Dust opacity after scaling
	96	REAL taudusttes(ngridmx) ! Dust opacity at IR ref. wav. as
	97	! "seen" by the GCM.
	98	REAL taucloudvis(ngridmx)! Cloud opacity at visible
	99	! reference wavelength
	100	REAL taucloudtes(ngridmx)! Cloud opacity at infrared
	101	! reference wavelength using
	102	! Qabs instead of Qext
	103	! (direct comparison with TES)
	104	REAL qdust(ngridmx,nlayermx) ! True dust mass mixing ratio
	105	REAL ccn(ngridmx,nlayermx) ! Cloud condensation nuclei
	106	! (particules kg-1)
	107	REAL qtot(ngridmx) ! Dust column (kg m-2)
	108	c
	109	c local saved variables
	110	c ---------------------
	111
	112	REAL topdust0(ngridmx)
	113	SAVE topdust0
	114	c Level under which the dust mixing ratio is held constant
	115	c when computing the dust opacity in each layer
	116	c (this applies when doubleq and active are true)
	117	INTEGER, PARAMETER :: cstdustlevel = 7
	118
	119	LOGICAL firstcall
	120	DATA firstcall/.true./
	121	SAVE firstcall
	122
	123	! indexes of water ice and dust tracers:
	124	INTEGER,SAVE :: nqdust(nqmx) ! to store the indexes of dust tracers
	125	INTEGER,SAVE :: i_ice=0 ! water ice
	126	CHARACTER(LEN=20) :: txt ! to temporarly store text
	127	CHARACTER(LEN=1) :: txt2 ! to temporarly store text
	128	! indexes of dust scatterers:
	129	INTEGER,SAVE :: iaerdust(naerkind)
	130	INTEGER,SAVE :: naerdust ! number of dust scatterers
	131
	132	tau(1:ngrid,1:naerkind)=0
	133
	134	! identify tracers
	135
	136	IF (firstcall) THEN
	137	! identify scatterers that are dust
	138	naerdust=0
	139	DO iaer=1,naerkind
	140	txt=name_iaer(iaer)
	141	IF (txt(1:4).eq."dust") THEN
	142	naerdust=naerdust+1
	143	iaerdust(naerdust)=iaer
	144	ENDIF
	145	ENDDO
	146	! identify tracers which are dust
	147	i=0
	148	DO iq=1,nq
	149	txt=noms(iq)
	150	IF (txt(1:4).eq."dust") THEN
	151	i=i+1
	152	nqdust(i)=iq
	153	ENDIF
	154	ENDDO
	155
	156	IF (water.AND.activice) THEN
	157	i_ice=igcm_h2o_ice
	158	write(,) "aeropacity: i_ice=",i_ice
	159	ENDIF
	160
	161	c altitude of the top of the aerosol layer (km) at Ls=2.76rad:
	162	c in the Viking year scenario
	163	DO ig=1,ngrid
	164	topdust0(ig)=60. -22.SIN(lati(ig))*2
	165	END DO
	166
	167	c typical profile of solsir and (1-w)^(-1):
	168	msolsir(1:nlayer,1:naerkind)=0
	169	mqextsqabs(1:nlayer,1:naerkind)=0
	170	WRITE(,) "Typical profiles of solsir and Qext/Qabs(IR):"
	171	DO iaer = 1, naerkind ! Loop on aerosol kind
	172	WRITE(,) "Aerosol # ",iaer
	173	DO l=1,nlayer
	174	DO ig=1,ngridmx
	175	msolsir(l,iaer)=msolsir(l,iaer)+
	176	& QREFvis3d(ig,l,iaer)/
	177	& QREFir3d(ig,l,iaer)
	178	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)+
	179	& (1.E0-omegaREFir3d(ig,l,iaer))**(-1)
	180	ENDDO
	181	msolsir(l,iaer)=msolsir(l,iaer)/REAL(ngridmx)
	182	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)/REAL(ngridmx)
	183	ENDDO
	184	WRITE(,) "solsir: ",msolsir(:,iaer)
	185	WRITE(,) "Qext/Qabs(IR): ",mqextsqabs(:,iaer)
	186	ENDDO
	187
	188	! load value of tauvis from callphys.def (if given there,
	189	! otherwise default value read from starfi.nc file will be used)
	190	call getin("tauvis",tauvis)
	191
	192	firstcall=.false.
	193
	194	END IF
	195
	196	c Vertical column optical depth at 700.Pa
	197	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	198	IF(iaervar.eq.1) THEN
	199	do ig=1, ngridmx
	200	tauref(ig)=max(tauvis,1.e-9) ! tauvis=cste (set in callphys.def
	201	! or read in starfi
	202	end do
	203	ELSE IF (iaervar.eq.2) THEN ! << "Viking" Scenario>>
	204
	205	tauref(1) = 0.7+.3cos(ls+80.pi/180.) ! like seen by VL1
	206	do ig=2,ngrid
	207	tauref(ig) = tauref(1)
	208	end do
	209
	210	ELSE IF (iaervar.eq.3) THEN ! << "MGS" scenario >>
	211
	212	taueq= 0.2 +(0.5-0.2) (cos(0.5(ls-4.363)))**14
	213	tauS= 0.1 +(0.5-0.1) (cos(0.5(ls-4.363)))**14
	214	tauN = 0.1
	215	c if (peri_day.eq.150) then
	216	c tauS=0.1
	217	c tauN=0.1 +(0.5-0.1) (cos(0.5(ls+pi-4.363)))**14
	218	c taueq= 0.2 +(0.5-0.2) (cos(0.5(ls+pi-4.363)))**14
	219	c endif
	220	do ig=1,ngrid/2 ! Northern hemisphere
	221	tauref(ig)= tauN +
	222	& (taueq-tauN)0.5(1+tanh((45-lati(ig)180./pi)6/60))
	223	end do
	224	do ig=ngrid/2+1, ngridmx ! Southern hemisphere
	225	tauref(ig)= tauS +
	226	& (taueq-tauS)0.5(1+tanh((45+lati(ig)180./pi)6/60))
	227	end do
	228	ELSE IF ((iaervar.eq.4).or.
	229	& ((iaervar.ge.24).and.(iaervar.le.26)))
	230	& THEN ! << "TES assimilated dust scenarios >>
	231	call readtesassim(ngrid,nlayer,zday,pplev,tauref)
	232
	233	ELSE IF (iaervar.eq.5) THEN ! << Escalier Scenario>>
	234	c tauref(1) = 0.2
	235	c if ((ls.ge.210.pi/180.).and.(ls.le.330.pi/180.))
	236	c & tauref(1) = 2.5
	237	tauref(1) = 2.5
	238	if ((ls.ge.30.pi/180.).and.(ls.le.150.pi/180.))
	239	& tauref(1) = .2
	240
	241	do ig=2,ngrid
	242	tauref(ig) = tauref(1)
	243	end do
	244	ELSE
	245	stop 'problem with iaervar in aeropacity.F'
	246	ENDIF
	247
	248	c -----------------------------------------------------------------
	249	c Computing the opacity in each layer
	250	c -----------------------------------------------------------------
	251
	252	DO iaer = 1, naerkind ! Loop on aerosol kind
	253	c --------------------------------------------
	254	aerkind: SELECT CASE (name_iaer(iaer))
	255	c==================================================================
	256	CASE("dust_conrath") aerkind ! Typical dust profile
	257	c==================================================================
	258
	259	c Altitude of the top of the dust layer
	260	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	261	zlsconst=SIN(ls-2.76)
	262	if (iddist.eq.1) then
	263	do ig=1,ngrid
	264	topdust(ig)=topdustref ! constant dust layer top
	265	end do
	266
	267	else if (iddist.eq.2) then ! "Viking" scenario
	268	do ig=1,ngrid
	269	topdust(ig)=topdust0(ig)+18.*zlsconst
	270	end do
	271
	272	else if(iddist.eq.3) then !"MGS" scenario
	273	do ig=1,ngrid
	274	topdust(ig)=60.+18.*zlsconst
	275	& -(32+18zlsconst)sin(lati(ig))**4
	276	& - 8zlsconst(sin(lati(ig)))**5
	277	end do
	278	endif
	279
	280	c Optical depth in each layer :
	281	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	282	if(iddist.ge.1) then
	283
	284	expfactor=0.
	285	DO l=1,nlayer
	286	DO ig=1,ngrid
	287	c Typical mixing ratio profile
	288	if(pplay(ig,l).gt.700.
	289	$ /(988.**(topdust(ig)/70.))) then
	290	zp=(700./pplay(ig,l))**(70./topdust(ig))
	291	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
	292	else
	293	expfactor=1.e-3
	294	endif
	295	c Vertical scaling function
	296	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) *
	297	& expfactor *
	298	& QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer)
	299	ENDDO
	300	ENDDO
	301
	302	else if(iddist.eq.0) then
	303	c old dust vertical distribution function (pollack90)
	304	DO l=1,nlayer
	305	DO ig=1,ngrid
	306	zp=700./pplay(ig,l)
	307	aerosol(ig,l,1)= tauref(ig)/700. *
	308	s (pplev(ig,l)-pplev(ig,l+1))
	309	s max( exp(.03(1.-max(zp,1.))) , 1.E-3 )
	310	ENDDO
	311	ENDDO
	312	end if
	313
	314	c==================================================================
	315	CASE("dust_doubleq") aerkind! Two-moment scheme for dust
	316	c (transport of mass and number mixing ratio)
	317	c==================================================================
	318
	319	DO l=1,nlayer
	320	IF (l.LE.cstdustlevel) THEN
	321	c Opacity in the first levels is held constant to
	322	c avoid unrealistic values due to constant lifting:
	323	DO ig=1,ngrid
	324	aerosol(ig,l,iaer) =
	325	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
	326	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
	327	& pq(ig,cstdustlevel,igcm_dust_mass) *
	328	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
	329	ENDDO
	330	ELSE
	331	DO ig=1,ngrid
	332	aerosol(ig,l,iaer) =
	333	& ( 0.75 * QREFvis3d(ig,l,iaer) /
	334	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
	335	& pq(ig,l,igcm_dust_mass) *
	336	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
	337	ENDDO
	338	ENDIF
	339	ENDDO
	340
	341	c==================================================================
	342	CASE("dust_submicron") aerkind ! Small dust population
	343	c==================================================================
	344
	345	DO l=1,nlayer
	346	IF (l.LE.cstdustlevel) THEN
	347	c Opacity in the first levels is held constant to
	348	c avoid unrealistic values due to constant lifting:
	349	DO ig=1,ngrid
	350	aerosol(ig,l,iaer) =
	351	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
	352	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
	353	& pq(ig,cstdustlevel,igcm_dust_submicron) *
	354	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
	355	ENDDO
	356	ELSE
	357	DO ig=1,ngrid
	358	aerosol(ig,l,iaer) =
	359	& ( 0.75 * QREFvis3d(ig,l,iaer) /
	360	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
	361	& pq(ig,l,igcm_dust_submicron) *
	362	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
	363	ENDDO
	364	ENDIF
	365	ENDDO
	366
	367	c==================================================================
	368	CASE("h2o_ice") aerkind ! Water ice crystals
	369	c==================================================================
	370
	371	c 1. Initialization
	372	aerosol(1:ngrid,1:nlayer,iaer) = 0.
	373	taucloudvis(1:ngrid) = 0.
	374	taucloudtes(1:ngrid) = 0.
	375	c 2. Opacity calculation
	376	DO ig=1, ngrid
	377	DO l=1,nlayer
	378	aerosol(ig,l,iaer) = max(1E-20,
	379	& ( 0.75 * QREFvis3d(ig,l,iaer) /
	380	& ( rho_ice * reffrad(ig,l,iaer) ) ) *
	381	& pq(ig,l,i_ice) *
	382	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
	383	& )
	384	taucloudvis(ig) = taucloudvis(ig) + aerosol(ig,l,iaer)
	385	taucloudtes(ig) = taucloudtes(ig) + aerosol(ig,l,iaer)*
	386	& QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer) *
	387	& ( 1.E0 - omegaREFir3d(ig,l,iaer) )
	388	ENDDO
	389	ENDDO
	390	c 3. Outputs
	391	IF (ngrid.NE.1) THEN
	392	CALL WRITEDIAGFI(ngridmx,'tauVIS','tauext VIS refwvl',
	393	& ' ',2,taucloudvis)
	394	CALL WRITEDIAGFI(ngridmx,'tauTES','tauabs IR refwvl',
	395	& ' ',2,taucloudtes)
	396	IF (callstats) THEN
	397	CALL wstats(ngridmx,'tauVIS','tauext VIS refwvl',
	398	& ' ',2,taucloudvis)
	399	CALL wstats(ngridmx,'tauTES','tauabs IR refwvl',
	400	& ' ',2,taucloudtes)
	401	ENDIF
	402	ELSE
	403	c CALL writeg1d(ngrid,1,taucloudtes,'tautes','NU')
	404	ENDIF
	405	c==================================================================
	406	END SELECT aerkind
	407	c -----------------------------------
	408	ENDDO ! iaer (loop on aerosol kind)
	409
	410	c -----------------------------------------------------------------
	411	c Rescaling each layer to reproduce the choosen (or assimilated)
	412	c dust extinction opacity at visible reference wavelength, which
	413	c is originally scaled to an equivalent 700Pa pressure surface.
	414	c -----------------------------------------------------------------
	415
	416	taudusttmp(1:ngrid)=0.
	417	DO iaer=1,naerdust
	418	DO l=1,nlayer
	419	DO ig=1,ngrid
	420	c Scaling factor
	421	taudusttmp(ig) = taudusttmp(ig) +
	422	& aerosol(ig,l,iaerdust(iaer))
	423	ENDDO
	424	ENDDO
	425	ENDDO
	426	DO iaer=1,naerdust
	427	DO l=1,nlayer
	428	DO ig=1,ngrid
	429	aerosol(ig,l,iaerdust(iaer)) = max(1E-20,
	430	& tauref(ig) *
	431	& pplev(ig,1) / 700.E0 *
	432	& aerosol(ig,l,iaerdust(iaer)) /
	433	& taudusttmp(ig)
	434	& )
	435	ENDDO
	436	ENDDO
	437	ENDDO
	438
	439	c -----------------------------------------------------------------
	440	c Computing the number of condensation nuclei
	441	c -----------------------------------------------------------------
	442	DO iaer = 1, naerkind ! Loop on aerosol kind
	443	c --------------------------------------------
	444	aerkind2: SELECT CASE (name_iaer(iaer))
	445	c==================================================================
	446	CASE("dust_conrath") aerkind2 ! Typical dust profile
	447	c==================================================================
	448	DO l=1,nlayer
	449	DO ig=1,ngrid
	450	ccn(ig,l) = max(aerosol(ig,l,iaer) /
	451	& pi / QREFvis3d(ig,l,iaer) *
	452	& (1.+nueffrad(ig,l,iaer))**3. /
	453	& reffrad(ig,l,iaer)*2. g /
	454	& (pplev(ig,l)-pplev(ig,l+1)),1e-30)
	455	ENDDO
	456	ENDDO
	457	c==================================================================
	458	CASE("dust_doubleq") aerkind2!Two-moment scheme for dust
	459	c (transport of mass and number mixing ratio)
	460	c==================================================================
	461	qtot(1:ngridmx) = 0.
	462	DO l=1,nlayer
	463	DO ig=1,ngrid
	464	qdust(ig,l) = pq(ig,l,igcm_dust_mass) * tauref(ig) *
	465	& pplev(ig,1) / 700.E0 / taudusttmp(ig)
	466	qtot(ig) = qtot(ig) + qdust(ig,l) *
	467	& (pplev(ig,l)-pplev(ig,l+1)) / g
	468	ccn(ig,l) = max( ( ref_r0 /
	469	& reffrad(ig,l,iaer) )*3.
	470	& r3n_q * qdust(ig,l) ,1e-30)
	471	ENDDO
	472	ENDDO
	473	c==================================================================
	474	END SELECT aerkind2
	475	c -----------------------------------
	476	ENDDO ! iaer (loop on aerosol kind)
	477
	478	c -----------------------------------------------------------------
	479	c Column integrated visible optical depth in each point
	480	c -----------------------------------------------------------------
	481	DO iaer=1,naerkind
	482	do l=1,nlayer
	483	do ig=1,ngrid
	484	tau(ig,iaer) = tau(ig,iaer) + aerosol(ig,l,iaer)
	485	end do
	486	end do
	487	ENDDO
	488	c -----------------------------------------------------------------
	489	c Density scaled opacity and column opacity output
	490	c -----------------------------------------------------------------
	491	dsodust(1:ngrid,1:nlayer) = 0.
	492	DO iaer=1,naerdust
	493	DO l=1,nlayermx
	494	DO ig=1,ngrid
	495	dsodust(ig,l) = dsodust(ig,l) +
	496	& aerosol(ig,l,iaerdust(iaer)) * g /
	497	& (pplev(ig,l) - pplev(ig,l+1))
	498	ENDDO
	499	ENDDO
	500	IF (ngrid.NE.1) THEN
	501	write(txt2,'(i1.1)') iaer
	502	call WRITEDIAGFI(ngridmx,'taudust'//txt2,
	503	& 'Dust col opacity',
	504	& ' ',2,tau(1,iaerdust(iaer)))
	505	IF (callstats) THEN
	506	CALL wstats(ngridmx,'taudust'//txt2,
	507	& 'Dust col opacity',
	508	& ' ',2,tau(1,iaerdust(iaer)))
	509	ENDIF
	510	ENDIF
	511	ENDDO
	512
	513	IF (ngrid.NE.1) THEN
	514	c CALL WRITEDIAGFI(ngridmx,'dsodust','tau*g/dp',
	515	c & 'm2.kg-1',3,dsodust)
	516	IF (callstats) THEN
	517	CALL wstats(ngridmx,'dsodust',
	518	& 'tau*g/dp',
	519	& 'm2.kg-1',3,dsodust)
	520	ENDIF
	521	c CALL WRITEDIAGFI(ngridmx,'ccn','Cond. nuclei',
	522	c & 'part kg-1',3,ccn)
	523	ELSE
	524	CALL writeg1d(ngrid,nlayer,dsodust,'dsodust','m2.kg-1')
	525	ENDIF
	526	c -----------------------------------------------------------------
	527	return
	528	end

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

Download in other formats:

Original Format