Context Navigation

aeropacity.F @ 1652

Last change on this file since 1652 was 1543, checked in by emillour, 9 years ago

All models: Further adaptations to keep up with changes in LMDZ5 concerning
physics/dynamics separation:

dyn3d:
adapted gcm.F so that all physics initializations are now done in iniphysiq.

dyn3dpar:
adapted gcm.F so that all physics initializations are now done in iniphysiq.
updated calfis_p.F to follow up with changes.
copied over updated "bands.F90" from LMDZ5.

dynphy_lonlat:
calfis_p.F90, mod_interface_dyn_phys.F90, follow up of changes in phy_common/mod_* routines

phy_common:
added "geometry_mod.F90" to store information about the grid (replaces phy*/comgeomphy.F90) and give variables friendlier names: rlond => longitude , rlatd => latitude, airephy => cell_area, cuphy => dx , cvphy => dy
added "physics_distribution_mod.F90"
updated "mod_grid_phy_lmdz.F90", "mod_phys_lmdz_mpi_data.F90", "mod_phys_lmdz_para.F90", "mod_phys_lmdz_mpi_transfert.F90", "mod_grid_phy_lmdz.F90", "mod_phys_lmdz_omp_data.F90", "mod_phys_lmdz_omp_transfert.F90", "write_field_phy.F90" and "ioipsl_getin_p_mod.F90" to LMDZ5 versions.

phy[venus/titan/mars/std]:
removed "init_phys_lmdz.F90", "comgeomphy.F90"; adapted routines to use geometry_mod (longitude, latitude, cell_area, etc.)

File size: 25.9 KB

Rev	Line
[38]	1	SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
[1353]	2	& pq,tauscaling,tauref,tau,taucloudtes,aerosol,dsodust,reffrad,
	3	& nueffrad,QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d)
[38]	4
	5	! to use 'getin'
[1036]	6	USE ioipsl_getincom, only: getin
	7	use tracer_mod, only: noms, igcm_h2o_ice, igcm_dust_mass,
[1224]	8	& igcm_dust_submicron, rho_dust, rho_ice,
	9	& nqdust
[1543]	10	use geometry_mod, only: latitude ! grid point latitudes (rad)
[1541]	11	use comgeomfi_h, only: sinlat ! sines of grid point latitudes
[1375]	12	#ifdef DUSTSTORM
[1543]	13	use geometry_mod, only: longitude
[1375]	14	use tracer_mod, only: r3n_q, ref_r0, igcm_dust_number
	15	#endif
[1226]	16	use planete_h
	17	USE comcstfi_h
[1246]	18	use dimradmars_mod, only: naerkind, name_iaer,
	19	& iaerdust,tauvis,
	20	& iaer_dust_conrath,iaer_dust_doubleq,
	21	& iaer_dust_submicron,iaer_h2o_ice
[38]	22	IMPLICIT NONE
	23	c=======================================================================
	24	c subject:
	25	c --------
	26	c Computing aerosol optical depth in each gridbox.
	27	c
	28	c author: F.Forget
	29	c ------
	30	c update F. Montmessin (water ice scheme)
	31	c and S. Lebonnois (12/06/2003) compatibility dust/ice/chemistry
	32	c update J.-B. Madeleine 2008-2009:
	33	c - added 3D scattering by aerosols;
	34	c - dustopacity transferred from physiq.F to callradite.F,
	35	c and renamed into aeropacity.F;
[607]	36	c update E. Millour, march 2012:
	37	c - reference pressure is now set to 610Pa (not 700Pa)
[38]	38	c
	39	c input:
	40	c -----
	41	c ngrid Number of gridpoint of horizontal grid
	42	c nlayer Number of layer
	43	c nq Number of tracer
	44	c zday Date (time since Ls=0, in martian days)
	45	c ls Solar longitude (Ls) , radian
	46	c pplay,pplev pressure (Pa) in the middle and boundary of each layer
	47	c pq Dust mixing ratio (used if tracer =T and active=T).
	48	c reffrad(ngrid,nlayer,naerkind) Aerosol effective radius
[1047]	49	c QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients
	50	c QREFir3d(ngrid,nlayer,naerkind) / at reference wavelengths;
	51	c omegaREFvis3d(ngrid,nlayer,naerkind) \ 3d single scat. albedo
	52	c omegaREFir3d(ngrid,nlayer,naerkind) / at reference wavelengths;
