Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

assim_aeropacity.F @ 43

Last change on this file since 43 was 42, checked in by aslmd, 15 years ago
LMD_MM_MARS: tests avec physique completement commune avec le GCM nouvelle physique
File size: 16.1 KB

Rev	Line
[42]	1	SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,pq,
	2	& tauref,tau,aerosol,reffrad,
	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
	56	c-----------------------------------------------------------------------
	57	c
	58	c Declarations :
	59	c --------------
	60	c
	61	c Input/Output
	62	c ------------
	63	INTEGER ngrid,nlayer,nq
	64
	65	REAL ls,zday,expfactor
	66	REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
	67	REAL pq(ngrid,nlayer,nq)
	68	REAL tauref(ngrid), tau(ngrid,naerkind)
	69	REAL aerosol(ngrid,nlayer,naerkind)
	70	REAL reffrad(ngrid,nlayer,naerkind)
	71	REAL QREFvis3d(ngridmx,nlayermx,naerkind)
	72	REAL QREFir3d(ngridmx,nlayermx,naerkind)
	73	REAL omegaREFvis3d(ngridmx,nlayermx,naerkind)
	74	REAL omegaREFir3d(ngridmx,nlayermx,naerkind)
	75	c
	76	c Local variables :
	77	c -----------------
	78	INTEGER l,ig,iq,i,j
	79	INTEGER iaer ! Aerosol index
	80	real topdust(ngridmx)
	81	real zlsconst, zp
	82	real taueq,tauS,tauN
	83	real r0,reff,coefsize
	84	c Mean Qext(vis)/Qext(ir) profile
	85	real msolsir(nlayermx,naerkind)
	86	c Mean Qext(ir)/Qabs(ir) profile
	87	real mqextsqabs(nlayermx,naerkind)
	88	c Variables used when multiple particle sizes are used
	89	c for dust or water ice particles in the radiative transfer
	90	c (see callradite.F for more information).
	91	REAL taudusttmp(ngridmx)! Temporary dust opacity
	92	! used before scaling
	93	REAL taudustvis(ngridmx) ! Dust opacity after scaling
	94	REAL taudusttes(ngridmx) ! Dust opacity at IR ref. wav. as
	95	! "seen" by the GCM.
	96	REAL taucloudvis(ngridmx)! Cloud opacity at visible
	97	! reference wavelength
	98	REAL taucloudtes(ngridmx)! Cloud opacity at infrared
	99	! reference wavelength using
	100	! Qabs instead of Qext
	101	! (direct comparison with TES)
	102	c
	103	c local saved variables
	104	c ---------------------
	105
	106	REAL topdust0(ngridmx)
	107	SAVE topdust0
	108
	109	LOGICAL firstcall
	110	DATA firstcall/.true./
	111	SAVE firstcall
	112
	113	! indexes of water ice and dust tracers:
	114	INTEGER,SAVE :: nqdust(nqmx) ! to store the indexes of dust tracers
	115	INTEGER,SAVE :: i_ice=0 ! water ice
	116	CHARACTER(LEN=20) :: tracername ! to temporarly store text
	117
	118	call zerophys(ngrid*naerkind,tau)
	119
	120	! identify tracers
	121
	122	IF (firstcall) THEN
	123	! identify tracers which are dust
	124	i=0
	125	DO iq=1,nq
	126	tracername=noms(iq)
	127	IF (tracername(1:4).eq."dust") THEN
	128	i=i+1
	129	nqdust(i)=iq
	130	ENDIF
	131	ENDDO
	132
	133	IF (water.AND.activice) THEN
	134	i_ice=igcm_h2o_ice
	135	write(,) "aeropacity: i_ice=",i_ice
	136	ENDIF
	137
	138	c altitude of the top of the aerosol layer (km) at Ls=2.76rad:
	139	c in the Viking year scenario
	140	DO ig=1,ngrid
	141	topdust0(ig)=60. -22.SIN(lati(ig))*2
	142	END DO
	143
	144	c typical profile of solsir and (1-w)^(-1):
	145	call zerophys(nlayer*naerkind,msolsir)
	146	call zerophys(nlayer*naerkind,mqextsqabs)
	147	WRITE(,) "Typical profiles of solsir and Qext/Qabs(IR):"
	148	DO iaer = 1, naerkind ! Loop on aerosol kind
	149	WRITE(,) "Aerosol # ",iaer
	150	DO l=1,nlayer
	151	DO ig=1,ngridmx
	152	msolsir(l,iaer)=msolsir(l,iaer)+
	153	& QREFvis3d(ig,l,iaer)/
	154	& QREFir3d(ig,l,iaer)
	155	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)+
	156	& (1.E0-omegaREFir3d(ig,l,iaer))**(-1)
	157	ENDDO
	158	msolsir(l,iaer)=msolsir(l,iaer)/REAL(ngridmx)
	159	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)/REAL(ngridmx)
	160	ENDDO
	161	WRITE(,) "solsir: ",msolsir(:,iaer)
	162	WRITE(,) "Qext/Qabs(IR): ",mqextsqabs(:,iaer)
	163	ENDDO
	164
	165	! load value of tauvis from callphys.def (if given there,
	166	! otherwise default value read from starfi.nc file will be used)
	167	call getin("tauvis",tauvis)
	168
	169	firstcall=.false.
	170
	171	END IF
	172
	173	DO iaer = 1, naerkind ! Loop on aerosol kind
	174	c --------------------------------------------
	175	aerkind: SELECT CASE (iaer)
	176	c==================================================================
	177	CASE(1) aerkind ! Dust (iaer=1)
	178	c==================================================================
	179
	180	c -------------------------------------------------------------
	181	c 1) Prescribed dust (if tracer=F or active=F)
	182	c -------------------------------------------------------------
	183	IF ((.not.tracer) .or. (.not.active)) THEN
	184
	185	c Vertical column optical depth at 700.Pa
	186	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	187	IF(iaervar.eq.1) THEN
	188	do ig=1, ngridmx
	189	tauref(ig)=max(tauvis,1.e-9) ! tauvis=cste (set in callphys.def
	190	! or read in starfi
	191	end do
	192	ELSE IF (iaervar.eq.2) THEN ! << "Viking" Scenario>>
	193
	194	tauref(1) = 0.7+.3cos(ls+80.pi/180.) ! like seen by VL1
	195	do ig=2,ngrid
	196	tauref(ig) = tauref(1)
	197	end do
	198
	199	ELSE IF (iaervar.eq.3) THEN ! << "MGS" scenario >>
	200
	201	taueq= 0.2 +(0.5-0.2) (cos(0.5(ls-4.363)))**14
	202	tauS= 0.1 +(0.5-0.1) (cos(0.5(ls-4.363)))**14
	203	tauN = 0.1
	204	c if (peri_day.eq.150) then
	205	c tauS=0.1
	206	c tauN=0.1 +(0.5-0.1) (cos(0.5(ls+pi-4.363)))**14
	207	c taueq= 0.2 +(0.5-0.2) (cos(0.5(ls+pi-4.363)))**14
	208	c endif
	209	do ig=1,ngrid/2 ! Northern hemisphere
	210	tauref(ig)= tauN +
	211	& (taueq-tauN)0.5(1+tanh((45-lati(ig)180./pi)6/60))
	212	end do
	213	do ig=ngrid/2+1, ngridmx ! Southern hemisphere
	214	tauref(ig)= tauS +
	215	& (taueq-tauS)0.5(1+tanh((45+lati(ig)180./pi)6/60))
	216	end do
	217
	218	ELSE IF ((iaervar.eq.4).or.
	219	& ((iaervar.ge.24).and.(iaervar.le.26)))
	220	& THEN ! << "TES assimilated dust scenarios >>
	221	call readtesassim(ngrid,nlayer,zday,pplev,tauref)
	222
	223	ELSE IF (iaervar.eq.5) THEN ! << Escalier Scenario>>
	224	c tauref(1) = 0.2
	225	c if ((ls.ge.210.pi/180.).and.(ls.le.330.pi/180.))
	226	c & tauref(1) = 2.5
	227	tauref(1) = 2.5
	228	if ((ls.ge.30.pi/180.).and.(ls.le.150.pi/180.))
	229	& tauref(1) = .2
	230	do ig=2,ngrid
	231	tauref(ig) = tauref(1)
	232	end do
	233
	234	ELSE IF (iaervar.gt.99) THEN ! << input netcdf file >>
	235	c*************WRF
	236	c
	237	c 2. customized dust opacity field
	238	c ex: from assimilation
	239	call meso_readtesassim(ngrid,nlayer,zday,pplev,tauref,
	240	. iaervar)
	241	c
	242	c*************WRF
	243
	244	ELSE
	245	stop 'problem with iaervar in aeropacity.F'
	246	ENDIF
	247
	248	c Altitude of the top of the dust layer
	249	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	250	zlsconst=SIN(ls-2.76)
	251	if (iddist.eq.1) then
	252	do ig=1,ngrid
	253	topdust(ig)=topdustref ! constant dust layer top
	254	end do
	255
	256	else if (iddist.eq.2) then ! "Viking" scenario
	257	do ig=1,ngrid
	258	topdust(ig)=topdust0(ig)+18.*zlsconst
	259	end do
	260
	261	else if(iddist.eq.3) then !"MGS" scenario
	262	do ig=1,ngrid
	263	topdust(ig)=60.+18.*zlsconst
	264	& -(32+18zlsconst)sin(lati(ig))**4
	265	& - 8zlsconst(sin(lati(ig)))**5
	266	end do
	267	endif
	268
	269
	270	c Optical depth in each layer :
	271	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	272	if(iddist.ge.1) then
	273	expfactor=0.
	274	CALL zerophys(ngrid,taudusttmp)
	275	DO l=1,nlayer
	276	DO ig=1,ngrid
	277	c Typical mixing ratio profile
	278	if(pplay(ig,l).gt.700.
	279	$ /(988.**(topdust(ig)/70.))) then
	280	zp=(700./pplay(ig,l))**(70./topdust(ig))
	281	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
	282	else
	283	expfactor=1.e-3
	284	endif
	285	c Vertical scaling function
	286	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) *
	287	& expfactor *
	288	& QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer)
	289	c Scaling factor
	290	taudusttmp(ig)=taudusttmp(ig)+aerosol(ig,l,iaer)
	291	ENDDO
	292	ENDDO
	293
	294	c Rescaling each layer to reproduce the choosen (or
	295	c assimilated) dust extinction opacity at visible
	296	c reference wavelength, which is originally scaled
	297	c to an equivalent 700Pa pressure surface.
	298	DO l=1,nlayer
	299	DO ig=1,ngrid
	300	aerosol(ig,l,iaer) = tauref(ig) *
	301	& pplev(ig,1) / 700.E0 *
	302	& aerosol(ig,l,iaer) / taudusttmp(ig)
	303	ENDDO
	304	ENDDO
	305
	306	CALL zerophys(ngrid,taudustvis)
	307	CALL zerophys(ngrid,taudusttes)
	308	DO l=1,nlayer
	309	DO ig=1,ngrid
	310	taudustvis(ig) = taudustvis(ig) + aerosol(ig,l,iaer)
	311	taudusttes(ig) = taudusttes(ig) + aerosol(ig,l,iaer)*
	312	& QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer)*
	313	& ( 1.E0 - omegaREFir3d(ig,l,iaer) )
	314	ENDDO
	315	ENDDO
	316
	317	c Outputs
	318	IF (ngrid.NE.1) THEN
	319	CALL WRITEDIAGFI(ngridmx,'taudustTES','dust abs IR',
	320	& ' ',2,taudusttes)
	321	CALL wstats(ngridmx,'taudustTES','dust abs IR',
	322	& ' ',2,taudusttes)
	323	ELSE
	324	CALL writeg1d(ngrid,1,taudusttes,'taudusttes','NU')
	325	ENDIF
	326
	327	c changement dans le calcul de la distribution verticale
	328	c dans le cas des scenarios de poussieres assimiles
	329	c if (iaervar.eq.4) THEN ! TES
	330	c call zerophys(ngrid*naerkind,tau)
	331	c
	332	c do l=1,nlayer
	333	c do ig=1,ngrid
	334	c tau(ig,1)=tau(ig,1)+ aerosol(ig,l,1)
	335	c end do
	336	c end do
	337	c do l=1,nlayer
	338	c do ig=1,ngrid
	339	c aerosol(ig,l,1)=aerosol(ig,l,1)*tauref(ig)/tau(ig,1)
	340	c $ *(pplev(ig,1)/700)
	341	c end do
	342	c end do
	343	c endif
	344	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	345	else if(iddist.eq.0) then
	346	c old dust vertical distribution function (pollack90)
	347	DO l=1,nlayer
	348	DO ig=1,ngrid
	349	zp=700./pplay(ig,l)
	350	aerosol(ig,l,1)= tauref(ig)/700. *
	351	s (pplev(ig,l)-pplev(ig,l+1))
	352	s max( exp(.03(1.-max(zp,1.))) , 1.E-3 )
	353	ENDDO
	354	ENDDO
	355	end if
	356
	357	c ---------------------------------------------------------------------
	358	c 2) Transported radiatively active dust (if tracer=T and active=T)
	359	c ----------------------------------------------------------------------
	360	ELSE IF ((tracer) .and. (active)) THEN
	361	c The dust opacity is computed from q
	362
	363	c a) "doubleq" technique (transport of mass and number mixing ratio)
	364	c ~~~~~~~~~~~~~~~~~~~
	365	if(doubleq) then
	366
	367	call zerophys(ngridnlayernaerkind,aerosol)
	368
	369	c Computing effective radius :
	370	do l=1,nlayer
	371	do ig=1, ngrid
	372	r0=
	373	& (r3n_qpq(ig,l,1)/max(pq(ig,l,2),0.01))*(1./3.)
	374	r0=min(max(r0,1.e-10),500.e-6)
	375	reff= ref_r0 * r0
	376	cc If reff is small, the transported dust mean Qext
	377	c is reduced from the reference dust Qext by a factor "coefsize"
	378
	379	coefsize=min(max(2.52e6*reff-0.043 ,0.) ,1.)
	380
	381	cc It is added 1.e-8 to pq to avoid low
	382
	383	aerosol(ig,l,1)=aerosol(ig,l,1)+ 1.E-8 +
	384	& ( 0.75Qext(1)coefsize/(rho_dust*reff))
	385	& * (pq(ig,l,nqdust(1)))*
	386	c only one dust bin to use with doubleq
	387	& (pplev(ig,l)-pplev(ig,l+1))/g
	388	end do
	389	end do
	390	call zerophys(ngrid,tauref)
	391
	392	c b) Size bin technique (each aerosol can contribute to opacity))
	393	c ~~~~~~~~~~~~~~~~~~
	394	else
	395	c The dust opacity is computed from q
	396	call zerophys(ngridnlayernaerkind,aerosol)
	397	do iq=1,dustbin
	398	do l=1,nlayer
	399	do ig=1,ngrid
	400	cc qextrhor(iq) is (3/4)Qext/(rhoreff)
	401	cc It is added 1.e-8 to pq to avoid low
	402	aerosol(ig,l,1)=aerosol(ig,l,1)+
	403	& qextrhor(nqdust(iq))(pq(ig,l,nqdust(iq))+1.e-8)
	404	& (pplev(ig,l)-pplev(ig,l+1))/g
	405	end do
	406	end do
	407	end do
	408	call zerophys(ngrid,tauref)
	409	end if ! (doubleq)
	410	END IF ! (dust scenario)
	411
	412
	413	c==================================================================
	414	CASE(2) aerkind ! Water ice crystals (iaer=2)
	415	c==================================================================
	416
	417	IF (water.AND.activice) THEN
	418	c 1. Initialization
	419	CALL zerophys(ngrid*nlayer,aerosol(1,1,iaer))
	420	CALL zerophys(ngrid,taucloudvis)
	421	CALL zerophys(ngrid,taucloudtes)
	422	c 2. Opacity calculation
	423	DO ig=1, ngrid
	424	DO l=1,nlayer
	425	aerosol(ig,l,iaer) =
	426	& ( 0.75 * QREFvis3d(ig,l,iaer) /
	427	& ( rho_ice * reffrad(ig,l,iaer) ) ) *
	428	& ( pq(ig,l,i_ice) + 1.E-8 ) *
	429	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
	430	taucloudvis(ig) = taucloudvis(ig) + aerosol(ig,l,iaer)
	431	taucloudtes(ig) = taucloudtes(ig) + aerosol(ig,l,iaer)*
	432	& QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer) *
	433	& ( 1.E0 - omegaREFir3d(ig,l,iaer) )
	434	ENDDO
	435	ENDDO
	436	c 3. Outputs
	437	IF (ngrid.NE.1) THEN
	438	CALL WRITEDIAGFI(ngridmx,'tauTES','tauabs IR refwvl',
	439	& ' ',2,taucloudtes)
	440	CALL wstats(ngridmx,'tauTES','tauabs IR refwvl',
	441	& ' ',2,taucloudtes)
	442	ELSE
	443	CALL writeg1d(ngrid,1,taucloudtes,'tautes','NU')
	444	ENDIF
	445	ENDIF
	446
	447	c==================================================================
	448	END SELECT aerkind
	449
	450	c -----------------------------------------------------------------
	451	c Column integrated visible optical depth in each point
	452	c -----------------------------------------------------------------
	453
	454	do l=1,nlayer
	455	do ig=1,ngrid
	456	tau(ig,iaer) = tau(ig,iaer) + aerosol(ig,l,iaer)
	457	end do
	458	end do
	459
	460	c -----------------------------------
	461	ENDDO ! iaer (loop on aerosol kind)
	462
	463	return
	464	end

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

Download in other formats:

Original Format