Context Navigation

rocketduststorm_mod.F90 @ 2228

Last change on this file since 2228 was 2226, checked in by emillour, 6 years ago
Mars GCM: Bug fix in rocketduststorm & topmons routines; "tauscaling" sould not be a local (moreover uninitialized) variable. Added it to the arguments of the routines. EM
File size: 30.4 KB

Rev	Line
[1974]	1	MODULE rocketduststorm_mod
	2
	3	IMPLICIT NONE
	4
	5	REAL, SAVE, ALLOCATABLE :: dustliftday(:) ! dust lifting rate (s-1)
	6
	7	CONTAINS
	8
	9	!=======================================================================
	10	! ROCKET DUST STORM - vertical transport and detrainment
	11	!=======================================================================
[2079]	12	! calculation of the vertical flux
[2201]	13	! call of van_leer : Van Leer transport scheme of the dust tracers
[2079]	14	! detrainement of stormdust in the background dust
[1974]	15	! -----------------------------------------------------------------------
[2079]	16	! Authors: M. Vals; C. Wang; F. Forget; T. Bertrand
[1974]	17	! Institution: Laboratoire de Meteorologie Dynamique (LMD) Paris, France
	18	! -----------------------------------------------------------------------
	19
	20	SUBROUTINE rocketduststorm(ngrid,nlayer,nq,ptime,ptimestep, &
	21	pq,pdqfi,pt,pdtfi,pplev,pplay,pzlev, &
	22	pzlay,pdtsw,pdtlw, &
	23	! input for radiative transfer
	24	clearatm,icount,zday,zls, &
	25	tsurf,igout,totstormfract, &
[2226]	26	tauscaling, &
[1974]	27	! input sub-grid scale cloud
	28	clearsky,totcloudfrac, &
[2199]	29	! input sub-grid scale topography
	30	nohmons,alpha_hmons, &
[2079]	31	! output
[2199]	32	pdqrds,wrad,dsodust,dsords, &
	33	tauref)
[1974]	34
[2079]	35	USE tracer_mod, only: igcm_stormdust_mass,igcm_stormdust_number &
[1974]	36	,igcm_dust_mass,igcm_dust_number &
	37	,rho_dust
	38	USE comcstfi_h, only: r,g,cpp,rcp
[2079]	39	USE dimradmars_mod, only: albedo,naerkind
	40	USE comsaison_h, only: dist_sol,mu0,fract
	41	USE surfdat_h, only: emis,co2ice,zmea, zstd, zsig, hmons
[2226]	42	USE callradite_mod, only: callradite
[2079]	43	IMPLICIT NONE
[1974]	44
[2160]	45	include "callkeys.h"
	46
[1974]	47	!--------------------------------------------------------
	48	! Input Variables
	49	!--------------------------------------------------------
	50
	51	INTEGER, INTENT(IN) :: ngrid ! number of horizontal grid points
	52	INTEGER, INTENT(IN) :: nlayer ! number of vertical grid points
	53	INTEGER, INTENT(IN) :: nq ! number of tracer species
	54	REAL, INTENT(IN) :: ptime
	55	REAL, INTENT(IN) :: ptimestep
	56
	57	REAL, INTENT(IN) :: pq(ngrid,nlayer,nq) ! advected field nq
	58	REAL, INTENT(IN) :: pdqfi(ngrid,nlayer,nq)! tendancy field pq
	59	REAL, INTENT(IN) :: pt(ngrid,nlayer) ! temperature at mid-layer (K)
	60	REAL, INTENT(IN) :: pdtfi(ngrid,nlayer) ! tendancy temperature at mid-layer (K)
	61
	62	REAL, INTENT(IN) :: pplay(ngrid,nlayer) ! pressure at middle of the layers
	63	REAL, INTENT(IN) :: pplev(ngrid,nlayer+1) ! pressure at intermediate levels
	64	REAL, INTENT(IN) :: pzlay(ngrid,nlayer) ! altitude at the middle of the layers
	65	REAL, INTENT(IN) :: pzlev(ngrid,nlayer+1) ! altitude at layer boundaries
	66
	67	REAL, INTENT(IN) :: pdtsw(ngrid,nlayer) ! (K/s) env
	68	REAL, INTENT(IN) :: pdtlw(ngrid,nlayer) ! (K/s) env
	69
	70	! input for second radiative transfer
[2079]	71	LOGICAL, INTENT(IN) :: clearatm
[1974]	72	INTEGER, INTENT(INOUT) :: icount
[2079]	73	REAL, INTENT(IN) :: zday
	74	REAL, INTENT(IN) :: zls
	75	REAL, INTENT(IN) :: tsurf(ngrid)
	76	INTEGER, INTENT(IN) :: igout
	77	REAL, INTENT(IN) :: totstormfract(ngrid)
