Context Navigation

aeropacity.F @ 706

Last change on this file since 706 was 677, checked in by emillour, 13 years ago

Mars GCM:

Added the possibility to use the "cold" (iaervar=6) or "warm" (iaervar=7) dust scenarios (derived from the 7 MY24 to MY30 scenarios) in aeropacity.F and read_dust_scenario.F

File size: 18.5 KB

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

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

Download in other formats:

Original Format