Context Navigation

rocketduststorm_mod.F90 @ 2642

Last change on this file since 2642 was 2634, checked in by abierjon, 3 years ago

Mars GCM:
Changes on dust_rad_adjust for GCM6:

put an upper limit on dust_rad_adjust to avoid skyrocketing values when there are strong diurnal variations of opacity
move the condition to compute dust_rad_adjust only once per time step before the whole call to compute_dust_rad_adjust, instead of repetitive local checks
some cleaning

File size: 31.0 KB

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

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

Download in other formats:

Original Format