[2226]	78	REAL, INTENT(INOUT) :: tauscaling(ngrid)
[2079]	79
[1974]	80	! sbgrid scale water ice clouds
	81	logical, intent(in) :: clearsky
[2199]	82	real, intent(in) :: totcloudfrac(ngrid)
	83
	84	! sbgrid scale topography
	85	LOGICAL, INTENT(IN) :: nohmons
	86	REAL, INTENT(IN) :: alpha_hmons(ngrid)
[1974]	87
	88	!--------------------------------------------------------
	89	! Output Variables
	90	!--------------------------------------------------------
	91
	92	REAL, INTENT(OUT) :: pdqrds(ngrid,nlayer,nq) ! tendancy field for dust when detraining
[2079]	93	REAL, INTENT(OUT) :: wrad(ngrid,nlayer+1) ! vertical speed within the rocket dust storm
[1974]	94	REAL, INTENT(OUT) :: dsodust(ngrid,nlayer) ! density scaled opacity of env. dust
	95	REAL, INTENT(OUT) :: dsords(ngrid,nlayer) ! density scaled opacity of storm dust
[2199]	96	REAL, INTENT(OUT) :: tauref(ngrid)
[1974]	97
	98	!--------------------------------------------------------
	99	! Local variables
	100	!--------------------------------------------------------
[2201]	101	INTEGER l,ig,iq,ll
[2079]	102	! local variables from callradite.F
[1974]	103	REAL zdtlw1(ngrid,nlayer) ! (K/s) storm
	104	REAL zdtsw1(ngrid,nlayer) ! (K/s) storm
	105	REAL zt(ngrid,nlayer) ! actual temperature at mid-layer (K)
[2091]	106	REAL zdtvert(ngrid,nlayer) ! dT/dz , lapse rate
	107	REAL ztlev(ngrid,nlayer) ! temperature at intermediate levels l+1/2 without last level
[1974]	108
[2079]	109	REAL zdtlw1_lev(nlayer),zdtsw1_lev(nlayer) ! rad. heating rate at intermediate levels l+1/2 for stormdust
	110	REAL zdtlw_lev(nlayer),zdtsw_lev(nlayer) ! rad. heating rate at intermediate levels l+1/2 for background dust
[1974]	111
[2079]	112	REAL zq_stormdust_mass(ngrid,nlayer) ! intermediate tracer stormdust mass
	113	REAL zq_stormdust_number(ngrid,nlayer) ! intermediate tracer stormdust number
	114	REAL zq_dust_mass(ngrid,nlayer) ! intermediate tracer dust mass
	115	REAL zq_dust_number(ngrid,nlayer) ! intermediate tracer dust number
[1974]	116
[2079]	117	REAL mr_stormdust_mass(ngrid,nlayer) ! intermediate mixing ratio to calculate van leer transport with the "real" concentration (stormdust mass)
	118	REAL mr_stormdust_number(ngrid,nlayer) ! intermediate mixing ratio to calculate van leer transport with the "real" concentration (stormdust number)
	119	REAL mr_dust_mass(ngrid,nlayer) ! intermediate mixing ratio to calculate van leer transport with the "real" concentration (dust mass)
	120	REAL mr_dust_number(ngrid,nlayer) ! intermediate mixing ratio to calculate van leer transport with the "real" concentration (sdust number)
	121
	122	REAL dqvl_stormdust_mass(ngrid,nlayer) ! tendancy of vertical transport (stormdust mass)
	123	REAL dqvl_stormdust_number(ngrid,nlayer) ! tendancy of vertical transport (stormdust number)
	124	REAL dqvl_dust_mass(ngrid,nlayer) ! tendancy of vertical transport (dust mass)
	125	REAL dqvl_dust_number(ngrid,nlayer) ! tendancy of vertical transport (dust number)
	126	REAL dqdet_stormdust_mass(ngrid,nlayer) ! tendancy of detrainement (stormdust mass)
	127	REAL dqdet_stormdust_number(ngrid,nlayer) ! tendancy of detrainement (stormdust number)
[1974]	128
[2090]	129	REAL masse_col(nlayer) ! mass of atmosphere (kg/m2)
[1974]	130	REAL zq(ngrid,nlayer,nq) ! updated tracers
	131
[2090]	132	REAL w(ngrid,nlayer) ! air mass flux (calculated with the vertical wind velocity profile) used as input in Van Leer (kgair/m2)
	133	REAL wqmass(ngrid,nlayer+1) ! tracer (dust_mass) mass flux in Van Leer (kg/m2)
	134	REAL wqnumber(ngrid,nlayer+1) ! tracer (dust_number) mass flux in Van Leer (kg/m2)