[38]	53	c
	54	c output:
	55	c -------
[607]	56	c tauref Prescribed mean column optical depth at 610 Pa
[38]	57	c tau Column total visible dust optical depth at each point
	58	c aerosol aerosol(ig,l,1) is the dust optical
	59	c depth in layer l, grid point ig
	60
	61	c
	62	c=======================================================================
	63	#include "callkeys.h"
	64
	65	c-----------------------------------------------------------------------
	66	c
	67	c Declarations :
	68	c --------------
	69	c
	70	c Input/Output
	71	c ------------
	72	INTEGER ngrid,nlayer,nq
	73
	74	REAL ls,zday,expfactor
	75	REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
	76	REAL pq(ngrid,nlayer,nq)
	77	REAL tauref(ngrid), tau(ngrid,naerkind)
	78	REAL aerosol(ngrid,nlayer,naerkind)
[1047]	79	REAL dsodust(ngrid,nlayer)
[38]	80	REAL reffrad(ngrid,nlayer,naerkind)
	81	REAL nueffrad(ngrid,nlayer,naerkind)
[1047]	82	REAL QREFvis3d(ngrid,nlayer,naerkind)
	83	REAL QREFir3d(ngrid,nlayer,naerkind)
	84	REAL omegaREFvis3d(ngrid,nlayer,naerkind)
	85	REAL omegaREFir3d(ngrid,nlayer,naerkind)
[38]	86	c
	87	c Local variables :
	88	c -----------------
	89	INTEGER l,ig,iq,i,j
	90	INTEGER iaer ! Aerosol index
[1047]	91	real topdust(ngrid)
[38]	92	real zlsconst, zp
	93	real taueq,tauS,tauN
	94	c Mean Qext(vis)/Qext(ir) profile
[1047]	95	real msolsir(nlayer,naerkind)
[38]	96	c Mean Qext(ir)/Qabs(ir) profile
[1047]	97	real mqextsqabs(nlayer,naerkind)
[38]	98	c Variables used when multiple particle sizes are used
	99	c for dust or water ice particles in the radiative transfer
	100	c (see callradite.F for more information).
[1047]	101	REAL taudusttmp(ngrid)! Temporary dust opacity
[38]	102	! used before scaling
[1047]	103	REAL tauscaling(ngrid) ! Scaling factor for qdust and Ndust
	104	REAL taudustvis(ngrid) ! Dust opacity after scaling
	105	REAL taudusttes(ngrid) ! Dust opacity at IR ref. wav. as
[38]	106	! "seen" by the GCM.
[1047]	107	REAL taucloudvis(ngrid)! Cloud opacity at visible
[38]	108	! reference wavelength
[1047]	109	REAL taucloudtes(ngrid)! Cloud opacity at infrared
[38]	110	! reference wavelength using
	111	! Qabs instead of Qext
	112	! (direct comparison with TES)
[1224]	113	REAL topdust0(ngrid)
[83]	114
[1375]	115	#ifdef DUSTSTORM
	116	!! Local dust storms
	117	logical localstorm ! =true to create a local dust storm
	118	real taulocref,ztoploc,radloc,lonloc,latloc ! local dust storm parameters
	119	real reffstorm, yeah
	120	REAL ray(ngrid) ! distance from dust storm center
	121	REAL tauuser(ngrid) ! opacity perturbation due to dust storm
	122	REAL more_dust(ngrid,nlayer,2) ! Mass mixing ratio perturbation due to the dust storm
	123	REAL int_factor(ngrid) ! useful factor to compute mmr perturbation
	124	real l_top ! layer of the storm's top
	125	REAL zalt(ngrid, nlayer) ! useful factor to compute l_top
	126	#endif
	127
[38]	128	c local saved variables
	129	c ---------------------
	130
	131	c Level under which the dust mixing ratio is held constant
	132	c when computing the dust opacity in each layer
	133	c (this applies when doubleq and active are true)
[1208]	134	INTEGER, PARAMETER :: cstdustlevel0 = 7
	135	INTEGER, SAVE :: cstdustlevel
[38]	136
[607]	137	LOGICAL,SAVE :: firstcall=.true.
[38]	138
	139	! indexes of water ice and dust tracers:
	140	INTEGER,SAVE :: i_ice=0 ! water ice
