Context Navigation

rocketduststorm_mod.F90 @ 2566

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

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

Download in other formats:

Original Format