[1974]	135
[2079]	136	LOGICAL storm(ngrid) ! true when there is a dust storm (if the opacity is high): trigger the rocket dust storm scheme
[2160]	137	REAL detrain(ngrid,nlayer) ! coefficient for detrainment : % of stormdust detrained
[2079]	138	INTEGER scheme(ngrid) ! triggered scheme
	139
	140	REAL,PARAMETER:: coefmin =0.025 ! 0<coefmin<1 Minimum fraction of stormdust detrained
[2201]	141	REAL,PARAMETER:: wmin =0.25 ! stormdust detrainment if wrad < wmin
[2079]	142	REAL,PARAMETER:: wmax =10. ! maximum vertical velocity of the rocket dust storms (m/s)
[2201]	143
	144	! subtimestep
	145	INTEGER tsub
	146	INTEGER nsubtimestep !number of subtimestep when calling van_leer
	147	REAL subtimestep !ptimestep/nsubtimestep
	148	REAL dtmax !considered time needed for dust to cross one layer
	149	REAL,PARAMETER:: secu=3.!3. !coefficient on wspeed to avoid dust crossing many layers during subtimestep
	150
[2079]	151	! diagnostics
	152	REAL lapserate(ngrid,nlayer)
	153	REAL deltahr(ngrid,nlayer+1)
[1974]	154
[2079]	155	LOGICAL,SAVE :: firstcall=.true.
	156
	157	! variables for the radiative transfer
[1974]	158	REAL fluxsurf_lw1(ngrid)
	159	REAL fluxsurf_sw1(ngrid,2)
	160	REAL fluxtop_lw1(ngrid)
	161	REAL fluxtop_sw1(ngrid,2)
	162	REAL tau(ngrid,naerkind)
	163	REAL aerosol(ngrid,nlayer,naerkind)
	164	REAL taucloudtes(ngrid)
	165	REAL rdust(ngrid,nlayer)
	166	REAL rstormdust(ngrid,nlayer)
[2199]	167	REAL rtopdust(ngrid,nlayer)
[1974]	168	REAL rice(ngrid,nlayer)
	169	REAL nuice(ngrid,nlayer)
	170
	171
	172	! **********************************************************************
	173	! **********************************************************************
[2079]	174	! Rocket dust storm parametrization to reproduce the detached dust layers
	175	! during the dust storm season:
[1974]	176	! The radiative warming due to the presence of storm dust is
	177	! balanced by the adiabatic cooling. The tracer "storm dust"
	178	! is transported by the upward/downward flow.
	179	! **********************************************************************
	180	! **********************************************************************
	181	!! 1. Radiative transfer in storm dust
	182	!! 2. Compute vertical velocity for storm dust
[2079]	183	!! case 1 storm = false: nothing to do
	184	!! case 2 rocket dust storm (storm=true)
	185	!! 3. Vertical transport (Van Leer)
[1974]	186	!! 4. Detrainment
	187	! **********************************************************************
	188	! **********************************************************************
	189
	190
	191	! **********************************************************************
[2079]	192	! Initializations
[1974]	193	! **********************************************************************
	194	storm(:)=.false.
	195	pdqrds(:,:,:) = 0.
[2079]	196	mr_dust_mass(:,:)=0.
	197	mr_dust_number(:,:)=0.
	198	mr_stormdust_mass(:,:)=0.
	199	mr_stormdust_number(:,:)=0.
	200	dqvl_dust_mass(:,:)=0.
	201	dqvl_dust_number(:,:)=0.
	202	dqvl_stormdust_mass(:,:)=0.
	203	dqvl_stormdust_number(:,:)=0.
	204	dqdet_stormdust_mass(:,:)=0.
	205	dqdet_stormdust_number(:,:)=0.
	206	wrad(:,:)=0.
[2090]	207	w(:,:)=0.
	208	wqmass(:,:)=0.
	209	wqnumber(:,:)=0.
[2091]	210	zdtvert(:,:)=0.
[1974]	211	lapserate(:,:)=0.
	212	deltahr(:,:)=0.
	213	scheme(:)=0
[2160]	214	detrain(:,:)=1.
[2083]	215
[1974]	216	!! no update for the stormdust tracer and temperature fields
	217	!! because previous callradite was for background dust
	218	zq(1:ngrid,1:nlayer,1:nq)=pq(1:ngrid,1:nlayer,1:nq)
	219	zt(1:ngrid,1:nlayer)=pt(1:ngrid,1:nlayer)
	220
	221
[2079]	222	zq_dust_mass(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_dust_mass)
	223	zq_dust_number(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_dust_number)
	224	zq_stormdust_mass(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_stormdust_mass)
	225	zq_stormdust_number(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_stormdust_number)
	226
	227	! *********************************************************************
	228	! 0. Check if there is a storm
	229	! *********************************************************************
	230	DO ig=1,ngrid
