Context Navigation

rocketduststorm_mod.F90 @ 3026

Last change on this file since 3026 was 2963, checked in by romain.vande, 20 months ago

Mars PEM :

Adapt PEM to the subslope PCM configuration, it is now fully compatible.

Create a PEM folder in deftank that contains:

run_pem1: a bash file that runs chained simulation of PEM as well as running a parameterizable number of PCM simulation in between.

It also takes care of reshaping XIOS output as well as renaming outputs etc… in the spirit of run_month1.

It is written for Irene machine and the header needs to be adapted for other machines.

run_PEM.def: A text file that shows the possible parameters to choose before a PEM simulation.

It should be included at the end of run.def just like callphys.def

ob_ex_lsp.asc: An ascii file containing the obliquity, eccentricity, ls_peri data from Laskar in Martian year.
README: A txt file explaining the content of the folder

Adapt field_def_physics_mars.xml to consider the case with 7 subslopes in the PCM.
Change context_lmdz_physics.xml to be able to output the file needed by the PEM.

Correct a few other minor bugs.

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

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

Download in other formats:

Original Format