[607]	141	real,parameter :: odpref=610. ! DOD reference pressure (Pa)
[38]	142	CHARACTER(LEN=20) :: txt ! to temporarly store text
	143	CHARACTER(LEN=1) :: txt2 ! to temporarly store text
	144	! indexes of dust scatterers:
	145	INTEGER,SAVE :: naerdust ! number of dust scatterers
	146
	147	tau(1:ngrid,1:naerkind)=0
	148
	149	! identify tracers
	150
	151	IF (firstcall) THEN
	152	! identify scatterers that are dust
	153	naerdust=0
	154	DO iaer=1,naerkind
	155	txt=name_iaer(iaer)
	156	IF (txt(1:4).eq."dust") THEN
	157	naerdust=naerdust+1
	158	iaerdust(naerdust)=iaer
	159	ENDIF
	160	ENDDO
	161	! identify tracers which are dust
	162	i=0
	163	DO iq=1,nq
	164	txt=noms(iq)
	165	IF (txt(1:4).eq."dust") THEN
	166	i=i+1
	167	nqdust(i)=iq
	168	ENDIF
	169	ENDDO
	170
	171	IF (water.AND.activice) THEN
	172	i_ice=igcm_h2o_ice
	173	write(,) "aeropacity: i_ice=",i_ice
	174	ENDIF
	175
	176	c typical profile of solsir and (1-w)^(-1):
	177	msolsir(1:nlayer,1:naerkind)=0
	178	mqextsqabs(1:nlayer,1:naerkind)=0
[222]	179	WRITE(,) "Typical profiles of Qext(vis)/Qext(IR)"
	180	WRITE(,) " and Qext(IR)/Qabs(IR):"
[38]	181	DO iaer = 1, naerkind ! Loop on aerosol kind
	182	WRITE(,) "Aerosol # ",iaer
	183	DO l=1,nlayer
[1047]	184	DO ig=1,ngrid
[38]	185	msolsir(l,iaer)=msolsir(l,iaer)+
	186	& QREFvis3d(ig,l,iaer)/
	187	& QREFir3d(ig,l,iaer)
	188	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)+
	189	& (1.E0-omegaREFir3d(ig,l,iaer))**(-1)
	190	ENDDO
[1047]	191	msolsir(l,iaer)=msolsir(l,iaer)/REAL(ngrid)
	192	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)/REAL(ngrid)
[38]	193	ENDDO
	194	WRITE(,) "solsir: ",msolsir(:,iaer)
	195	WRITE(,) "Qext/Qabs(IR): ",mqextsqabs(:,iaer)
	196	ENDDO
	197
	198	! load value of tauvis from callphys.def (if given there,
	199	! otherwise default value read from starfi.nc file will be used)
	200	call getin("tauvis",tauvis)
	201
[1208]	202	IF (freedust) THEN
	203	cstdustlevel = 1
	204	ELSE
	205	cstdustlevel = cstdustlevel0
	206	ENDIF
	207
	208
[1375]	209	#ifndef DUSTSTORM
[38]	210	firstcall=.false.
[1375]	211	#endif
[38]	212
	213	END IF
	214
[607]	215	c Vertical column optical depth at "odpref" Pa
	216	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1088]	217	IF(freedust) THEN
	218	tauref(:) = 0. ! tauref is computed after, instead of being forced
	219
	220	ELSE IF(iaervar.eq.1) THEN