[1974]	231	storm(ig)=.false.
[2079]	232	DO l=1,nlayer
	233	IF (zq(ig,l,igcm_stormdust_mass) &
	234	.gt. zq(ig,l,igcm_dust_mass)*(1.E-4)) THEN
[1974]	235	storm(ig)=.true.
[2079]	236	EXIT
	237	ENDIF
	238	ENDDO
	239	ENDDO
[1974]	240
	241	! *********************************************************************
	242	! 1. Call the second radiative transfer for stormdust, obtain the extra heating
	243	! *********************************************************************
[2199]	244	CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo, &
	245	emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, &
	246	zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
	247	fluxtop_sw1,tauref,tau,aerosol,dsodust,tauscaling, &
	248	taucloudtes,rdust,rice,nuice,co2ice,rstormdust,rtopdust, &
	249	totstormfract,clearatm,dsords,alpha_hmons,nohmons, &
[1974]	250	clearsky,totcloudfrac)
	251
	252	! **********************************************************************
	253	! 2. Compute vertical velocity for storm dust
[2079]	254	! **********************************************************************
[1974]	255	!! **********************************************************************
[2079]	256	!! 2.1 Nothing to do when no storm
	257	!! no storm
	258	DO ig=1,ngrid
	259	IF (.NOT.(storm(ig))) then
	260	scheme(ig)=1
	261	cycle
	262	ENDIF ! IF (storm(ig))
	263	ENDDO ! DO ig=1,ngrid
[1974]	264
[2079]	265	!! **********************************************************************
	266	!! 2.2 Calculation of the extra heating : computing heating rates
	267	!! gradient at boundaries of each layer, start from surface
	268	DO ig=1,ngrid
[2091]	269	IF (storm(ig)) THEN
	270
	271	scheme(ig)=2
[1974]	272
[2091]	273	!! computing heating rates gradient at boundraies of each layer
	274	!! start from surface
	275	zdtlw1_lev(1)=0.
	276	zdtsw1_lev(1)=0.
	277	zdtlw_lev(1)=0.
	278	zdtsw_lev(1)=0.
	279	ztlev(ig,1)=zt(ig,1)
