Context Navigation

rocketduststorm_mod.F90 @ 2199

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

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

Download in other formats:

Original Format