[1047]	221	do ig=1, ngrid
[38]	222	tauref(ig)=max(tauvis,1.e-9) ! tauvis=cste (set in callphys.def
	223	! or read in starfi
	224	end do
	225	ELSE IF (iaervar.eq.2) THEN ! << "Viking" Scenario>>
	226
	227	tauref(1) = 0.7+.3cos(ls+80.pi/180.) ! like seen by VL1
	228	do ig=2,ngrid
	229	tauref(ig) = tauref(1)
	230	end do
	231
	232	ELSE IF (iaervar.eq.3) THEN ! << "MGS" scenario >>
	233
	234	taueq= 0.2 +(0.5-0.2) (cos(0.5(ls-4.363)))**14
	235	tauS= 0.1 +(0.5-0.1) (cos(0.5(ls-4.363)))**14
	236	tauN = 0.1
	237	c if (peri_day.eq.150) then
	238	c tauS=0.1
	239	c tauN=0.1 +(0.5-0.1) (cos(0.5(ls+pi-4.363)))**14
	240	c taueq= 0.2 +(0.5-0.2) (cos(0.5(ls+pi-4.363)))**14
	241	c endif
[1047]	242	do ig=1,ngrid
[1541]	243	if (latitude(ig).ge.0) then
[1047]	244	! Northern hemisphere
	245	tauref(ig)= tauN +
[1541]	246	& (taueq-tauN)0.5(1+tanh((45-latitude(ig)180./pi)6/60))
[1047]	247	else
	248	! Southern hemisphere
	249	tauref(ig)= tauS +
[1541]	250	& (taueq-tauS)0.5(1+tanh((45+latitude(ig)180./pi)6/60))
[1047]	251	endif
	252	enddo ! of do ig=1,ngrid
[38]	253	ELSE IF (iaervar.eq.5) THEN ! << Escalier Scenario>>
	254	c tauref(1) = 0.2
	255	c if ((ls.ge.210.pi/180.).and.(ls.le.330.pi/180.))
	256	c & tauref(1) = 2.5
	257	tauref(1) = 2.5
	258	if ((ls.ge.30.pi/180.).and.(ls.le.150.pi/180.))
	259	& tauref(1) = .2
	260
	261	do ig=2,ngrid
	262	tauref(ig) = tauref(1)
	263	end do
[1278]	264	ELSE IF ((iaervar.ge.6).and.(iaervar.le.8)) THEN
	265	! clim, cold or warm synthetic scenarios
[677]	266	call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
[1502]	267	ELSE IF ((iaervar.ge.24).and.(iaervar.le.32))
[607]	268	& THEN ! << MY... dust scenarios >>
	269	call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
	270	ELSE IF ((iaervar.eq.4).or.
	271	& ((iaervar.ge.124).and.(iaervar.le.126))) THEN
	272	! "old" TES assimation dust scenario (values at 700Pa in files!)
	273	call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
[38]	274	ELSE
	275	stop 'problem with iaervar in aeropacity.F'
	276	ENDIF
	277
	278	c -----------------------------------------------------------------
	279	c Computing the opacity in each layer
	280	c -----------------------------------------------------------------
	281
	282	DO iaer = 1, naerkind ! Loop on aerosol kind
	283	c --------------------------------------------
	284	aerkind: SELECT CASE (name_iaer(iaer))
	285	c==================================================================
	286	CASE("dust_conrath") aerkind ! Typical dust profile
	287	c==================================================================
	288
	289	c Altitude of the top of the dust layer
	290	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	291	zlsconst=SIN(ls-2.76)
	292	if (iddist.eq.1) then
	293	do ig=1,ngrid
	294	topdust(ig)=topdustref ! constant dust layer top
	295	end do
	296
	297	else if (iddist.eq.2) then ! "Viking" scenario
	298	do ig=1,ngrid
[1224]	299	! altitude of the top of the aerosol layer (km) at Ls=2.76rad:
	300	! in the Viking year scenario
	301	topdust0(ig)=60. -22.sinlat(ig)*2
[38]	302	topdust(ig)=topdust0(ig)+18.*zlsconst
	303	end do
	304
	305	else if(iddist.eq.3) then !"MGS" scenario
	306	do ig=1,ngrid
	307	topdust(ig)=60.+18.*zlsconst
[1541]	308	& -(32+18zlsconst)sin(latitude(ig))**4
	309	& - 8zlsconst(sin(latitude(ig)))**5
[38]	310	end do
	311	endif
	312
	313	c Optical depth in each layer :
	314	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	315	if(iddist.ge.1) then
	316
	317	expfactor=0.
	318	DO l=1,nlayer
	319	DO ig=1,ngrid
	320	c Typical mixing ratio profile
[607]	321	if(pplay(ig,l).gt.odpref
[38]	322	$ /(988.**(topdust(ig)/70.))) then
[607]	323	zp=(odpref/pplay(ig,l))**(70./topdust(ig))
[38]	324	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
	325	else
	326	expfactor=1.e-3
	327	endif
	328	c Vertical scaling function
	329	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) *
	330	& expfactor *
	331	& QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer)
	332	ENDDO
	333	ENDDO
	334
	335	else if(iddist.eq.0) then
	336	c old dust vertical distribution function (pollack90)
	337	DO l=1,nlayer
	338	DO ig=1,ngrid