[1974]	280
[2091]	281	DO l=1,nlayer-1
	282	!! Calculation for the dust storm fraction
	283	zdtlw1_lev(l+1)=(zdtlw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
	284	zdtlw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / &
	285	(pzlay(ig,l+1)-pzlay(ig,l))
	286
	287	zdtsw1_lev(l+1)=(zdtsw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
	288	zdtsw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / &
	289	(pzlay(ig,l+1)-pzlay(ig,l))
	290	!! Calculation for the background dust fraction
	291	zdtlw_lev(l+1)=(pdtlw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
	292	pdtlw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / &
	293	(pzlay(ig,l+1)-pzlay(ig,l))
[2079]	294
[2091]	295	zdtsw_lev(l+1)=(pdtsw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
	296	pdtsw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / &
	297	(pzlay(ig,l+1)-pzlay(ig,l))
[2079]	298
[2091]	299	ztlev(ig,l+1)=(zt(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
	300	zt(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / &
	301	(pzlay(ig,l+1)-pzlay(ig,l))
[2201]	302	ENDDO ! DO l=1,nlayer-1
[1974]	303
[2201]	304	!! This is the env. lapse rate
	305	zdtvert(ig,1)=0.
	306	DO l=1,nlayer-1
	307	zdtvert(ig,l+1)=(ztlev(ig,l+1)-ztlev(ig,l))/(pzlay(ig,l+1)-pzlay(ig,l))
	308	lapserate(ig,l+1)=zdtvert(ig,l+1)
	309	ENDDO
	310
[2091]	311	!! **********************************************************************
	312	!! 2.3 Calculation of the vertical velocity : extra heating
	313	!! balanced by adiabatic cooling
	314
[2079]	315	DO l=1,nlayer
	316	deltahr(ig,l)=(zdtlw1_lev(l)+zdtsw1_lev(l)) &
	317	-(zdtlw_lev(l)+zdtsw_lev(l))
	318	wrad(ig,l)=-deltahr(ig,l)/(g/cpp+ &
[2091]	319	max(zdtvert(ig,l),-0.99*g/cpp))
[2079]	320	!! Limit vertical wind in case of lapse rate close to adiabatic
	321	wrad(ig,l)=max(wrad(ig,l),-wmax)
	322	wrad(ig,l)=min(wrad(ig,l),wmax)
	323	ENDDO
[2091]	324
[2079]	325	ENDIF ! IF (storm(ig))
	326	ENDDO ! DO ig=1,ngrid
[1974]	327
	328	! **********************************************************************
	329	! 3. Vertical transport
	330	! **********************************************************************
[2079]	331	!! **********************************************************************
	332	!! 3.1 Transport of background dust + storm dust (concentrated)
	333	DO ig=1,ngrid
	334	IF (storm(ig)) THEN
	335	DO l=1,nlayer
	336	mr_dust_mass(ig,l) = zq_dust_mass(ig,l)
	337	mr_dust_number(ig,l) = zq_dust_number(ig,l)
	338	mr_stormdust_mass(ig,l) = zq_dust_mass(ig,l)+ &
	339	zq_stormdust_mass(ig,l)/totstormfract(ig)
	340	mr_stormdust_number(ig,l) = zq_dust_number(ig,l)+ &
	341	zq_stormdust_number(ig,l)/totstormfract(ig)
	342	ENDDO ! DO l=1,nlayer
	343	ENDIF ! IF (storm(ig))
	344	ENDDO ! DO ig=1,ngrid
[1974]	345
[2079]	346	DO ig=1,ngrid
	347	IF (storm(ig)) THEN
[2201]	348	!! **********************************************************************
	349	!! 3.2 Compute the subtimestep to conserve the mass in the Van Leer transport
	350	dtmax=ptimestep
	351	DO l=2,nlayer
	352	IF (wrad(ig,l).lt.0.) THEN
	353	dtmax=min(dtmax,(pzlev(ig,l)-pzlev(ig,l-1))/ &
	354	(secu*abs(wrad(ig,l))))
	355	ELSE IF (wrad(ig,l).gt.0.) then
	356	dtmax=min(dtmax,(pzlev(ig,l+1)-pzlev(ig,l))/ &
	357	(secu*abs(wrad(ig,l))))
	358	ENDIF
	359	ENDDO
	360	nsubtimestep= int(ptimestep/dtmax)
	361	subtimestep=ptimestep/float(nsubtimestep)
	362	!! Mass flux generated by wup in kg/m2
	363	DO l=1,nlayer
	364	w(ig,l)=wrad(ig,l)pplev(ig,l)/(rztlev(ig,l)) &
	365	*subtimestep
	366	ENDDO ! l=1,nlayer
	367
	368	!! **********************************************************************
	369	!! 3.3 Vertical transport by a Van Leer scheme
[2079]	370	!! Mass of atmosphere in the layer
	371	DO l=1,nlayer
[2090]	372	masse_col(l)=(pplev(ig,l)-pplev(ig,l+1))/g
[2079]	373	ENDDO
[2201]	374	!! Mass flux in kg/m2 if you are not using the subtimestep
	375	!DO l=1,nlayer
	376	! w(ig,l)=wrad(ig,l)(pplev(ig,l)/(rztlev(ig,l)))*ptimestep
	377	!ENDDO
	378	!! Loop over the subtimestep
	379	DO tsub=1,nsubtimestep
	380	!! Van Leer scheme
	381	wqmass(ig,:)=0.
	382	wqnumber(ig,:)=0.
	383	CALL van_leer(nlayer,mr_stormdust_mass(ig,:),2., &
	384	masse_col,w(ig,:),wqmass(ig,:))
	385	CALL van_leer(nlayer,mr_stormdust_number(ig,:),2., &
	386	masse_col,w(ig,:),wqnumber(ig,:))
	387	ENDDO !tsub=...
	388
[2079]	389	ENDIF ! IF storm(ig)
	390	ENDDO ! DO ig=1,ngrid
[1974]	391
[2079]	392	!! **********************************************************************
[2201]	393	!! 3.4 Re-calculation of the mixing ratios after vertical transport
[2079]	394	DO ig=1,ngrid
	395	IF (storm(ig)) THEN
	396	DO l=1,nlayer
	397
	398	!! General and "healthy" case
	399	IF (mr_stormdust_mass(ig,l).ge.mr_dust_mass(ig,l)) THEN
	400	zq_dust_mass(ig,l) = mr_dust_mass(ig,l)
	401	zq_dust_number(ig,l) = mr_dust_number(ig,l)
	402	zq_stormdust_mass(ig,l) = totstormfract(ig)*(mr_stormdust_mass(ig,l)-mr_dust_mass(ig,l))
	403	zq_stormdust_number(ig,l) = totstormfract(ig)*(mr_stormdust_number(ig,l)-mr_dust_number(ig,l))
	404	!! Particular case: there is less than initial dust mixing ratio after the vertical transport
	405	ELSE
	406	zq_dust_mass(ig,l) = (1.-totstormfract(ig))mr_dust_mass(ig,l)+totstormfract(ig)mr_stormdust_mass(ig,l)
	407	zq_dust_number(ig,l) = (1.-totstormfract(ig))mr_dust_number(ig,l)+totstormfract(ig)mr_stormdust_number(ig,l)
	408	zq_stormdust_mass(ig,l) = 0.
	409	zq_stormdust_number(ig,l) = 0.
	410	ENDIF
	411
	412	ENDDO ! DO l=1,nlayer
	413	ENDIF ! IF storm(ig)
	414	ENDDO ! DO ig=1,ngrid
[1974]	415
[2079]	416	!! **********************************************************************
[2201]	417	!! 3.5 Calculation of the tendencies of the vertical transport
[2079]	418	DO ig=1,ngrid
	419	IF (storm(ig)) THEN
	420	DO l=1,nlayer
	421	dqvl_stormdust_mass(ig,l) = (zq_stormdust_mass(ig,l)- &
	422	zq(ig,l,igcm_stormdust_mass)) /ptimestep
	423	dqvl_stormdust_number(ig,l) = (zq_stormdust_number(ig,l)- &
	424	zq(ig,l,igcm_stormdust_number)) /ptimestep
	425	dqvl_dust_mass(ig,l) = (zq_dust_mass(ig,l)-zq(ig,l,igcm_dust_mass)) /ptimestep
	426	dqvl_dust_number(ig,l) = (zq_dust_number(ig,l)-zq(ig,l,igcm_dust_number)) /ptimestep
	427	ENDDO
	428	ENDIF ! IF storm(ig)
	429	ENDDO ! DO ig=1,ngrid
[1974]	430
	431	! **********************************************************************
[2079]	432	! 4. Detrainment: stormdust is converted to background dust
[1974]	433	! **********************************************************************
[2079]	434	!! **********************************************************************
	435	!! 4.1 Compute the coefficient of detrainmen
	436	DO ig=1,ngrid
	437	DO l=1,nlayer
	438	IF ((max(abs(wrad(ig,l)),abs(wrad(ig,l+1))) .lt. &
	439	wmin) .or. (zq_dust_mass(ig,l) .gt. &
	440	10000.*zq_stormdust_mass(ig,l))) THEN
[2160]	441	detrain(ig,l)=1.
[2079]	442	ELSE IF (max(abs(wrad(ig,l)),abs(wrad(ig,l+1))) &
	443	.le. wmax) THEN
[2160]	444	detrain(ig,l)=coeff_detrainment* &
	445	(((1-coefmin)/(wmin-wmax)*2) &
	446	(max(abs(wrad(ig,l)),abs(wrad(ig,l+1)))-wmax)**2 &
	447	+coefmin)
[2079]	448	ELSE IF (max(abs(wrad(ig,l)),abs(wrad(ig,l+1))).gt. wmax ) THEN
[2160]	449	detrain(ig,l)=coefmin
[2079]	450	ELSE
[2160]	451	detrain(ig,l)=coefmin
[2079]	452	ENDIF
	453	ENDDO ! DO l=1,nlayer
	454	ENDDO ! DO ig=1,ngrid
	455
	456	!! **********************************************************************
	457	!! 4.2 Calculation of the tendencies of the detrainment
	458	DO ig=1,ngrid
	459	DO l=1,nlayer
[2160]	460	dqdet_stormdust_mass(ig,l)=-detrain(ig,l)*zq_stormdust_mass(ig,l) &
[1974]	461	/ptimestep
[2160]	462	dqdet_stormdust_number(ig,l)=-detrain(ig,l)*zq_stormdust_number(ig,l) &
[1974]	463	/ptimestep
[2079]	464	ENDDO ! DO l=1,nlayer
	465	ENDDO ! DO ig=1,ngrid
	466
	467	! **********************************************************************
	468	! 5. Final tendencies: vertical transport + detrainment
	469	! **********************************************************************
	470	DO ig=1,ngrid
	471	DO l=1,nlayer
	472	pdqrds(ig,l,igcm_stormdust_mass)=dqdet_stormdust_mass(ig,l) &
	473	+dqvl_stormdust_mass(ig,l)
	474	pdqrds(ig,l,igcm_stormdust_number)=dqdet_stormdust_number(ig,l) &
	475	+dqvl_stormdust_number(ig,l)
	476	pdqrds(ig,l,igcm_dust_mass)= -dqdet_stormdust_mass(ig,l) &
	477	+dqvl_dust_mass(ig,l)
	478	pdqrds(ig,l,igcm_dust_number)= -dqdet_stormdust_number(ig,l) &
	479	+dqvl_dust_number(ig,l)
	480	ENDDO ! DO l=1,nlayer
	481	ENDDO ! DO ig=1,ngrid
[1974]	482
[2201]	483	! ! **********************************************************************
	484	! ! 6. To prevent from negative values
	485	! ! **********************************************************************
	486	! DO ig=1,ngrid
	487	! DO l=1,nlayer
	488	! IF ((pq(ig,l,igcm_stormdust_mass) &
	489	! +pdqrds(ig,l,igcm_stormdust_mass)*ptimestep .le. 0.) .or. &
	490	! (pq(ig,l,igcm_stormdust_number) &
	491	! +pdqrds(ig,l,igcm_stormdust_number)*ptimestep .le. 0.)) THEN
	492	! pdqrds(ig,l,igcm_stormdust_mass)=-pq(ig,l,igcm_stormdust_mass)/ptimestep
	493	! pdqrds(ig,l,igcm_stormdust_number)=-pq(ig,l,igcm_stormdust_number)/ptimestep
	494	! ENDIF
	495	! ENDDO ! nlayer
	496	! ENDDO ! DO ig=1,ngrid
	497	!
	498	! DO ig=1,ngrid
	499	! DO l=1,nlayer
	500	! IF ((pq(ig,l,igcm_dust_mass) &
	501	! +pdqrds(ig,l,igcm_dust_mass)*ptimestep .le. 0.) .or. &
	502	! (pq(ig,l,igcm_dust_number) &
	503	! +pdqrds(ig,l,igcm_dust_number)*ptimestep .le. 0.)) THEN
	504	! pdqrds(ig,l,igcm_dust_mass)=-pq(ig,l,igcm_dust_mass)/ptimestep
	505	! pdqrds(ig,l,igcm_dust_number)=-pq(ig,l,igcm_dust_number)/ptimestep
	506	! ENDIF
	507	! ENDDO ! nlayer
	508	! ENDDO ! DO ig=1,ngrid
[2079]	509
	510	!=======================================================================
	511	! WRITEDIAGFI
	512
[1974]	513	call WRITEDIAGFI(ngrid,'lapserate','lapse rate in the storm', &
	514	& 'k/m',3,lapserate)
	515	call WRITEDIAGFI(ngrid,'deltahr','extra heating rates', &
	516	& 'K/s',3,deltahr)
	517	call writediagfi(ngrid,'scheme','which scheme',&
	518	' ',2,real(scheme))
[2079]	519
[1974]	520	END SUBROUTINE rocketduststorm
	521
[2201]	522	!=======================================================================
	523	! **********************************************************************
	524	! Subroutine to determine the vertical transport with
	525	! a Van Leer advection scheme (copied from the sedimentation scheme --> see vlz_fi.F)
	526	!***********************************************************************
	527	SUBROUTINE van_leer(nlay,q,pente_max,masse,w,wq)
[2079]	528
[1974]	529	IMPLICIT NONE
	530
[2201]	531	!--------------------------------------------------------
	532	! Input/Output Variables
	533	!--------------------------------------------------------
	534	INTEGER,INTENT(IN) :: nlay ! number of atmospheric layers
	535	REAL,INTENT(IN) :: masse(nlay) ! mass of the layer Dp/g
	536	REAL,INTENT(IN) :: pente_max != 2 conseillee
	537	REAL,INTENT(INOUT) :: q(nlay) ! mixing ratio (kg/kg)
	538	REAL,INTENT(INOUT) :: w(nlay) ! atmospheric mass "transferred" at each timestep (kg.m-2)
	539	REAL,INTENT(INOUT) :: wq(nlay+1)
	540
	541	!--------------------------------------------------------
	542	! Local Variables
	543	!--------------------------------------------------------
	544
[1974]	545	INTEGER i,l,j,ii
[2079]	546	REAL dzq(nlay),dzqw(nlay),adzqw(nlay),dzqmax
	547	REAL newmasse
	548	REAL sigw, Mtot, MQtot
	549	INTEGER m
[1974]	550
[2201]	551	! **********************************************************************
	552	! Mixing ratio vertical gradient at the levels
	553	! **********************************************************************
[1974]	554	do l=2,nlay
	555	dzqw(l)=q(l-1)-q(l)
	556	adzqw(l)=abs(dzqw(l))
	557	enddo
	558
	559	do l=2,nlay-1
	560	if(dzqw(l)*dzqw(l+1).gt.0.) then
	561	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
	562	else
	563	dzq(l)=0.
	564	endif
	565	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
	566	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
	567	enddo
	568
	569	dzq(1)=0.
	570	dzq(nlay)=0.
[2079]	571
[2201]	572	! **********************************************************************
	573	! Vertical advection
	574	! **********************************************************************
[1974]	575
[2201]	576	!! No flux at the model top:
[1974]	577	wq(nlay+1)=0.
	578
[2201]	579	!! Surface flux up:
[2119]	580	if(w(1).lt.0.) wq(1)=0. ! warning : not always valid
[1974]	581
[2201]	582	do l = 1,nlay-1
[2119]	583
[2201]	584	! 1) Compute wq where w < 0 (up) (UPWARD TRANSPORT)
[2119]	585	! ===============================
	586
[1974]	587	if(w(l+1).le.0)then
	588	! Regular scheme (transfered mass < 1 layer)
	589	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	590	if(-w(l+1).le.masse(l))then
	591	sigw=w(l+1)/masse(l)
	592	wq(l+1)=w(l+1)(q(l)-0.5(1.+sigw)*dzq(l))
[2201]	593	!!-------------------------------------------------------
	594	! The following part should not be needed in the
	595	! case of an integration over subtimesteps
	596	!!-------------------------------------------------------
[1974]	597	! Extended scheme (transfered mass > 1 layer)
	598	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	599	else
	600	m = l-1
	601	Mtot = masse(m+1)
	602	MQtot = masse(m+1)*q(m+1)
	603	if (m.le.0)goto 77
	604	do while(-w(l+1).gt.(Mtot+masse(m)))
	605	! do while(-w(l+1).gt.Mtot)
	606	m=m-1
	607	Mtot = Mtot + masse(m+1)
	608	MQtot = MQtot + masse(m+1)*q(m+1)
	609	if (m.le.0)goto 77
	610	end do
	611	77 continue
	612
	613	if (m.gt.0) then
	614	sigw=(w(l+1)+Mtot)/masse(m)
[2079]	615	wq(l+1)= -(MQtot + (-w(l+1)-Mtot)* &
[1974]	616	(q(m)-0.5(1.+sigw)dzq(m)) )
	617	else
	618	w(l+1) = -Mtot
	619	wq(l+1) = -MQtot
	620	end if
[2119]	621	endif ! (-w(l+1).le.masse(l))
[1974]	622
[2201]	623	! 2) Compute wq where w > 0 (down) (DOWNWARD TRANSPORT)
[1974]	624	! ===============================
	625
[2119]	626	else if(w(l).gt.0.)then
[1974]	627
	628	! Regular scheme (transfered mass < 1 layer)
	629	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	630	if(w(l).le.masse(l))then
	631	sigw=w(l)/masse(l)
[2201]	632	wq(l)=w(l)(q(l)+0.5(1.-sigw)*dzq(l))
	633	!!-------------------------------------------------------
	634	! The following part should not be needed in the
	635	! case of an integration over subtimesteps
	636	!!-------------------------------------------------------
[1974]	637	! Extended scheme (transfered mass > 1 layer)
	638	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	639	else
	640	m=l
	641	Mtot = masse(m)
	642	MQtot = masse(m)*q(m)
	643	if(m.ge.nlay)goto 88
	644	do while(w(l).gt.(Mtot+masse(m+1)))
	645	m=m+1
	646	Mtot = Mtot + masse(m)
	647	MQtot = MQtot + masse(m)*q(m)
	648	if(m.ge.nlay)goto 88
	649	end do
	650	88 continue
	651	if (m.lt.nlay) then
	652	sigw=(w(l)-Mtot)/masse(m+1)
	653	wq(l)=(MQtot + (w(l)-Mtot)* &
	654	(q(m+1)+0.5(1.-sigw)dzq(m+1)) )
	655	else
	656	w(l) = Mtot
	657	wq(l) = MQtot
	658	end if
[2119]	659	end if
[2100]	660
[2119]	661	end if ! w<0 (up)
[2100]	662
[2119]	663	enddo ! l = 1,nlay-1
[1974]	664
[2201]	665	do l = 1,nlay
[1974]	666
[2119]	667	! it cannot entrain more than available mass !
	668	if ( (wq(l+1)-wq(l)) .lt. -(masse(l)*q(l)) ) then
	669	wq(l+1) = wq(l)-masse(l)*q(l)
	670	end if
	671
[1974]	672	q(l)=q(l) + (wq(l+1)-wq(l))/masse(l)
	673
	674	enddo
[2079]	675
[2201]	676	END SUBROUTINE van_leer
[1974]	677
[2079]	678	!=======================================================================
	679	! Initialization of the module variables
[1974]	680
	681	subroutine ini_rocketduststorm_mod(ngrid)
	682
	683	implicit none
	684
	685	integer, intent(in) :: ngrid
	686
	687	allocate(dustliftday(ngrid))
	688
	689	end subroutine ini_rocketduststorm_mod
	690
	691	subroutine end_rocketduststorm_mod
	692
	693	implicit none
	694
	695	if (allocated(dustliftday)) deallocate(dustliftday)
	696
	697	end subroutine end_rocketduststorm_mod
	698
	699	END MODULE rocketduststorm_mod

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

Download in other formats:

Original Format