Context Navigation

rocketduststorm_mod.F90 @ 2167

Last change on this file since 2167 was 2160, checked in by mvals, 6 years ago

Mars GCM:
Set adpatable parameters for the rocket dust storm scheme (parameters included in callkeys.h, and adaptable according to the callphys.def with the function "call getin" in conf_phys.F):

ti_injection, tf_injection (by default: ti_injection=10. and tf_injection=12., impacted files: compute_dtau_mod.F90, vdifc_mod.F)
coeff_detrainment (by default: coeff_detrainment=0., impacted files: rocketduststorm_mod.F90)

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

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

Download in other formats:

Original Format