Context Navigation

rocketduststorm_mod.F90 @ 2929

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

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

Download in other formats:

Original Format