Context Navigation

rocketduststorm_mod.F90 @ 2325

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

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

Download in other formats:

Original Format