[607]	339	zp=odpref/pplay(ig,l)
	340	aerosol(ig,l,1)= tauref(ig)/odpref *
[38]	341	s (pplev(ig,l)-pplev(ig,l+1))
	342	s max( exp(.03(1.-max(zp,1.))) , 1.E-3 )
	343	ENDDO
	344	ENDDO
	345	end if
	346
	347	c==================================================================
	348	CASE("dust_doubleq") aerkind! Two-moment scheme for dust
	349	c (transport of mass and number mixing ratio)
	350	c==================================================================
	351
	352	DO l=1,nlayer
	353	IF (l.LE.cstdustlevel) THEN
	354	c Opacity in the first levels is held constant to
	355	c avoid unrealistic values due to constant lifting:
	356	DO ig=1,ngrid
	357	aerosol(ig,l,iaer) =
	358	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
	359	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
	360	& pq(ig,cstdustlevel,igcm_dust_mass) *
	361	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
[1353]	362	! DENSITY SCALED OPACITY IN INFRARED:
	363	dsodust(ig,l) =
	364	& ( 0.75 * QREFir3d(ig,cstdustlevel,iaer) /
	365	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
	366	& pq(ig,cstdustlevel,igcm_dust_mass)
[38]	367	ENDDO
	368	ELSE
	369	DO ig=1,ngrid
	370	aerosol(ig,l,iaer) =
	371	& ( 0.75 * QREFvis3d(ig,l,iaer) /
	372	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
	373	& pq(ig,l,igcm_dust_mass) *
	374	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
[1353]	375	! DENSITY SCALED OPACITY IN INFRARED:
	376	dsodust(ig,l) =
	377	& ( 0.75 * QREFir3d(ig,l,iaer) /
	378	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
	379	& pq(ig,l,igcm_dust_mass)
[38]	380	ENDDO
	381	ENDIF
	382	ENDDO
	383
	384	c==================================================================
	385	CASE("dust_submicron") aerkind ! Small dust population
	386	c==================================================================
	387
	388	DO l=1,nlayer
	389	IF (l.LE.cstdustlevel) THEN
	390	c Opacity in the first levels is held constant to
	391	c avoid unrealistic values due to constant lifting:
	392	DO ig=1,ngrid
	393	aerosol(ig,l,iaer) =
	394	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
	395	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
	396	& pq(ig,cstdustlevel,igcm_dust_submicron) *
	397	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
	398	ENDDO
	399	ELSE
	400	DO ig=1,ngrid
	401	aerosol(ig,l,iaer) =
	402	& ( 0.75 * QREFvis3d(ig,l,iaer) /
	403	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
	404	& pq(ig,l,igcm_dust_submicron) *
	405	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
	406	ENDDO
	407	ENDIF
	408	ENDDO
	409
	410	c==================================================================
	411	CASE("h2o_ice") aerkind ! Water ice crystals
	412	c==================================================================
	413
	414	c 1. Initialization
	415	aerosol(1:ngrid,1:nlayer,iaer) = 0.
	416	taucloudvis(1:ngrid) = 0.
	417	taucloudtes(1:ngrid) = 0.
	418	c 2. Opacity calculation
	419	DO ig=1, ngrid
	420	DO l=1,nlayer
	421	aerosol(ig,l,iaer) = max(1E-20,
	422	& ( 0.75 * QREFvis3d(ig,l,iaer) /
	423	& ( rho_ice * reffrad(ig,l,iaer) ) ) *
	424	& pq(ig,l,i_ice) *
	425	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
	426	& )
	427	taucloudvis(ig) = taucloudvis(ig) + aerosol(ig,l,iaer)
	428	taucloudtes(ig) = taucloudtes(ig) + aerosol(ig,l,iaer)*
	429	& QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer) *
	430	& ( 1.E0 - omegaREFir3d(ig,l,iaer) )
	431	ENDDO
	432	ENDDO
[520]	433	c 3. Outputs -- Now done in physiq.F
	434	! IF (ngrid.NE.1) THEN
[1047]	435	! CALL WRITEDIAGFI(ngrid,'tauVIS','tauext VIS refwvl',
[520]	436	! & ' ',2,taucloudvis)
[1047]	437	! CALL WRITEDIAGFI(ngrid,'tauTES','tauabs IR refwvl',
[520]	438	! & ' ',2,taucloudtes)
	439	! IF (callstats) THEN
[1047]	440	! CALL wstats(ngrid,'tauVIS','tauext VIS refwvl',
[520]	441	! & ' ',2,taucloudvis)
[1047]	442	! CALL wstats(ngrid,'tauTES','tauabs IR refwvl',
[520]	443	! & ' ',2,taucloudtes)
	444	! ENDIF
	445	! ELSE
	446	! CALL writeg1d(ngrid,1,taucloudtes,'tautes','NU')
	447	! ENDIF
[38]	448	c==================================================================
	449	END SELECT aerkind
	450	c -----------------------------------
	451	ENDDO ! iaer (loop on aerosol kind)
	452
	453	c -----------------------------------------------------------------
	454	c Rescaling each layer to reproduce the choosen (or assimilated)
	455	c dust extinction opacity at visible reference wavelength, which
[607]	456	c is originally scaled to an equivalent odpref Pa pressure surface.
[38]	457	c -----------------------------------------------------------------
	458
[1375]	459	#ifdef DUSTSTORM
	460	c -----------------------------------------------------------------
[1410]	461	! Calculate reference opacity without perturbation
[1375]	462	c -----------------------------------------------------------------
	463	IF (firstcall) THEN
	464	DO iaer=1,naerdust
	465	DO l=1,nlayer
	466	DO ig=1,ngrid
	467	tauref(ig) = tauref(ig) +
	468	& aerosol(ig,l,iaerdust(iaer))
	469	ENDDO
	470	ENDDO
	471	ENDDO
	472	tauref(:) = tauref(:) * odpref / pplev(:,1)
[1410]	473
[1375]	474	c--------------------------------------------------
[1410]	475	c Get parameters of the opacity perturbation
[1375]	476	c--------------------------------------------------
[1410]	477	iaer=1 ! just change dust
[1375]	478
	479	write(,) "Add a local storm ?"
	480	localstorm=.true. ! default value
	481	call getin("localstorm",localstorm)
	482	write(,) " localstorm = ",localstorm
	483
	484	IF (localstorm) THEN
	485	WRITE(,) "********************"
	486	WRITE(,) "ADDING A LOCAL STORM"
	487	WRITE(,) "********************"
	488
	489	write(,) "ref opacity of local dust storm"
	490	taulocref = 4.25 ! default value
	491	call getin("taulocref",taulocref)
	492	write(,) " taulocref = ",taulocref
	493
	494	write(,) "target altitude of local storm (km)"
	495	ztoploc = 10.0 ! default value
	496	call getin("ztoploc",ztoploc)
	497	write(,) " ztoploc = ",ztoploc
	498
	499	write(,) "radius of dust storm (degree)"
	500	radloc = 0.5 ! default value
	501	call getin("radloc",radloc)
	502	write(,) " radloc = ",radloc
	503
	504	write(,) "center longitude of storm (deg)"
	505	lonloc = 25.0 ! default value
	506	call getin("lonloc",lonloc)
	507	write(,) " lonloc = ",lonloc
	508
	509	write(,) "center latitude of storm (deg)"
	510	latloc = -2.5 ! default value
	511	call getin("latloc",latloc)
	512	write(,) " latloc = ",latloc
	513
	514	write(,) "reff storm (mic) 0. for background"
	515	reffstorm = 0.0 ! default value
	516	call getin("reffstorm",reffstorm)
	517	write(,) " reffstorm = ",reffstorm
	518
[1410]	519	!! LOOP: modify opacity
[1375]	520	DO ig=1,ngrid
	521
[1410]	522	!! distance to the center:
[1541]	523	ray(ig)=SQRT((latitude(ig)180./pi-latloc)*2 +
	524	& (longitude(ig)180./pi -lonloc)*2)
[1375]	525
	526	!! transition factor for storm
[1410]	527	!! factor is hardcoded -- increase it to steepen
[1375]	528	yeah = (TANH(2.+(radloc-ray(ig))*10.)+1.)/2.
	529
[1410]	530	!! new opacity field
	531	!! -- add an opacity set to taulocref
	532	!! -- the additional reference opacity will
	533	!! thus be taulocref*odpref/pplev
	534	tauuser(ig)=max( tauref(ig) * pplev(ig,1) /odpref ,
	535	& taulocref * yeah )
[1375]	536
[1410]	537	!! compute l_top
[1375]	538	DO l=1,nlayer
	539	zalt(ig,l) = LOG( pplev(ig,1)/pplev(ig,l) )
	540	& / g / 44.01
	541	& * 8.31 * 210.
	542	IF ( (ztoploc .lt. zalt(ig,l) )
	543	& .and. (ztoploc .gt. zalt(ig,l-1)) ) l_top=l-1
	544	ENDDO
	545
[1410]	546	!! change reffrad if ever needed
[1375]	547	IF (reffstorm .gt. 0.) THEN
	548	DO l=1,nlayer
	549	IF (l .lt. l_top+1) THEN
	550	reffrad(ig,l,iaer) = max( reffrad(ig,l,iaer), reffstorm
	551	& * 1.e-6 * yeah )
	552	ENDIF
	553	ENDDO
	554	ENDIF
	555
[1410]	556	ENDDO
	557	!! END LOOP
[1375]	558
[1410]	559	!! compute perturbation in each layer (equation 8 in Spiga et al. JGR 2013)
[1375]	560	DO ig=1,ngrid
	561	int_factor(ig)=0.
	562	DO l=1,nlayer
	563	IF (l .lt. l_top+1) THEN
	564	int_factor(ig) =
	565	& int_factor(ig) +
	566	& ( 0.75 * QREFvis3d(ig,l,iaer) /
	567	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
	568	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
	569	ENDIF
	570	ENDDO
	571	DO l=1, nlayer
[1410]	572	!! Mass mixing ratio perturbation due to local dust storm in each layer
[1375]	573	more_dust(ig,l,1)=
	574	& (tauuser(ig)-(tauref(ig)
	575	& * pplev(ig,1) /odpref)) /
	576	& int_factor(ig)
	577	more_dust(ig,l,2)=
	578	& (tauuser(ig)-(tauref(ig) *
	579	& pplev(ig,1) /odpref))
	580	& / int_factor(ig) *
	581	& ((ref_r0/reffrad(ig,l,iaer))**3)
	582	& * r3n_q
	583	ENDDO
	584	ENDDO
	585
[1410]	586	!! quantity of dust for each layer with the addition of the perturbation
	587	DO l=1, l_top
[1376]	588	pq(:,l,igcm_dust_mass)= pq(:,l,igcm_dust_mass)
[1410]	589	. + more_dust(:,l,1)
[1376]	590	pq(:,l,igcm_dust_number)= pq(:,l,igcm_dust_number)
[1410]	591	. + more_dust(:,l,2)
	592	ENDDO
	593	ENDIF !! IF (localstorm)
	594	tauref(:)=0.
	595	ENDIF !! IF (firstcall)
[1375]	596	#endif
	597
[1088]	598	IF (freedust) THEN
	599	tauscaling(:) = 1.
	600
	601	ELSE
	602	c Temporary scaling factor
	603	taudusttmp(1:ngrid)=0.
	604	DO iaer=1,naerdust
	605	DO l=1,nlayer
	606	DO ig=1,ngrid
	607	c Scaling factor
	608	taudusttmp(ig) = taudusttmp(ig) +
	609	& aerosol(ig,l,iaerdust(iaer))
	610	ENDDO
[38]	611	ENDDO
	612	ENDDO
[358]	613
[1088]	614	c Saved scaling factor
	615	DO ig=1,ngrid
	616	tauscaling(ig) = tauref(ig) *
	617	& pplev(ig,1) / odpref / taudusttmp(ig)
	618	ENDDO
[358]	619
[1410]	620	ENDIF ! IF (freedust)
[1088]	621
[358]	622	c Opacity computation
[38]	623	DO iaer=1,naerdust
	624	DO l=1,nlayer
	625	DO ig=1,ngrid
	626	aerosol(ig,l,iaerdust(iaer)) = max(1E-20,
[411]	627	& aerosol(ig,l,iaerdust(iaer))* tauscaling(ig))
[38]	628	ENDDO
	629	ENDDO
	630	ENDDO
[1088]	631
[1353]	632	DO l=1,nlayer
	633	DO ig=1,ngrid
	634	dsodust(ig,l) = max(1E-20,dsodust(ig,l)* tauscaling(ig))
	635	ENDDO
	636	ENDDO
	637
[1088]	638	IF (freedust) THEN
[1208]	639	! tauref has been initialized to 0 before.
	640	DO iaer=1,naerdust
	641	DO l=1,nlayer
	642	DO ig=1,ngrid
[1410]	643	#ifdef DUSTSTORM
	644	!! recalculate opacity because storm perturbation has been added
	645	IF (firstcall) THEN
	646	aerosol(ig,l,iaer) =
	647	& ( 0.75 * QREFvis3d(ig,l,iaer) /
	648	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
	649	& pq(ig,l,igcm_dust_mass) *
	650	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
	651	ENDIF
	652	#endif
[1208]	653	tauref(ig) = tauref(ig) +
	654	& aerosol(ig,l,iaerdust(iaer))
	655	ENDDO
[1088]	656	ENDDO
	657	ENDDO
[1208]	658	tauref(:) = tauref(:) * odpref / pplev(:,1)
[1088]	659	ENDIF
[411]	660
	661	c output for debug
[420]	662	c IF (ngrid.NE.1) THEN
[1047]	663	c CALL WRITEDIAGFI(ngrid,'taudusttmp','virtual tau dust',
[420]	664	c & '#',2,taudusttmp)
[1047]	665	c CALL WRITEDIAGFI(ngrid,'tausca','tauscaling',
[420]	666	c & '#',2,tauscaling)
	667	c ELSE
[1047]	668	c CALL WRITEDIAGFI(ngrid,'taudusttmp','virtual tau dust',
[420]	669	c & '#',0,taudusttmp)
[1047]	670	c CALL WRITEDIAGFI(ngrid,'tausca','tauscaling',
[420]	671	c & '#',0,tauscaling)
	672	c ENDIF
[38]	673	c -----------------------------------------------------------------
	674	c Column integrated visible optical depth in each point
	675	c -----------------------------------------------------------------
	676	DO iaer=1,naerkind
	677	do l=1,nlayer
	678	do ig=1,ngrid
	679	tau(ig,iaer) = tau(ig,iaer) + aerosol(ig,l,iaer)
	680	end do
	681	end do
	682	ENDDO
[1375]	683
	684	#ifdef DUSTSTORM
	685	IF (firstcall) THEN
	686	firstcall=.false.
	687	ENDIF
	688	#endif
	689
[38]	690	c -----------------------------------------------------------------
	691	c Density scaled opacity and column opacity output
	692	c -----------------------------------------------------------------
[420]	693	c dsodust(1:ngrid,1:nlayer) = 0.
	694	c DO iaer=1,naerdust
[1047]	695	c DO l=1,nlayer
[420]	696	c DO ig=1,ngrid
	697	c dsodust(ig,l) = dsodust(ig,l) +
	698	c & aerosol(ig,l,iaerdust(iaer)) * g /
	699	c & (pplev(ig,l) - pplev(ig,l+1))
	700	c ENDDO
	701	c ENDDO
	702	c IF (ngrid.NE.1) THEN
	703	c write(txt2,'(i1.1)') iaer
[1047]	704	c call WRITEDIAGFI(ngrid,'taudust'//txt2,
[420]	705	c & 'Dust col opacity',
	706	c & ' ',2,tau(1,iaerdust(iaer)))
	707	c IF (callstats) THEN
[1047]	708	c CALL wstats(ngrid,'taudust'//txt2,
[420]	709	c & 'Dust col opacity',
	710	c & ' ',2,tau(1,iaerdust(iaer)))
	711	c ENDIF
	712	c ENDIF
	713	c ENDDO
[38]	714
[420]	715	c IF (ngrid.NE.1) THEN
[1047]	716	c CALL WRITEDIAGFI(ngrid,'dsodust','tau*g/dp',
[38]	717	c & 'm2.kg-1',3,dsodust)
[420]	718	c IF (callstats) THEN
[1047]	719	c CALL wstats(ngrid,'dsodust',
[420]	720	c & 'tau*g/dp',
	721	c & 'm2.kg-1',3,dsodust)
	722	c ENDIF
	723	c ELSE
	724	c CALL WRITEDIAGFI(ngrid,"dsodust","dsodust","m2.kg-1",1,
	725	c & dsodust)
	726	c ENDIF ! of IF (ngrid.NE.1)
[38]	727	c -----------------------------------------------------------------
	728	return
	729	end

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

Download in other formats:

Original Format