Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

Iondiff_red.F @ 3567

Last change on this file since 3567 was 3035, checked in by amartinez, 16 months ago

======= VENUS PCM ===============

COMMIT BY ANTOINE MARTINEZ

Implementation of vertical ambipolar diffusion in physics

=================================

NEW KEYWORD OF PHYSIQ.DEF

=================================

NEW VERSION OF "physiq.def"

deftank/physiq-96x96x90-chemistry-IONOSPHERE-IONDIFF.def

ok_iondiff: keyword supposed to activate ion ambipolar diffusion

nbapp_chem: replaces nbapp_chim, in order to characterize the Number of calls to the chemistry routines (per Venusian day)

================

phyvenus

================

Iondiff_red.F

Calculation of the Ion Ambipolar Diffusion for 13 ions on 14:

O2+, O+, C+, N+, CO2+,
NO+, CO+, H+, H2O+, H3O+,

HCO+, N2+, OH+

The ion temperature is fixed as the half of the electron temperature
calculated in ion_h.F for the stability of the program and because the
ion temperature is lower than electron temperature.

plasmaphys_venus_mod.F90

parameters of the ambipolar diffusion scheme:

parameter (Pdiffion = 7.0D-4) ! pressure in Pa below which ion diffusion is computed
parameter (naltvert = 168) ! number of level on the altimetric subgrid
parameter (nsubvert = 24000) ! ptimephysiq/nsubvert - minimum sub-timestep allowed
parameter ( nsubmin = 40) ! ptimephysiq/nsubmin - maximum sub-timestep allowed

physic_mod.F

nbapp_chem is not fixed anymore here but deftank/in physiq.def
Ambipolar diffusion activated if (ok_iondiff) is true

conf_phys.f90

add ok_iondiff as parameters to read from physiq.def file set to .false. by default
add nbapp_chem as parameters to characterize the Number of calls to the chemistry routines (per Venusian day) instead of be fixed in physic_mod.F

to read from physiq.def file set to 24000 by default

cleshphys.h

added ok_iondiff & nbapp_chem in COMMON/clesphys_l/
removed nbapp_chim

phytrac_chemistry.F

Added security in the calculation of the sza_local in order to avoid the rare case where the |range| is above 1

File size: 64.4 KB

Rev	Line
[3035]	1	! SUBROUTINE plasmaphys(q,rhoN,P,pk,teta,mmean,mudyn,
	2	! & tempk,tetak,tempe,tetae,wcontk,wcovk,Efieldz,wphysk,alt)
	3
	4	SUBROUTINE iondiff_red(ngrid,nlayer,nqmx,pplay,pplev,
	5	! & ptneutre,ptelect,ption,pq,
	6	& ptneutre,pq,pphis,
	7	& gmtime,lat,lon,
	8	& ptimestep,pdteuv,pdtconduc,pdqdiff,
	9	& pdqiondiff)
	10
	11	!****************************************************************
	12	!
	13	! Ion ambipolar diffusion routine
	14	!
	15	! Authors: J-Y. Chaufray
	16	! -------
	17	!
	18	! J-Y. Chaufray
	19	!
	20	! Vitesse verticale et ses effets sur nk, Tk, uk, vk, etc...
	21	! Seule les effets du terme (wk - wn) sont pris en compte ici
	22	! Pour l instant seul les effets sur la densite nk sont pris en compte
	23	!
	24	! Version: 14/11/2020
	25	!
	26	! -------
	27	!
	28	! 2022/10/22: adapted and adding by Antoine Martinez
	29	!
	30	!*****************************************************************
	31
	32	USE chemparam_mod
	33	USE infotrac_phy, only: tname
	34	!USE dimphy
	35
	36	use plasmaphys_venus_mod
	37	use iono_h, only: temp_elect, temp_ion
	38
	39	IMPLICIT NONE
	40
	41	! Solve the vertical dynamics equation
	42	! The equation is solved on a altitude refined grid after interpolation
	43	! final results are reinterpolated on the atmospheric pressure grids
	44	! Upper boundary
	45	! The velocity at upper altitude should be parametrized (escape)
	46	! Here I assume w2up (nlayer)=0
	47	! Lower boundaru
	48	! W2(1) = O (The ion velocity is equal to neutral velocity)
	49
	50	!#include "dimensions.h"
	51	!#include "paramet.h"
	52	!#include "plasmadyn.h"
	53	!#include "comgeom.h"
	54
	55	! =========================
	56	! ==== INPUT / OUTPUT =====
	57	! =========================
	58
	59	! ====== INPUT ======
	60
	61	INTEGER,intent(in) :: ngrid ! number of atmospheric columns
	62	INTEGER,intent(in) :: nlayer ! number of atmospheric layers
	63	INTEGER,intent(in) :: nqmx ! number of advected tracers
	64	REAL,intent(in) :: ptimestep ! physical timestep (s)
	65	REAL,intent(in) :: pphis(ngrid) ! geopotential of the surface (m2/s)
	66
	67	REAL,intent(in) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
	68	REAL,intent(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
	69	REAL,intent(in) :: pq(ngrid,nlayer,nqmx) ! mid-layer mass mixing ratio tracers (kg/kg)
	70
	71	REAL,intent(in) :: ptneutre(ngrid,nlayer) ! mid-layer NEUTRAL temperature (K)
	72	! REAL pdt(ngrid,nlayer)
	73	! REAL,intent(in) :: ptelect(ngrid,nlayer) ! mid-layer ELECTRON temperature (K)
	74	! REAL,intent(in) :: ption(ngrid,nlayer) ! mid-layer ION temperature (K)
	75
	76	REAL,intent(in) :: pdteuv(ngrid,nlayer) ! Temperature EUV heating tendancy (K/s)
	77	REAL,intent(in) :: pdtconduc(ngrid,nlayer) ! Temperature conduction tendancy (K/s)
	78	REAL,intent(in) :: pdqdiff(ngrid,nlayer,nqmx) ! neutral mmr molecular diffusion tendancy (kg/kg/s)
	79	! real pdq(ngrid,nlayer,nq)
	80
	81	REAL,intent(in) :: gmtime ! day fraction ratio [from 0 to 1]
	82	REAL,intent(in) :: lat(ngrid) ! grid latitude [degree]
	83	REAL,intent(in) :: lon(ngrid) ! grid longitude [degree]
	84
	85	! ====== OUTPUT ======
	86
	87	REAL :: pdqiondiff(ngrid,nlayer,nqmx) ! ion mmr tendancy (kg/kg/s)
	88
	89	! =========================
	90	! ======== LOCAL ==========
	91	! =========================
	92
	93	integer, parameter :: t_elec_origin= 2
	94	integer, parameter :: t_ion_origin = 1
	95	!Electronic temperature. Index 1 -> Theis et al. 1980 - model data ; Index 2-> Theis et al. 1984 - model data
	96	!Ion temperature. Index 1 -> from VIRA Model based on the model of Miller et al. 1980 & 1984.
	97
	98	REAL :: lon_sun, sza_local,szad_local,mu_local
	99	REAL :: lon_local(ngrid), lat_local(ngrid)
	100
	101	REAL :: ptelect(ngrid,nlayer) ! mid-layer ELECTRON temperature (K)
	102	REAL :: ption(ngrid,nlayer) ! mid-layer ION temperature (K)
	103
	104	REAL :: qnew(ngrid,nlayer,nqmx), qold(ngrid,nlayer,nqmx)
	105
	106	REAL :: mmean_new(ngrid,nlayer), vmr_new(ngrid,nlayer,nqmx)
	107	REAL :: vmr_ion(ngrid,nlayer)
	108
	109	!REAL wcontk(ip1jmp1,llp+1,ndiffions),Efieldz(ip1jmp1,llp+1)
	110	!REAL wcovk(ip1jmp1,llp+1,ndiffions),wphysk(ip1jmp1,llp+1,ndiffions)
	111
	112	! 1D field I pass to double precision
	113	REAL(kind=kind(1.D0)) tdiff1,tdiff2
	114	REAL(kind=kind(1.D0)),dimension(nlayer+1) :: P1D,Ez1d
	115	REAL(kind=kind(1.D0)),dimension(nlayer) :: Mmean1D,RhoN1D,qe1D
	116	REAL(kind=kind(1.D0)),dimension(nlayer) :: Pnlay1D,ZZ1D,TempN1d
	117	REAL(kind=kind(1.D0)),dimension(nlayer) :: TempE1D,TempE1dint
	118	REAL(kind=kind(1.D0)),dimension(nlayer) :: FacTe,TempE1dnew
	119
	120	REAL(kind=kind(1.D0)),dimension(naltvert) :: Zraf,Praf,Tnraf,Mraf
	121	REAL(kind=kind(1.D0)),dimension(naltvert) :: Rraf,QEraf,REraf
	122	REAL(kind=kind(1.D0)),dimension(naltvert) :: RErafold,TEraf,Ez
	123
	124	!REAL(kind=kind(1.D0)) X(lld),Y(lld),a(lld),b(lld),c(lld)
	125	!REAL(kind=kind(1.D0)) X2(lld),Y2(lld),a2(lld),b2(lld),c2(lld)
	126
	127	INTEGER,dimension(1) :: indZ
	128
	129	REAL(kind=kind(1.D0)) :: Rho0,Nu0,T0,H0,D0,dt,dtad,time0
	130
	131	REAL(kind=kind(1.D0)) :: tdiff3,tmin,Wmax,ttot,tsim,tdiff
	132
	133	REAL(kind=kind(1.D0)) :: pphis1D ! geopotential of the surface (m2/s)
	134
	135	REAL(kind=kind(1.D0)),save :: Wlim
	136
	137	REAL(kind=kind(1.D0)),dimension(naltvert) :: RIad,Tnad,TIad,Tead
	138	REAL(kind=kind(1.D0)),dimension(naltvert) :: RElad,Dad,Zad
	139	REAL(kind=kind(1.D0)) :: Dzraf,Dz,FacEsc,Time
	140	REAL(kind=kind(1.D0)),dimension(naltvert) :: alpha,beta,delta
	141	REAL(kind=kind(1.D0)),dimension(naltvert) :: eps,gama,dzeta,eta
	142	REAL(kind=kind(1.D0)),dimension(naltvert) :: Atri,Btri,Ctri
	143	REAL(kind=kind(1.D0)),dimension(naltvert) :: Dtri,Xtri,Ytri
	144	INTEGER :: ig,k,l,ln,kn,ij0,l0,ln0,iter,ld,iq,istep,ke,niter,nn
	145	INTEGER :: iqelec
	146	INTEGER :: nsubreal
	147	LOGICAL,SAVE :: firstcall=.true.
	148
	149	REAL(kind=kind(1.D0)) masse,invsgmu,Konst
	150	REAL(kind=kind(1.D0)) masseU,kBolt,g,RgazP,Rvenus,Grav,Mvenus,PI
	151
	152	integer,save :: ntracers
	153	integer,dimension(:),allocatable,save :: gcmind ! array of GCM indexes
	154	real,dimension(:),allocatable,save :: mol_tr ! array of mol mass traceurs
	155	real,dimension(:),allocatable,save :: Charge_tr ! array of electric charge traceurs
	156
	157	integer,dimension(nqmx) :: indic_tracers,indic_iondiff
	158	integer,parameter :: nb_espIon_diff=24
	159	character(len=20),dimension(nb_espIon_diff) :: ListeIonDiff
	160	character(len=20) :: electrans(1)
	161
	162	integer,dimension(:),allocatable,save :: itrans ! array of sub index gcmind of fluid ions
	163	integer,dimension(1),save :: etrans ! index of gcmind electron
	164
	165	INTEGER,save :: iz_vertplasma ! vertical index above which ion diffusion is computed
	166	INTEGER,save :: ndiffions ! nombre d'ions decrit
	167
	168	! ---------- ALLOCATABLE ARRAY ------------------------------
	169
	170	! (ngrid,nlayer,ntracers)
	171	!REAL(kind=kind(1.D0)),dimension(:,:),allocatable :: qnew
	172	!REAL(kind=kind(1.D0)),dimension(:,:),allocatable :: qold
	173
	174	! (nlayer,ntracers)
	175	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: q1D
	176	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: RhoK1D
	177
	178	! (nlayer,ndiffions)
	179	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: qk1d
	180	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: qk1dnew
	181	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: qk1dint
	182	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: TempI1D
	183	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: TempI1din
	184	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: TempI1dnw
	185	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: FacQ
	186	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: FacTI
	187
	188	! (nlayer+1,ndiffions)
	189	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Wk1d
	190	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Wk1d2
	191	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Wk1d3
	192
	193	! (ngrid,ndiffions)
	194	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: ueff
	195
	196	! (ndiffions)
	197	REAL(kind=kind(1.D0)),dimension(:),allocatable,save :: Mtot
	198	REAL(kind=kind(1.D0)),dimension(:),allocatable,save :: Mtot2
	199
	200	! ---------- ALLOCATABLE ARRAY ------------------------------
	201
	202	! (naltvert,ntracers)
	203	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Qraf
	204	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Rkraf
	205
	206	! (naltvert,ndiffions)
	207	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: QIraf
	208	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: RIraf
	209	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: RIrafold
	210	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: TIraf
	211	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: FIraf
	212	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: DIraf
	213	REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Rhock
	214
	215	! (ndiffions)
	216	REAL(kind=kind(1.D0)),dimension(:),allocatable,save :: MtotZ1
	217	REAL(kind=kind(1.D0)),dimension(:),allocatable,save :: MtotZ2
	218
	219	! ==========================================================
	220	! ==========================================================
	221	! ========== FIRST CALL AND MAIN PARAMETERS ================
	222	! ==========================================================
	223	! ==========================================================
	224
	225	! ======= Initializations at first call ==========
	226
	227	if (firstcall) then
	228
	229	! === On elimine la microphysique ==========
	230
	231	ntracers=0
	232	indic_tracers(1:nqmx)=0
	233
	234	! type_tr ==> 1: neutral gas, 2: ion gas, 3: liquid, 10: others
	235	do iq=1,nqmx
	236	if ((type_tr(iq) .le. 2.) .and. (type_tr(iq) .gt. 0.)) then
	237	indic_tracers(iq)=1
	238	ntracers=ntracers+1
	239	endif
	240	enddo
	241	allocate(gcmind(ntracers))
	242	allocate(mol_tr(ntracers))
	243	allocate(Charge_tr(ntracers))
	244
	245	! Store gcm indexes in gcmind and molar masses in mol_tr
	246	nn=0
	247	do iq=1,nqmx
	248	if (indic_tracers(iq) .eq. 1) then
	249	nn=nn+1
	250	gcmind(nn) = iq
	251	mol_tr(nn) = M_tr(iq)
	252	Charge_tr(nn) = 0.D0
	253	endif
	254	enddo
	255
	256	! === On selectionne les especes ioniques ==========
	257
	258	ndiffions = 0
	259	indic_iondiff(1:nqmx) = 0
	260
	261	! define used ions
	262	ListeIonDiff(1)="o2plus"
	263	ListeIonDiff(2)="oplus"
	264	ListeIonDiff(3)="cplus"
	265	ListeIonDiff(4)="nplus"
	266	ListeIonDiff(5)="co2plus"
	267	ListeIonDiff(6)="noplus"
	268	ListeIonDiff(7)="coplus"
	269	ListeIonDiff(8)="hplus"
	270	ListeIonDiff(9)="h2oplus"
	271	ListeIonDiff(10)="h3oplus"
	272	ListeIonDiff(11)="hcoplus"
	273	ListeIonDiff(12)="n2plus"
	274	ListeIonDiff(13)="ohplus"
	275
	276	! define electron
	277	electrans(1)="elec"
	278
	279	! seek the index and number of ion/electron
	280	do iq=1,ntracers
	281	do nn=1,nb_espIon_diff
	282	if (ListeIonDiff(nn) .eq. tname(gcmind(iq))) then
	283	indic_iondiff(iq)=iq ! 1
	284	endif
	285	! we don't diffuse electron but we readapt electron density at the end
	286	if (electrans(1) .eq. tname(gcmind(iq))) then
	287	etrans(1)=iq
	288	endif
	289	enddo
	290	if (indic_iondiff(iq) .ge. 1) then
	291	print*,'Ion specie ', tname(gcmind(iq)),
	292	& 'diffused in Iondiff_red'
	293	ndiffions=ndiffions+1
	294	endif
	295	enddo
	296	print*,'number of diffused ion:',ndiffions
	297	allocate(itrans(ndiffions))
	298
	299	! Store gcm ion indexes in itrans
	300	nn=0
	301	do iq=1,ntracers
	302	if (indic_iondiff(iq) .ge. 1) then
	303	nn=nn+1
	304	itrans(nn)=indic_iondiff(iq) ! iq
	305	Charge_tr(itrans(nn))=1.!qcharge
	306	endif
	307	enddo
	308	Charge_tr(etrans(1))=-1.!qcharge*(-1.)
	309
	310	! find vertical index above which ion diffusion is computed
	311	do l=1,nlayer
	312	if (pplay(1,l) .gt. Pdiffion) then
	313	iz_vertplasma = l
	314	endif
	315	enddo
	316	iz_vertplasma = iz_vertplasma+1 ! seuil de la diffusion verticale
	317	print*,'vertical index for ion diffusion',
	318	& iz_vertplasma,pplay(1,iz_vertplasma)
	319	! iz_plasma = iz_vertplasma-3 ! seuil de la dynamique horizontale
	320
	321	! allocatation des tableaux dependants du nombre d especes diffusees
	322	allocate(ueff(ngrid,ndiffions))
	323
	324	! ----------------------------------- !
	325
	326	!allocate-(qnew(ngrid,nlayer,ntracers))
	327	!allocate-(qold(ngrid,nlayer,ntracers))
	328
	329	allocate(q1D(nlayer,ntracers))
	330	allocate(RhoK1D(nlayer,ntracers))
	331
	332	allocate(qk1d(nlayer,ndiffions))
	333	allocate(qk1dnew(nlayer,ndiffions))
	334	allocate(qk1dint(nlayer,ndiffions))
	335	allocate(TempI1D(nlayer,ndiffions))
	336	allocate(TempI1din(nlayer,ndiffions))
	337	allocate(TempI1dnw(nlayer,ndiffions))
	338	allocate(FacQ(nlayer,ndiffions))
	339	allocate(FacTI(nlayer,ndiffions))
	340
	341	allocate(Wk1d(nlayer+1,ndiffions))
	342	allocate(Wk1d2(nlayer+1,ndiffions))
	343	allocate(Wk1d3(nlayer+1,ndiffions))
	344
	345	allocate(Mtot(ndiffions))
	346	allocate(Mtot2(ndiffions))
	347
	348	! ----------------------------------- !
	349	allocate(Qraf(naltvert,ntracers))
	350	allocate(Rkraf(naltvert,ntracers))
	351
	352	allocate(QIraf(naltvert,ndiffions))
	353	allocate(RIraf(naltvert,ndiffions))
	354	allocate(RIrafold(naltvert,ndiffions))
	355	allocate(TIraf(naltvert,ndiffions))
	356	allocate(FIraf(naltvert,ndiffions))
	357	allocate(DIraf(naltvert,ndiffions))
	358	allocate(Rhock(naltvert,ndiffions))
	359
	360	allocate(MtotZ1(ndiffions))
	361	allocate(MtotZ2(ndiffions))
	362
	363	! ----------------------------------- !
	364
	365	! Conditions aux limites
	366
	367	Wlim=500. ! Maximal absolute value of vertical velocity (m/s) (to avoid unstabilities in wdu/dz term)
	368	ueff(1:ngrid,1:ndiffions)=0. ! Effusion velocity at top in m/s
	369	!ueff(1:ngrid,1:ndiffions)=200.
	370
	371	firstcall=.FALSE.
	372	endif ! firstcall
	373
	374	masseU = 1.660538782d-27
	375	kBolt = 1.3806504d-23
	376	Konst = kBolt/masseU
	377	RgazP = 8.314472
	378	PI = 2.*ASIN(1.) ! 3.141592653
	379	g = 8.87d0
	380	Rvenus = 6051800d0 ! Used to compute escape parameter no need to be more accurate
	381	Grav = 6.67d-11
	382	Mvenus = 4.86d24
	383	ij0 = 6000 ! For test
	384	invsgmu= 1d0/g/masseU
	385
	386	ke=etrans(1)
	387
	388	alpha(1:naltvert)=0D0
	389	beta(1:naltvert)=0D0
	390	delta(1:naltvert)=0D0
	391	eps(1:naltvert)=0D0
	392	dzeta(1:naltvert)=0D0
	393	eta(1:naltvert)=0D0
	394	!! Open Vertical file to check it
	395	! if (firstcall) then
	396	! open(unit=301,file='Plasma_vertical.dat',status='unknown')
	397	! firstcall=.FALSE.
	398	! endif
	399	!! First we derive only parameters useful for diffusion
	400	!! Extrapolation of paramaters to all altitudes (we solve the dynamics equation on larger altitude range)
	401	!!
	402	!! ===> ON A PAS BESOIN DE CETTE PARTIE CAR ON EST DEJA DANS LA PHYSIQUE
	403	!! CALL DYN2DIFFK(Pk,teta,P,mmean,rhoN,Preff,cpp,kappa,
	404	!! & q,Mtraceur,Tempk,tetak,TempE,TetaE,Preffplasma,itrans,etrans,
	405	!! & llm,nqtot,llp,ndiffions,ip1jmp1,iz_plasma,kappak,kappae,
	406	!! & PNlay,TempN,Rhok,TempI,TetaI,rhotetaI,TempF,TetaF,rhotetaF)
	407	!!
	408
	409	! ==========================================================
	410	! ==========================================================
	411	! ========== Main Loop on horizontal grids =================
	412	! ==========================================================
	413	! ==========================================================
	414
	415	! ========== Initialization =================
	416
	417	dt=dble(ptimestep) !dtplasma*iappvert
	418	nsubreal=nsubvert
	419	niter=2
	420	qold(1:ngrid,1:nlayer,1:nqmx)=pq(1:ngrid,1:nlayer,1:nqmx)
	421	qnew(1:ngrid,1:nlayer,1:nqmx)=pq(1:ngrid,1:nlayer,1:nqmx)
	422	! tdiff=dt/nsubreal
	423	! print*,'dt',dt
	424
	425	lon_sun = (0.5 - gmtime)2.PI
	426	lon_local(:) = lon(:)*PI/180.
	427	lat_local(:) = lat(:)*PI/180.
	428
	429	! ========== Horizontal loop start ==========
	430
	431	DO ig=1,ngrid
	432
	433	sza_local = cos(lat_local(ig))*cos(lon_local(ig))
	434	& *cos(lon_sun) + cos(lat_local(ig))
	435	& sin(lon_local(ig))sin(lon_sun)
	436	sza_local = min(sza_local,1.)
	437	sza_local = max(-1.,sza_local)
	438	sza_local = acos(sza_local) ! RAD
	439	szad_local = sza_local * 180. / PI ! Degree
	440	mu_local = cos(sza_local)
	441
	442	ttot = dt
	443	tsim = 0.
	444	tdiff1 = dt/nsubvert !
	445	tdiff2 = dt/nsubmin !
	446	!if (nsubreal .ne. nsubvert) print*,'ig',ig,tdiff
	447
	448	ZZ1D(1:nlayer) = 0D0
	449	P1D(1:nlayer) = 0D0
	450	Pnlay1D(1:nlayer) = 0D0
	451	TempN1D(1:nlayer) = 0D0
	452	Mmean1D(1:nlayer) = 0D0
	453	RHON1D(1:nlayer) = 0D0
	454	RHOK1D(1:nlayer,1:ntracers)=0D0
	455	Q1D(1:nlayer,1:ntracers)=0D0
	456	TEMPI1D(1:nlayer,1:ndiffions)=0D0
	457	!TETAI1D(1:nlayer,1:ndiffions)=0D0
	458	TEMPE1D(1:nlayer)=0D0
	459	!TETAE1D(1:nlayer)=0D0
	460	Zraf(1:naltvert)=0D0
	461	Praf(1:naltvert)=0D0
	462	Tnraf(1:naltvert)=0D0
	463	Mraf(1:naltvert)=0D0
	464	Qraf(1:naltvert,1:ntracers)=0D0
	465	Rkraf(1:naltvert,1:ntracers)=0D0
	466	REraf(1:naltvert)=0D0
	467	QIraf(1:naltvert,1:ndiffions)=0D0
	468	TIraf(1:naltvert,1:ndiffions)=0D0
	469	FIraf(1:naltvert,1:ndiffions)=0D0
	470	QEraf(1:naltvert)=0D0
	471	TEraf(1:naltvert)=0D0
	472	rhock(1:naltvert,1:ndiffions)=0D0
	473
	474	! print*,'ig',ig,llm
	475	! First we compute 1d fields in double precision
	476	! The mixing ratios slightly change here (< 0.1%)
	477
	478	pphis1D = pphis(ig)
	479	P1D(nlayer+1) = pplev(ig,nlayer+1)
	480	do l=1,nlayer
	481	P1D(l) = pplev(ig,l)
	482	Pnlay1D(l) = pplay(ig,l)
	483	TempN1D(l) = ptneutre(ig,l)
	484	do iq=1,ntracers
	485	q1d(l,iq)=qold(ig,l,gcmind(iq))
	486	enddo
	487	! compute the mean molecular mass
	488	Mmean1d(l) = (1D0/sum(q1d(l,:)/dble(mol_tr(:))))
	489
	490	! Compute the total mass density (kg/m3)
	491	rhok1d(l,:) = 0D0
	492	RhoN1D(l) = Pnlay1D(l)*Mmean1d(l)/TempN1D(l)/Konst
	493	do iq=1,ntracers
	494	rhok1d(l,iq)=q1d(l,iq)*RhoN1D(l)
	495	enddo
	496
	497	do k=1,ndiffions
	498	kn=itrans(k)
	499	qk1d(l,k)=q1D(l,kn)
	500	!TempI1D(l,k)=ption(ig,l)
	501	!TetaI1D(l,k)=TetaI(ig,l,k)
	502	!rhotetaI1D(l,k)=rhotetaI(ig,l,k)
	503	enddo
	504
	505	!ke=etrans(1)
	506	qe1D(l)=q1D(l,ke)
	507	!TempE1D(l)=ptelect(ig,l)
	508	!TetaE1D(l)=TetaF(ig,l)
	509	!rhoTetaE1D(l)=rhoTetaf(ig,l)
	510	enddo ! l (vertical levels)
	511
	512	! Compute colonne mass of each dynamic ion
	513
	514	Mtot(1:ndiffions)=0D0
	515
	516	! do k=1,ndiffions
	517	! kn=itrans(k)
	518	! do l=1,lld
	519	! ln=l+iz_vertplasma
	520	! print*,ln,llm
	521	! Mtot(k)=Mtot(k)+qk1d(ln,k)*rhoN1d(ln)
	522	! enddo
	523	! enddo
	524
	525	! print*,'Mtot before',Mtot
	526
	527	! Compute altitude of the scalar grid pressure level
	528
	529	!if (isnan(TempN1d(1))) then
	530	! write (,) 'PROBLEME NAN: TITI'
	531	!endif
	532
	533	CALL ZVERTK(Pnlay1D,TempN1D,mmean1d,ZZ1d,nlayer,pphis1D)
	534
	535	! Compute lectron and Ion temperature of the scalar grid pressure level
	536
	537	do l=1,nlayer
	538	ptelect(ig,l) = temp_elect(ZZ1d(l)/1000.,TempN1d(l),
	539	& szad_local,t_elec_origin)
	540	ption(ig,l) = temp_ion(ZZ1d(l)/1000.,TempN1d(l),
	541	& szad_local,t_ion_origin)
	542
	543	TempE1D(l) = ptelect(ig,l)
	544	!!! The ion temperature is fixed as the half of the electron temperature
	545	!!! calculated in ion_h.F for the stability of the program and because the
	546	!!! ion temperature is lower than electron temperature.
	547	do k=1,ndiffions
	548	!kn=itrans(k)
	549	TempI1D(l,k)=max(0.5*ptelect(ig,l),TempN1d(l))
	550	enddo
	551	enddo ! l (vertical levels)
	552	!write(,) szad_local,ZZ1d(:)/1000.,ptelect(ig,:), ption(ig,:)
	553	! do l=1,llm
	554	! print*,'Pnlay1D, Plev',Pnlay1D(l),P1D(l),rhok1D(l,1:ntracers)
	555	! enddo
	556	! I use a better vertial resolution above iz_vertplasma and altitude grid interpolation
	557
	558	indZ(1)=0
	559	CALL UPPER_RESOLK(ZZ1D,Pnlay1D,TempN1d,mmean1d,rhok1d,TempI1d,
	560	& TempE1D,mol_tr,nlayer,ndiffions,ntracers,naltvert,iz_vertplasma,
	561	& itrans,etrans,Zraf,Tnraf,Praf,Mraf,Qraf,Rkraf,
	562	& QIraf,TIraf,FIraf,Qeraf,Teraf,indZ,gcmind)
	563
	564	! print*,'z last level, P last level',Zraf(:),Praf(:)
	565
	566	DZraf=Zraf(2)-Zraf(1)
	567
	568	do l=1,naltvert
	569	do k=1,ndiffions
	570	kn=itrans(k)
	571	RIraf(l,k)=Rkraf(l,kn)
	572	enddo
	573	!ke=etrans(1)
	574	REraf(l)=Rkraf(l,ke)
	575	enddo
	576
	577	! Test
	578	! do l=1,naltvert
	579	! print*,'l',l,Zraf(l),Praf(l),Tnraf(l),Qiraf(l,1),Tiraf(l,1),
	580	! & Qeraf(l),Teraf(l),Riraf(l,1),Reraf(l)
	581	! enddo
	582
	583
	584	! Before solving the diffusion equation, we reddo interpolation on GCM grid
	585	! to derive mass correction factor due to interpolation (for mass conservation purpose)
	586
	587	CALL GCMGRID_PK(Zraf,Praf,TNraf,Mraf,
	588	& QIraf,RIraf,TIraf,Teraf,
	589	& P1D,Pnlay1D,TempN1d,Mmean1d,
	590	& qk1d,qk1dint,TempI1D,TempI1Din,TempE1d,TempE1dint,
	591	& naltvert,ndiffions,nlayer,iz_vertplasma,itrans,etrans)
	592
	593	!DO ln=1,nlayer
	594	!IF (k.eq.1) THEN
	595	! write(,) real(qk1Dint(ln,:)), qold(ig,ln,gcmind(ke))
	596	!ENDIF
	597	!ENDDO
	598
	599	FacQ(1:nlayer,1:ndiffions)=1D0
	600	FacTI(1:nlayer,1:ndiffions)=1D0
	601	FacTE(1:nlayer)=1D0
	602
	603	CALL CORRFACK(qk1d,qk1dint,TempI1d,TempI1din,TempE1d,
	604	& TempE1dint,FacQ,FacTI,FacTE,nlayer,ndiffions)
	605
	606	CALL DCOEFFK(TIraf,FIraf,mol_tr,DIraf,
	607	& naltvert,ndiffions,ntracers,itrans)
	608
	609	Rho0=1d-20
	610
	611	T0=Tnraf(naltvert)!TIraf(naltvert,1)
	612
	613	do ln=1,naltvert
	614	Tnad(ln)=TNraf(ln)/T0
	615	enddo
	616
	617	! if(abs(lon(ig)-(48.75)).le.0.1) then
	618	! if(abs(lat(ig)-(-80.625)).le.0.1) then
	619	! write(,) 'ZZ1d= ', ZZ1d(iz_vertplasma:90)
	620	! write(,) 'Zraf= ', Zraf(1:naltvert)
	621	! write(,) 'Pnlay1D= ', Pnlay1D(iz_vertplasma:90)
	622	! write(,) 'Praf= ', Praf(1:naltvert)
	623	! DO k=1,ndiffions
	624	! kn=itrans(k)
	625	! if((tname(gcmind(kn)).eq.'hplus').or.
	626	! & (tname(gcmind(kn)).eq.'1oplus').or.
	627	! & (tname(gcmind(kn)).eq.'1noplus')) then
	628	! write(,) tname(gcmind(kn))//' DEBUT'
	629	! write(,) 'avant= ',
	630	! & qk1d(iz_vertplasma:90,k)
	631	! write(,) 'apres= ',
	632	! & qk1dint(iz_vertplasma:90,k)
	633	! write(,) 'interm= ',
	634	! & Rkraf(1:naltvert,kn)
	635	! write(,) 'corrmass= ',
	636	! & real(FacQ(iz_vertplasma:90,k))
	637	! endif
	638	! ENDDO
	639	!
	640	! endif
	641	! endif
	642
	643	! do l=1,llm
	644	! print,'effect interp',l,qk1d(l,2),qk1dint(l,2)facQ(l,2),
	645	! & TempI1d(l,2),TempI1din(l,2)*facTI(l,2),TempE1D(l),
	646	! & TempE1Dint(l)*facTE(l),facQ(l,:),FacTI(l,:),FacTe(l)
	647	! enddo
	648	! stop
	649
	650	! Compute vertical velocity to derive subtime step
	651
	652	Tdiff=1.E-9
	653	DO k=1,ndiffions
	654	kn=itrans(k)
	655	H0=kbolt*T0/mol_tr(kn)/g/masseU
	656	D0=DIraf(naltvert,k)
	657	Time0=H0*H0/D0
	658	Time=Tdiff/Time0
	659	do ln=1,naltvert
	660	RIad(ln)=RIraf(ln,k)/Rho0
	661	RElad(ln)=Reraf(ln)/rho0
	662	TIad(ln)=TIraf(ln,k)/T0
	663	TEad(ln)=TEraf(ln)/T0
	664	Zad(ln)=Zraf(ln)/H0
	665	Dad(ln)=DIraf(ln,k)/D0
	666	enddo
	667
	668	DZ=Zad(2)-Zad(1)
	669
	670	CALL DIFFPARAMK(Tiad,Tead,RElad,Dad,DZ,
	671	& alpha,beta,gama,delta,eps,dzeta,eta,naltvert,Wmax)
	672
	673	! if (ig .eq. ij0) then
	674	! print*,'zeta',dzeta(:)
	675	! print*,'eta',eta(:)
	676	! endif
	677
	678	FacEsc=exp(-ueff(ig,k)Time0/H0DZ/Dad(naltvert-1))
	679
	680	CALL MATCOEFFK(alpha,beta,gama,delta,eps,dzeta,eta,Dad,Riad,
	681	& FacEsc,dZ,time,Atri,Btri,Ctri,Dtri,Tiad,Tead,naltvert)
	682
	683	rhock(1,k)=0D0
	684	rhock(2,k)=0D0
	685	do ln=3,naltvert-1
	686
	687	if (ln .lt. naltvert-1) then
	688	rhock(ln,k)=-Dad(ln)(-RIraf(ln+2,k)+8.RIraf(ln+1,k)
	689	& -8.*RIraf(ln-1,k)+RIraf(ln-2,k))/RIraf(ln,k)/12D0/DZ
	690	& -Dad(ln)*(alpha(ln)+beta(ln)+dzeta(ln))
	691	endif
	692
	693	if (ln .eq. naltvert-1) then
	694	rhock(ln,k)=-Dad(ln)*
	695	& (log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz
	696	& -Dad(ln)*(alpha(ln)+beta(ln)+dzeta(ln))
	697	endif
	698
	699	rhock(ln,k)=D0/H0*rhock(ln,k)
	700	enddo
	701
	702	rhock(naltvert,k)=ueff(ig,k)
	703	enddo
	704
	705
	706	Wmax=maxval(abs(rhock))
	707	if (mu_local .ge. 0.3) tdiff3=h0*Dz/Wmax/300d0 !20d0 !/100d0
	708	if (mu_local .le. -0.3) tdiff3=h0*dz/Wmax/300d0 !300d0
	709	if (abs(mu_local) .lt. 0.3) tdiff3=h0*dz/Wmax/400d0 !/1000d0
	710
	711	tmin=max(tdiff1,tdiff3)
	712	tmin=min(tmin,tdiff2)
	713
	714	nsubreal=anint(dt/tmin)
	715
	716	tdiff=tmin
	717
	718	!! This is a weird factor to optimize time calculation
	719	!! Here, it is to be sure that the first dtime is very low
	720	tmin = 20.45(DZH0)*2.
	721	tmin = tmin/(maxval(abs(DIraf(:,:)))+(DZH0)Wmax)
	722	if (tdiff .ge. tmin) then
	723	tdiff = tmin
	724	endif
	725	! if (ig .eq. ij0) then
	726	! print*,'nreal',nsubreal,nsubvert,nsubmin,
	727	! & Wmax,tmin,tdiff,tdiff1,tdiff2,tdiff3,h0*dZ,tsim,ttot,h0,dZ
	728	! write(301,*) nsubreal
	729	! endif
	730
	731	!tdiff=8.75D-5 !tmin
	732	tdiff=tmin
	733
	734	! Fix a minimal value for tdiff
	735
	736	! ========== TEMPORAL Loop =================
	737
	738	! Loop over dynamic ions to solve the diffusion equation
	739
	740	! Compute total mass on the Z grid before diffusion
	741
	742	! MtotZ1(:)=0D0
	743	! do k=1,ndiffions
	744	! do l=1,naltvert
	745	! MtotZ1(k)=MtotZ1(k)+RIraf(l,k)*Dzraf
	746	! enddo
	747	! enddo
	748	! print*,'MtotZ',MtotZ1
	749
	750	!if (szad_local.ge.93.) tsim = ttot
	751	! DO istep=1,nsubreal
	752	DO WHILE (tsim .lt. ttot)
	753	! if (ig .eq. ij0) then
	754	! print*,'tsim',tsim
	755	! do l=1,naltvert
	756	! print*,'density1',l,REraf(l)/Mtraceur(31)/masseU/1d6,
	757	! & Riraf(l,2)/Mtraceur(22)/masseU/1d6
	758	! enddo
	759	! endif
	760	RErafold=REraf
	761	RIrafold=RIraf
	762	DO iter=1,niter ! needed to use approximatively electron density (t+dt)
	763	DO ln=1,naltvert
	764	RElad(ln)=REraf(ln)/Rho0
	765	REraf(ln)=RErafold(ln)
	766
	767	! if (ig .eq. ij0) then
	768	! print,'density',RElad(ln)/Mtraceur(31)/masseU/1d6rho0
	769	! & ,Riraf(ln,2)/Mtraceur(22)/masseU/1d6
	770	! endif
	771
	772	ENDDO ! ln
	773
	774	DO k=1,ndiffions
	775	kn=itrans(k)
	776	H0=kbolt*T0/mol_tr(kn)/g/masseU
	777	D0=DIraf(naltvert,k)
	778	Time0=H0*H0/D0
	779	Time=Tdiff/Time0
	780
	781	! Compute the diffusion coefficient using the collion frequency
	782
	783	! print*,'now adimension'
	784
	785	! Adimension all parameters by value at iz_vertplasma+1
	786
	787	do ln=1,naltvert
	788	RIad(ln)=RIraf(ln,k)/Rho0
	789	! RElad(ln)=REraf(ln)/Rho0
	790	TIad(ln)=TIraf(ln,k)/T0
	791	TEad(ln)=TEraf(ln)/T0
	792	Zad(ln)=Zraf(ln)/H0
	793	Dad(ln)=DIraf(ln,k)/D0
	794	enddo ! ln
	795
	796	DZ=Zad(2)-Zad(1)
	797
	798	! Tead(:)=0D0
	799
	800	CALL DIFFPARAMK(Tiad,Tead,RElad,Dad,DZ,
	801	& alpha,beta,gama,delta,eps,dzeta,eta,naltvert,Wmax)
	802
	803	! if (ig .eq. ij0) then
	804	! print*,'zeta',dzeta(:)
	805	! print*,'eta',eta(:)
	806	! endif
	807
	808	FacEsc=exp(-ueff(ig,k)Time0/H0DZ/Dad(naltvert-1))
	809	! & *Tiad(ln-1)/(Tiad(ln-1)+Tead(ln-1)))
	810
	811	! eta(:)=0D0
	812	! dzeta(:)=0D0
	813
	814	! do l=1,naltvert
	815	! print*,l,alpha(l),dzeta(l),eta(l),RElad(l),Riad(l)
	816	! enddo
	817	! stop
	818
	819
	820	CALL MATCOEFFK(alpha,beta,gama,delta,eps,dzeta,eta,Dad,
	821	& Riad,FacEsc,dZ,time,Atri,Btri,Ctri,Dtri,Tiad,Tead,naltvert)
	822
	823	Xtri(:)=0D0
	824
	825	Xtri=Riad
	826	Ytri=Xtri
	827
	828	! Solve system
	829	!write(,) tsim, iter
	830	CALL TridagDP(Atri,Btri,Ctri,Dtri,Xtri,naltvert)
	831
	832	!STOP
	833	! Check mass conservation
	834
	835	! CALL CheckmassK(RIad,Xtri,naltvert)
	836
	837	! if (ig .eq. 1653) then
	838	! do l=1,naltvert
	839	! print,'lphysnew',l,rho0Xtri(l),k
	840	! enddo
	841	! endif
	842
	843	do l=1,naltvert
	844	if (iter .eq. niter) RIraf(l,k)=Rho0*Xtri(l)
	845	if (RIraf(l,k) .lt. 0. .or. isnan(RIraf(l,k))) then
	846	! print*,'phys q <0',ig,l,k,iter,istep,Xtri(l),Ytri(l),Riraf(l,k),
	847	! & REraf(l),Tead(l)T0,Tiad(l)T0
	848	RIraf(l,k)=1D-24
	849	! RIraf(l,k)=RIraf(l-1,k)exp(-DZ(alpha(l-1)+beta(l-1)+dzeta(l-1)))
	850	! nsubreal=nsubreal+50
	851	! if (nsubreal .gt. 2000) then
	852	! print*,'oula nsubreal',ig,nsubreal
	853	! endif
	854	! goto 1
	855	endif
	856	enddo ! l
	857
	858	! Update electron mass density for next iteration
	859
	860	DO l=1,naltvert
	861	REraf(l)=REraf(l)+0.5Charge_tr(kn)mol_tr(ke)/mol_tr(kn)
	862	& * rho0*(Xtri(l)-Ytri(l))
	863	ENDDO ! l
	864
	865	! if (ig .eq. ij0 .and. k .eq. 2) then
	866	! do l=1,naltvert
	867	! print*,'l1',l,Reraf(l)/mol_tr(ke)/masseU/1d6,
	868	! & Xtri(l)*rho0/mol_tr(kn)/masseU/1d6,
	869	! & Ytri(l)*rho0/mol_tr(kn)/masseU/1d6
	870	! enddo ! l
	871	! endif
	872
	873	! Compute vertical velocity in atmospheric frame
	874
	875	IF (iter .eq. niter) THEN
	876	rhock(1,k)=0D0
	877	rhock(2,k)=0D0
	878	do ln=3,naltvert-1
	879	! rhock(ln,k)=-Dad(ln)*((RIraf(ln+1,k)-RIraf(ln-1,k))/2D0/Dz
	880	! & +RIraf(ln,k)*(alpha(ln)+beta(ln)+dzeta(ln)))
	881	! rhock(ln,k)=D0/H0*rhock(ln,k)/RIraf(ln,k)
	882
	883
	884	! rhock(ln,k)=-Dad(ln)*
	885	! & ((log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz
	886	! & +(alpha(ln)+beta(ln))*(Tiad(ln)/(Tead(ln)+Tiad(ln))))
	887	! rhock(ln,k)=D0/H0*rhock(ln,k)
	888
	889	if (ln .lt. naltvert-1) then
	890	! rhock(ln,k)=-Dad(ln)(-log(RIraf(ln+2,k))+8.log(RIraf(ln+1,k))
	891	! & -8.*log(RIraf(ln-1,k))+log(RIraf(ln-2,k)))
	892	! & /12D0/Dz
	893	rhock(ln,k)=-Dad(ln)(-RIraf(ln+2,k)+8.RIraf(ln+1,k)
	894	& -8.*RIraf(ln-1,k)+RIraf(ln-2,k))/RIraf(ln,k)/12D0/DZ
	895	& -Dad(ln)*(alpha(ln)+beta(ln)+dzeta(ln))
	896	endif
	897
	898	if (ln .eq. naltvert-1) then
	899	! rhock(ln,k)=-Dad(ln)((Tiad(ln)+Tead(ln))/Tiad(ln)
	900	! & (log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz)
	901	! & -Dad(ln)*(alpha(ln)+beta(ln))
	902
	903	rhock(ln,k)=-Dad(ln)*
	904	& (log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz
	905	& -Dad(ln)*(alpha(ln)+beta(ln)+dzeta(ln))
	906
	907	endif
	908
	909	rhock(ln,k)=D0/H0*rhock(ln,k)
	910
	911	! if (RIraf(ln,k) .eq. 1D-24 .or. RIraf(ln+1,k) .eq. 1D-24)
	912	! & rhock(ln,k)=0.
	913	! I limit the vertical velocity to 1km/s to avoid instabilities in velocity advection (TO BE IMPROVED)
	914
	915	! print*,'rhock',istep,ln,k,Zraf(ln),ln,rhock(ln,k),RIraf(ln,k),
	916	! & Facesc
	917	! & RIraf(ln,k)exp(-DZ(alpha(ln)+beta(ln)+dzeta(ln))),
	918	! & (alpha(ln)+beta(ln)+dzeta(ln)),
	919	! & (log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz,Atri(ln+1),
	920	! & exp(-DZ*(alpha(ln)+beta(ln)+dzeta(ln)))
	921
	922	! if (abs(rhock(ln,k)) .gt. 1000. .or. isnan(rhock(ln,k))) then
	923	! print*,ig,istep
	924	! print*,'rhock',Zraf(ln),ln,rhock(ln,k),RIraf(ln+1,k)
	925	! & RIraf(ln,k)exp(-DZ(alpha(ln)+beta(ln)+dzeta(ln))),
	926	! & (alpha(ln)+beta(ln)+dzeta(ln)),
	927	! & (log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz,Atri(ln+1),
	928	! & exp(-DZ*(alpha(ln)+beta(ln)+dzeta(ln)))
	929	! if (abs(rhock(ln,k)) .gt. 10000.) stop
	930	! rhock(ln,k) =1000.*rhock(ln,k)/abs(rhock(ln,k))
	931	! endif
	932
	933	enddo ! ln
	934	rhock(naltvert,k)=ueff(ig,k)
	935	ENDIF ! (iter .eq. niter)
	936	ENDDO ! dynamic ions
	937
	938	! Upadte the electron density to assume neutral equilibrium
	939
	940	IF (iter .eq. niter) THEN
	941	RERaf=Rerafold
	942	!ke=etrans(1)
	943	DO l=1,naltvert
	944	DO k=1,ndiffions
	945	kn=itrans(k)
	946	REraf(l)=REraf(l)+Charge_tr(kn)*mol_tr(ke)/mol_tr(kn)
	947	& * (RIraf(l,k)-Rirafold(l,k))
	948	ENDDO
	949	! if (ig .eq. ij0) then
	950	! DO k=1,ndiffions
	951	! kn=itrans(k)
	952	! print*,'fin cycle',l,Rerafold(l)/Mtraceur(ke)/masseU/1d6,
	953	! & Reraf(l)/Mtraceur(ke)/masseU/1d6,
	954	! & Rirafold(l,k)/Mtraceur(kn)/masseU/1d6,
	955	! & Riraf(l,k)/Mtraceur(kn)/masseU/1d6
	956	! enddo
	957	! endif
	958	ENDDO
	959	ENDIF ! (iter .eq. niter)
	960
	961	! Compute vertical component of electric field
	962
	963	! DZ=Zraf(2)-Zraf(1)
	964	Ez(1)=0D0
	965	DO ln=2,naltvert-1
	966	Ez(ln)=(log(Teraf(ln+1))-log(Teraf(ln-1)))+
	967	& (log(Reraf(ln+1))-log(Reraf(ln-1)))
	968	Ez(ln)=-kbolt/qchargeTeraf(ln)/2D0/Dz/H0Ez(ln)
	969	ENDDO
	970	Ez(naltvert)=3D0log(Reraf(naltvert))-4D0log(Reraf(naltvert-1))
	971	& +log(Reraf(naltvert-2))+
	972	& 3D0log(Teraf(naltvert))-4D0log(Teraf(naltvert-1))
	973	& +log(Teraf(naltvert-2))
	974	Ez(naltvert)=-kbolt/qchargeTeraf(naltvert)/2D0/Dz/H0Ez(naltvert)
	975	! Ez(naltvert)=0D0
	976
	977	! if (ig .eq. ij0) then
	978	! do ln=1,naltvert
	979	! ! print*,'reraf write',ln,Reraf(ln)/Mtraceur(31)/masseU/1d6
	980	! write(301,*) ig,istep,ln,Zraf(ln),RIraf(ln,:),rhock(ln,:),
	981	! & REraf(ln),TIraf(ln,:),Teraf(ln),Tnraf(ln)
	982	! enddo
	983	! endif
	984
	985	ENDDO ! iteration
	986
	987	tsim=tsim+tdiff
	988	! compute new time step
	989	Wmax=maxval(abs(rhock))
	990	! if (mu_local .ge. 0.) tdiff3=Dz/Wmax/300d0
	991	! if (mu_local .lt. 0.) tdiff3=dz/Wmax/500d0
	992	!if (mu_local .ge. 0.5) tdiff3=h0*Dz/Wmax/20d0!20d0 !20d0 !/100d0
	993	!if (mu_local .le. -0.5) tdiff3=2d016d0h0*dz/Wmax/20d0!20d0 !300d0
	994	!if (abs(mu_local) .lt. 0.5) tdiff3=2d016d0h0*dz/Wmax/20d0!20d0 !400d0 !/1000d0
	995
	996	!! This is a weird factor to optimize time calculation
	997	if (mu_local .ge. 0.3) tdiff3=2d016d0h0*Dz/Wmax/20d0!20d0 !20d0 !/100d0
	998	if (mu_local .le. -0.3) tdiff3=2d016d0h0*dz/Wmax/20d0!20d0 !300d0
	999	if (abs(mu_local) .lt. 0.3) tdiff3=2d016d0h0*dz/Wmax/20d0!20d0 !400d0 !/1000d0
	1000
	1001	tmin=max(tdiff1,tdiff3)
	1002	tmin=min(tmin,tdiff2)
	1003	tdiff=tmin
	1004	!! This is a weird factor to optimize time calculation
	1005	!!tmin = 20.45(DZH0)*2./maxval(abs(DIraf(:,:)))
	1006	tmin = 3.8.4.20020.45(DZH0)*2.
	1007	tmin = tmin/(maxval(abs(DIraf(:,:)))+(DZH0)Wmax)
	1008	if (tdiff .ge. tmin) then
	1009	tdiff = tmin
	1010	endif
	1011	tdiff=max(tdiff1,tdiff)
	1012	!if (szad_local .le. 20.) tdiff = tdiff/20.
	1013	!!if (tsim .le. 8.75D-2) tdiff = 8.75D-5
	1014	if (tsim+tdiff .gt. ttot) tdiff=ttot-tsim
	1015
	1016	! if (ig .eq. ij0) then
	1017	! print*,'dt',Wmax,tdiff1,tdiff2,tdiff3,tdiff,tsim,ttot,
	1018	! & alpha(naltvert-1),beta(naltvert-1),dzeta(naltvert-1),
	1019	! & Atri(naltvert),Xtri(naltvert),Xtri(naltvert-1)*Atri(naltvert),
	1020	! & h0,dz,h0*dz,mu_local
	1021	! endif
	1022
	1023	ENDDO ! time step ou while
	1024	! nsubreal=nsubvert
	1025	! end time step here
	1026	! if (ig .eq. ij0) close(301)
	1027
	1028	MtotZ2(:)=0D0
	1029
	1030	! Limites pour la vitesse verticale pour eviter problemes numeriques
	1031
	1032	DO k=1,ndiffions
	1033	DO ln=1,naltvert
	1034	IF (abs(rhock(ln,k)).gt.Wlim .or. isnan(rhock(ln,k))) THEN
	1035	rhock(ln,k) = Wlim*rhock(ln,k)/abs(rhock(ln,k))
	1036	ENDIF
	1037	!IF (k.eq.1) THEN
	1038	! write(,) real(qk1Dnew(ln,:)), qold(ig,ln,gcmind(ke))
	1039	!ENDIF
	1040	ENDDO ! ln
	1041	ENDDO ! k
	1042
	1043	! print*,'rhock',l,k,Zraf(l),l,rhock(l,k),RIraf(l,k),Facesc
	1044	! MtotZ2(k)=MtotZ2(k)+RIraf(l,k)*Dzraf
	1045	! enddo
	1046	! enddo
	1047
	1048	! We reinterpolate results on the Pressure grid and correct with FacQ
	1049
	1050	CALL GCMGRID_P2K(Zraf,Praf,TNraf,Mraf,QIraf,RIraf,TIraf,TEraf,
	1051	& rhock,Ez,P1D,Pnlay1D,TempN1d,Mmean1d,qk1d,qk1dnew,
	1052	& TempI1D,TempI1Dnw,TempE1D,TempE1dnew,Wk1D,Wk1D2,Wk1d3,
	1053	& Ez1d,ZZ1d,naltvert,ndiffions,nlayer,
	1054	& iz_vertplasma,itrans,etrans,FacQ,FacTI,FacTe)
	1055
	1056	!iq = itrans(1)
	1057	!write(,) tname(gcmind(iq)),ptimestep,qk1dnew(:,1)
	1058
	1059
	1060	! do l=1,llm+1
	1061	! print*,'Ez,W',P1D(l),Ez1D(l),Wk1D(l,:)
	1062	! enddo
	1063
	1064	! Compute the mixing ratio and temperature & Update potential temperature
	1065	! Here I neglect advection of temperature (later)
	1066
	1067	! do ln=1,llm
	1068	! l=ln-iz_plasma
	1069	! Alt(ig,ln)=ZZ1d(ln)
	1070	! do k=1,ndiffions
	1071	! kn=itrans(k)
	1072	!! if (ig .eq. ij0) print*,'qold',ln,q(ig,ln,kn),qk1d(ln,k)
	1073	! q(ig,ln,kn)=real(qk1Dnew(ln,k))
	1074	!! if (ig .eq. ij0) print*,'qnew',q(ig,ln,kn),qk1dnew(ln,k)
	1075	! if (l .gt. 0) then
	1076	! Tetak(ig,l,k)=real(tempI1Dnw(ln,k))*(1D0-kappak(k))
	1077	! & (masseUMtraceur(kn)/kboltPreffplasma/RhoN(ig,ln)/q(ig,ln,kn))
	1078	! & **(kappak(k))
	1079	! wcontk(ig,l,k)=real(Wk1D(ln,k))
	1080	! wcovk(ig,l,k)=real(Wk1D2(ln,k))
	1081	! wphysk(ig,l,k)=real(Wk1D3(ln,k))
	1082	! EfieldZ(ig,l)=real(Ez1D(ln))
	1083	! endif
	1084	!
	1085	! if (q(ig,ln,kn) .lt. 0. .or. isnan(q(ig,ln,kn))) then
	1086	! print*,'aie q < 0',q(ig,ln,kn),ln,l,ig,k,kn
	1087	! print*,qk1dnew(:,k),qk1d(:,k)
	1088	! stop
	1089	! endif
	1090	!
	1091	! if (l .ge. 1) then
	1092	! if (isnan(Wcontk(ig,l,k))) then
	1093	! Wcontk(ig,l,k)=0.
	1094	! Wcovk(ig,l,k)=0.
	1095	! ! print*,'aie Wcontk',l,ln,ig,Wcontk(ig,l,k),Wk1D(ln,k)
	1096	! ! do l=1,naltvert
	1097	! ! print*,'rhock',l,rhock(l,:),RIraf(l,:),Dad(l),alpha(l),beta(l),
	1098	! ! & dzeta(l),Praf(l)
	1099	! ! enddo
	1100	! endif
	1101	! if (isnan(Efieldz(ig,l))) then
	1102	! Efieldz(ig,l)=0.
	1103	! endif
	1104	! endif
	1105	!
	1106	! enddo
	1107	! enddo
	1108	! do k=1,ndiffions
	1109	! Wcontk(ig,llp+1,k)=real(Wk1D(llm+1,k))
	1110	! Wcovk(ig,llp+1,k)=real(Wk1D2(llm+1,k))
	1111	! wphysk(ig,llp+1,k)=real(Wk1D3(llm+1,k))
	1112	! enddo
	1113	! EfieldZ(ig,llp+1)=Ez1D(llm+1)
	1114
	1115	DO ln=1,nlayer
	1116	l=ln-iz_vertplasma
	1117	DO k=1,ndiffions
	1118	iq=gcmind(itrans(k))
	1119	qnew(ig,ln,iq)=real(qk1Dnew(ln,k))
	1120	ENDDO
	1121	ENDDO
	1122
	1123	ENDDO ! END ig - Main Loop on horizontal grids - plan horizontal
	1124
	1125	! ==========================================================
	1126	! ===== Compute the diffusion trends du to diffusion =======
	1127	! ================== for the ions ==========================
	1128	! ==========================================================
	1129
	1130	!ke = etrans(1)
	1131	DO l=1,nlayer
	1132	DO k=1,ndiffions
	1133	iq = gcmind(itrans(k))
	1134	pdqiondiff(:,l,iq) = (qnew(:,l,iq)-qold(:,l,iq))/ptimestep
	1135	! CE CALCUL EST FAUX, CE N'EST PAS UNE CONSERVATION MMR QU'IL FAUT MAIS UNE CONSERVATION VMR !!!!!!!!!!
	1136	! pdqiondiff(:,l,iqelec)= pdqiondiff(:,l,iqelec)
	1137	! & + pdqiondiff(:,l,iq)
	1138	!pdqiondiff(:,l,iq) =(qnew(:,l,kn)-qold(:,l,kn))/ptimestep
	1139	!pdqiondiff(:,l,iqelec)= pdqiondiff(:,l,ke)+pdqiondiff(:,l,iq)
	1140	ENDDO
	1141	ENDDO
	1142	!iq = gcmind(itrans(1))
	1143	!write(,) tname(iq),iq,ptimestep,pdqiondiff(1,:,iq),qold(1,:,iq)
	1144	!write(,) maxval(pdqiondiff(:,:,iq))
	1145	! ==========================================================
	1146	! ===== Compute the diffusion trends du to diffusion =======
	1147	! ================== for the electron ======================
	1148	! ==========================================================
	1149
	1150	! ------ update mmean ------ !
	1151	mmean_new(:,:) = 0.
	1152	do iq = 1,ntracers
	1153	mmean_new(:,:) = mmean_new(:,:) +
	1154	& qnew(:,:,gcmind(iq))/M_tr(gcmind(iq))
	1155	enddo
	1156	mmean_new(:,:) = 1./mmean_new(:,:)
	1157
	1158	! ------ conversion mmr into vmr ------ !
	1159	do iq = 1,ntracers
	1160	vmr_new(:,:,gcmind(iq)) = qnew(:,:,gcmind(iq)) *
	1161	& mmean_new(:,:)/M_tr(gcmind(iq))
	1162	enddo
	1163
	1164	iqelec = gcmind(ke)
	1165	! ------ vmr ion density ------ !
	1166	vmr_ion(:,:) = 0.
	1167	do iq = 1,ntracers
	1168	if ((type_tr(gcmind(iq)) .eq. 2) .and.
	1169	& ( tname(gcmind(iq)) .ne. tname(iqelec))) then
	1170	vmr_ion(:,:) = vmr_ion(:,:) + vmr_new(:,:,gcmind(iq))
	1171	endif
	1172	enddo
	1173
	1174	! ------ force charge neutrality ------ !
	1175	! ------ vmr(ion) = vmr(elec) ------ !
	1176	DO ig=1,ngrid
	1177	DO l=1,nlayer
	1178	if (vmr_new(ig,l,iqelec) .ne. vmr_ion(ig,l)) then
	1179	vmr_new(ig,l,iqelec) = vmr_ion(ig,l)
	1180	! DO k=1,ndiffions
	1181	! vmr_new(ig,l,iqelec) = vmr_new(ig,l,iqelec) +
	1182	! & vmr_new(ig,l,gcmind(itrans(k)))
	1183	! ENDDO
	1184	endif
	1185	ENDDO
	1186	ENDDO
	1187
	1188	! ------ conversion vmr into mmr for electron ------ !
	1189	qnew(:,:,iqelec) = vmr_new(:,:,iqelec)*m_tr(iqelec)/mmean_new(:,:)
	1190
	1191	! ------ Compute the diffusion trends du to diffusion ------ !*
	1192	! ------ for the electron ------ !
	1193	DO l=1,nlayer
	1194	pdqiondiff(:,l,iqelec) =
	1195	& (qnew(:,l,iqelec)-qold(:,l,iqelec))/ptimestep
	1196	ENDDO
	1197
	1198
	1199	! write(,) 'TATA EST PLUS LA'
	1200	! do ig=1,ip1jm+1,iip1
	1201	! do l=1,llp
	1202	! ln=l+iz_plasma
	1203	! do k=1,ndiffions
	1204	! kn=itrans(k)
	1205	! q(ig+iim,ln,kn)=q(ig,ln,kn)
	1206	! tetak(ig+iim,l,k)=tetak(ig,l,k)
	1207	! wcontk(ig+iim,l,k)=wcontk(ig,l,k)
	1208	! if (ig .eq. ij0 .and. l .eq. 2 .and. k .eq.1)
	1209	! & print*,ij0+iim,wcontk(ij0,2,1),wcontk(ij0+iim,2,1)
	1210	! wcovk(ig+iim,l,k)=wcovk(ig,l,k)
	1211	! wphysk(ig+iim,l,k)=wphysk(ig,l,k)
	1212	! Efieldz(ig+iim,l)=Efieldz(ig,l)
	1213	! tetae(ig+iim,l)=tetae(ig,l)
	1214	! enddo
	1215	! enddo
	1216	! enddo
	1217	! print*,'wphys3',ij0,ij0+iim,wcontk(ij0,2,1),wcontk(ij0+iim,2,1)
	1218
	1219	RETURN
	1220	END
	1221
	1222	! ********************************************************************
	1223	! ********************************************************************
	1224	! ********************************************************************
	1225
	1226	SUBROUTINE TridagDP(a,b,c,r,u,n)
	1227	! USE mod_phys_lmdz_para, only: mpi_rank
	1228	! parameter (nmax=5000)
	1229	! dimension gam(nmax),a(n),b(n),c(n),r(n),u(n)
	1230	real(kind=kind(1.D0)) gam(n),a(n),b(n),c(n),r(n),u(n)
	1231	if(b(1).eq.0.)then
	1232	stop 'tridag: error: b(1)=0 !!! '
	1233	endif
	1234	bet=b(1)
	1235	u(1)=r(1)/bet
	1236	do 11 j=2,n
	1237	gam(j)=c(j-1)/bet
	1238	bet=b(j)-a(j)*gam(j)
	1239	if(bet.eq.0.) then
	1240	stop 'tridag: error: bet=0 !!! '
	1241	endif
	1242	u(j)=(r(j)-a(j)*u(j-1))/bet
	1243	11 continue
	1244	do 12 j=n-1,1,-1
	1245	u(j)=u(j)-gam(j+1)*u(j+1)
	1246	12 continue
	1247	return
	1248	END
	1249
	1250
	1251	! SUBROUTINE DYN2DIFFK(Pk,teta,P,mmean,rhoN,Preff,cpp,kappa,
	1252	! & q,Mtraceur,Tempk,tetak,TempE,tetaE,Preffplasma,itrans,etrans,
	1253	! & llm,nqtot,llp,ndiffions,ip1jmp1,iz_plasma,kappak,kappae,
	1254	! & Pnlay,TempN,Rhok,TempI,TetaI,rhotetaI,TempF,TetaF,rhotetaF)
	1255	!
	1256	! USE infotrac, only: tname
	1257	!
	1258	! IMPLICIT NONE
	1259	!
	1260	! INTEGER :: llm,llp,nqtot,ndiffions,ip1jmp1,l,ig,k,ln,kn,iz_plasma
	1261	! REAL :: Pk(ip1jmp1,llm),teta(ip1jmp1,llm)
	1262	! REAL :: mmean(ip1jmp1,llm),rhoN(ip1jmp1,llm)
	1263	! REAL :: P(ip1jmp1,llm+1),Pnlay(ip1jmp1,llm)
	1264	! REAL :: q(ip1jmp1,llm,nqtot)
	1265	! REAL :: Tempk(ip1jmp1,llp,ndiffions),tetak(ip1jmp1,llp,ndiffions)
	1266	! REAL :: TempE(ip1jmp1,llp),TetaE(ip1jmp1,llp)
	1267	! REAL :: Mtraceur(nqtot)
	1268	! INTEGER :: itrans(ndiffions),etrans(1)
	1269	! REAL :: Preff,Cpp,Preffplasma,kappa,unskappa
	1270	! REAL :: kappak(ndiffions)
	1271	! REAL :: kappae
	1272	!
	1273	! ! Output
	1274	! REAL :: TempN(ip1jmp1,llm)
	1275	! REAL :: Rhok(ip1jmp1,llm,nqtot)
	1276	! REAL :: TempI(ip1jmp1,llm,ndiffions)
	1277	! REAL :: TetaI(ip1jmp1,llm,ndiffions)
	1278	! REAL :: rhotetaI(ip1jmp1,llm,ndiffions)
	1279	! REAL :: TempF(ip1jmp1,llm)
	1280	! REAL :: TetaF(ip1jmp1,llm)
	1281	! REAL :: rhotetaF(ip1jmp1,llm)
	1282	!
	1283	! REAL :: masseU,kbolt
	1284	! masseU=1.660538782d-27
	1285	! kbolt=1.3806504d-23
	1286	!
	1287	! unskappa=1./kappa
	1288	!
	1289	! Rhok(:,:,:)=0.
	1290	! DO ig=1,ip1jmp1
	1291	! DO ln=1,llm
	1292	! l=ln-iz_plasma
	1293	! TempN(ig,ln)=Teta(ig,ln)*pk(ig,ln)/Cpp
	1294	! Pnlay(ig,ln)=preff(pk(ig,ln)/cpp)*unskappa
	1295	! do k=1,nqtot
	1296	! if (tname(k) .ne. "dust_number" .and. tname(k) .ne. "dust_mass"
	1297	! & .and. tname(k).ne."ccn_number".and.tname(k).ne."ccn_mass") then
	1298	! if (q(ig,ln,k) .le. 0.) q(ig,ln,k)=1d-30
	1299	! Rhok(ig,ln,k)=q(ig,ln,k)*RhoN(ig,ln)
	1300	!! if (q(ig,ln,k) .le. 0. .and. l .le. 0) Rhok(ig,ln,k)=1d-30*rhoN(ig,ln)
	1301	! endif
	1302	! enddo
	1303	!
	1304	!! Ions parameters
	1305	! do k=1,ndiffions
	1306	! kn=itrans(k)
	1307	! if (ln .le. iz_plasma) TempI(ig,ln,k)=TempN(ig,ln)
	1308	! if (ln .gt. iz_plasma) TempI(ig,ln,k)=Tempk(ig,l,k)
	1309	!
	1310	! if (ln .le. iz_plasma) TetaI(ig,ln,k)=teta(ig,ln)*
	1311	! & (Preffplasmakappak(k))/(Preffkappa)*
	1312	!! & (Mtraceur(kn)/mmean(ig,ln)/q(ig,ln,kn))**kappa
	1313	! & (rhoN(ig,ln)/mmean(ig,ln))*kappa
	1314	! & 1D0/(rhok(ig,ln,kn)/Mtraceur(kn))*kappak(k)
	1315	! & (kboltTempN(ig,ln)/masseu)*(kappa-kappak(k))
	1316	! if (ln .gt. iz_plasma) TetaI(ig,ln,k)=Tetak(ig,l,k)
	1317	!
	1318	! rhotetaI(ig,ln,k)=rhok(ig,ln,kn)*TetaI(ig,ln,k)
	1319	! enddo
	1320	!
	1321	!! Electron parameters
	1322	!
	1323	! kn=etrans(1)
	1324	! if (ln .le. iz_plasma) TempF(ig,ln)=TempN(ig,ln)
	1325	! if (ln .gt. iz_plasma) TempF(ig,ln)=TempE(ig,l)
	1326	!
	1327	! if (ln .le. iz_plasma) TetaF(ig,ln)=teta(ig,ln)*
	1328	! & (Preffplasmakappae)/(Preffkappa)*
	1329	!! & (Mtraceur(kn)/mmean(ig,l)/q(ig,ln,kn))**kappa
	1330	! & (rhoN(ig,ln)/mmean(ig,ln))*kappa
	1331	! & 1D0/(rhok(ig,ln,kn)/Mtraceur(kn))*kappae
	1332	! & (kboltTempN(ig,ln)/masseU)*(kappa-kappae)
	1333	! if (ln .gt. iz_plasma) TetaF(ig,ln)=TetaE(ig,l)
	1334	!
	1335	! rhotetaF(ig,ln)=rhok(ig,ln,kn)*TetaF(ig,ln)
	1336	!
	1337	!
	1338	! ENDDO ! l or l=ln-iz_plasma
	1339	! ENDDO ! ig
	1340	! END
	1341
	1342	SUBROUTINE ZVERTK(P,T,M,Z,nl,gsurf)
	1343	IMPLICIT NONE
	1344	#include "YOMCST.h"
	1345	integer :: nl,l
	1346	REAL(kind=kind(1.D0)):: P(nl),T(nl),M(nl),Z(nl)
	1347	REAL(kind=kind(1.D0)):: masseU,kbolt,Konst,g,z0,Hpm
	1348	REAL(kind=kind(1.D0)):: Tmean,Mmean,gsurf
	1349	masseU=1.e-3/RNAVO
	1350	kbolt=RKBOL
	1351	Konst=kbolt/masseU
	1352	g=RG
	1353
	1354	Z0=gsurf/g!0d0
	1355	Z(1)=z0
	1356	Hpm=Konst*T(1)/g/M(1)
	1357
	1358	do l=2,nl
	1359	Tmean=T(l-1)
	1360	Mmean=M(l-1)
	1361	Hpm=Konst*Tmean/g/Mmean
	1362	Z(l)=Z(l-1)-Hpm*log(P(l)/P(l-1))
	1363	! print*,'z',l,Z(l)
	1364	enddo
	1365
	1366	END
	1367
	1368	SUBROUTINE UPPER_RESOLK(ZZ1D,P1D,TN1d,M1d,rhok1d,TI1d,Te1d,
	1369	& mol_tr,nlayer,ndiffions,ntracers,nlraf,idiffZ,
	1370	& itrans,etrans,Zraf,TNraf,Praf,Mraf,Qraf,rhokraf,
	1371	& QIraf,TIraf,FIraf,QEraf,Teraf,IndZ,gcmind)
	1372	USE infotrac_phy, only: tname
	1373	IMPLICIT NONE
	1374	#include "YOMCST.h"
	1375	INTEGER :: nlayer,ndiffions,ntracers,nlraf,idiffZ,iq,k,kn,ke,iz,l
	1376	INTEGER :: indZ(1)
	1377	INTEGER :: itrans(ndiffions),etrans(1),gcmind(ntracers)
	1378	REAL :: mol_tr(ntracers)
	1379	REAL(kind=kind(1.D0)),dimension(nlayer) :: P1D,TN1D,M1D,ZZ1D,Te1D
	1380	REAL(kind=kind(1.D0)),dimension(nlayer,ntracers) :: rhok1d
	1381	REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: TI1d
	1382	REAL(kind=kind(1.D0)),dimension(nlraf) :: Zraf,Praf,TNraf
	1383	REAL(kind=kind(1.D0)),dimension(nlraf) :: Mraf
	1384	REAL(kind=kind(1.D0)),dimension(nlraf,ntracers) :: Qraf,Rhokraf
	1385	REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: QIraf,FIraf
	1386	REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: TIraf
	1387	REAL(kind=kind(1.D0)),dimension(nlraf) :: Qeraf,TEraf
	1388	REAL(kind=kind(1.D0)) :: kbolt,masseU,Konst,g,a0,apol
	1389	REAL(kind=kind(1.D0)) :: mu,nlay,facZ,dZ
	1390	masseU=1.e-3/RNAVO
	1391	kbolt=RKBOL
	1392	Konst=kbolt/masseU
	1393	g=RG
	1394
	1395	a0=2.9D-15 ! absolute constant (when densities expressed in /m3)
	1396	!apol=2.911D0 ! in 10^-24 cm3 (Withers 2009 for CO2)
	1397
	1398	Zraf(1) = ZZ1D(idiffZ)
	1399	Praf(1) = P1D(idiffZ)
	1400	TNraf(1) = TN1D(idiffZ)
	1401	Qraf(:,:) = 0d0
	1402	Rhokraf(1:nlraf,1:ntracers)=0d0
	1403	do iq=1,ntracers
	1404	! if (tname(iq) .ne. "dust_number" .and. tname(iq) .ne. "dust_mass"
	1405	! & .and. tname(iq).ne."ccn_number".and.tname(iq).ne."ccn_mass") then
	1406	Rhokraf(1,iq)=rhok1d(idiffZ,iq)
	1407	Qraf(1,iq)=rhok1d(idiffZ,iq)/sum(rhok1d(idiffZ,1:ntracers))
	1408	! endif
	1409	enddo
	1410	Mraf(1)=1D0/sum(Qraf(1,1:ntracers)/mol_tr(1:ntracers))
	1411	nlay=Praf(1)/kbolt/TNraf(1)
	1412	do k=1,ndiffions
	1413	kn=itrans(k)
	1414	TIraf(1,k)=TI1D(idiffZ,k)
	1415	QIraf(1,k)=Qraf(1,kn)
	1416	! mu=mol_tr(kn)*Mraf(1)/(mol_tr(kn)+Mraf(1))
	1417	! FIraf(1,k)=nlaya0Mraf(1)/(Mraf(1)+mol_tr(kn))*sqrt(apol/mu)
	1418	enddo
	1419	ke=etrans(1)
	1420	Teraf(1)=TE1D(idiffZ)
	1421	Qeraf(1)=Qraf(1,ke)
	1422
	1423	Zraf(nlraf)=ZZ1d(nlayer)
	1424	do l=2,nlraf-1
	1425	Zraf(l)=Zraf(1)+(Zraf(nlraf)-Zraf(1))/DBLE(nlraf-1)*DBLE(l-1)
	1426
	1427	iZ=1
	1428	do while (ZZ1D(iz) .le. Zraf(l))
	1429	iZ=iZ+1
	1430	enddo
	1431	Iz=iz-1
	1432
	1433	! indZ=maxloc(ZZ1D,MASK=ZZ1d < Zraf(l))
	1434	! print*,'indZ',indZ
	1435	! iZ=indZ(1)
	1436	dZ=Zraf(l)-Zraf(l-1)
	1437	facZ=(Zraf(l)-zz1d(iz))/(zz1d(iZ+1)-ZZ1d(iz))
	1438	TNraf(l)=TN1D(iz)(TN1D(iz+1)/TN1d(iz))*facZ
	1439	do iq=1,ntracers
	1440	if (Rhok1d(iz,iq) .gt. 0.) then
	1441	Rhokraf(l,iq)=Rhok1d(iz,iq) *
	1442	& (Rhok1d(iz+1,iq)/Rhok1d(iz,iq))**facZ
	1443	endif
	1444	enddo
	1445	do iq=1,ntracers
	1446	Qraf(l,iq)=Rhokraf(l,iq)/sum(Rhokraf(l,1:ntracers))
	1447	enddo
	1448	Mraf(l)=1D0/sum(Qraf(l,1:ntracers)/mol_tr(1:ntracers))
	1449	nlay=sum(Rhokraf(l,1:ntracers))/Mraf(l)/masseU
	1450	Praf(l)=nlaykboltTNraf(l)
	1451	do k=1,ndiffions
	1452	kn=itrans(k)
	1453	TIraf(l,k)=TI1D(iz,k)(TI1D(iz+1,k)/TI1D(iz,k))*facZ
	1454	Qiraf(l,k)=Qraf(l,kn)
	1455	! mu=mol_tr(kn)*Mraf(l)/(mol_tr(kn)+Mraf(l))
	1456	! FIraf(l,k)=nlaya0Mraf(l)/(Mraf(l)+mol_tr(kn))*sqrt(apol/mu)
	1457	enddo
	1458	Teraf(l)=TE1D(iz)(TE1D(iz+1)/TE1D(iz))*facZ
	1459	Qeraf(l)=Qraf(l,ke)
	1460	enddo
	1461
	1462	Zraf(nlraf)=ZZ1d(nlayer)
	1463	TNraf(nlraf)=TN1D(nlayer)
	1464	! Qraf(:,:)=0d0
	1465	do iq=1,ntracers
	1466	Rhokraf(nlraf,iq)=rhok1d(nlayer,iq)
	1467	Qraf(nlraf,iq)=rhok1d(nlayer,iq)/sum(rhok1d(nlayer,1:ntracers))
	1468	enddo
	1469	Mraf(nlraf)=1D0/sum(Qraf(nlraf,1:ntracers)/mol_tr(1:ntracers))
	1470	nlay=sum(rhokraf(nlraf,1:ntracers))/Mraf(nlraf)/masseU
	1471	Praf(nlraf)=P1D(nlayer)
	1472	do k=1,ndiffions
	1473	kn=itrans(k)
	1474	TIraf(nlraf,k)=TI1D(nlayer,k)
	1475	QIraf(nlraf,k)=Qraf(nlraf,kn)
	1476	! mu=mol_tr(kn)*Mraf(nlraf)/(mol_tr(kn)+Mraf(nlraf))
	1477	! nlay=Praf(nlraf)/kbolt/TNraf(nlraf)
	1478	! FIraf(nlraf,k)=nlaya0Mraf(nlraf)/(Mraf(nlraf)+mol_tr(kn))
	1479	! & *sqrt(apol/mu)
	1480	enddo
	1481	Teraf(nlraf)=TE1D(nlayer)
	1482	QEraf(nlraf)=Qraf(nlraf,ke)
	1483
	1484	apol=0D0
	1485	do iq=1,ntracers
	1486
	1487	! REF: CRC Handbook CHEMISTRY and Physics - 95th EDITION 2014-2015
	1488	if (tname(gcmind(iq)).eq.'co2') apol=2.911d0 ! in 10^-24 cm3 (Withers 2009 for CO2)
	1489	if (tname(gcmind(iq)).eq.'co') apol=1.95d0 ! in 10^-24 cm3
	1490	if (tname(gcmind(iq)).eq.'o2') apol=1.5689d0 ! in 10^-24 cm3
	1491	if (tname(gcmind(iq)).eq.'n2') apol=1.7403d0 ! in 10^-24 cm3
	1492	if (tname(gcmind(iq)).eq.'h2') apol=0.8042d0 ! in 10^-24 cm3
	1493	if (tname(gcmind(iq)).eq.'he') apol=0.20550522d0 ! in 10^-24 cm3
	1494	if (tname(gcmind(iq)).eq.'h') apol=0.666793d0 ! in 10^-24 cm3
	1495	if (tname(gcmind(iq)).eq.'n') apol=1.10d0 ! in 10^-24 cm3
	1496	if (tname(gcmind(iq)).eq.'o') apol=0.802d0 ! in 10^-24 cm3
	1497
	1498	if(apol.ne.0d0) then
	1499
	1500	do l=1,nlraf
	1501	nlay=Rhokraf(l,iq)/Mraf(l)/masseU
	1502	do k=1,ndiffions
	1503	kn=itrans(k)
	1504	mu=mol_tr(kn)*mol_tr(iq)/(mol_tr(kn)+mol_tr(iq))
	1505	FIraf(l,k) = FIraf(l,k) + nlay * a0
	1506	& * mol_tr(iq)/(mol_tr(iq)+mol_tr(kn)) * sqrt(apol/mu)
	1507	enddo
	1508	enddo
	1509
	1510	endif
	1511	apol=0d0
	1512
	1513	enddo
	1514
	1515	END
	1516
	1517	SUBROUTINE GCMGRID_PK(Zraf,Praf,Traf,Mraf,QIraf,RIraf,
	1518	& TIraf,TEraf,P1di,P1D,Tn1d,M1d,qk1d,qk1dnew,TI1D,TI1Dnew,
	1519	& TE1D,TE1Dnew,nlraf,ndiffions,nlayer,izv,itrans,etrans)
	1520
	1521	IMPLICIT NONE
	1522	#include "YOMCST.h"
	1523	INTEGER :: nlraf,ndiffions,nlayer,l,iP,k,izv
	1524	INTEGER,dimension(ndiffions) :: itrans
	1525	INTEGER,dimension(1) :: etrans
	1526	REAL(kind=kind(1.D0)),dimension(nlayer) :: P1d,Tn1d,M1d
	1527	REAL(kind=kind(1.D0)),dimension(nlayer+1) :: P1di
	1528	REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: Qk1d,Qk1dnew
	1529	REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: TI1d,TI1dnew
	1530	REAL(kind=kind(1.D0)),dimension(nlayer) :: TE1d,TE1Dnew
	1531	REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: Rknew
	1532	REAL(kind=kind(1.D0)),dimension(nlraf) :: Zraf,Praf,Traf,Mraf
	1533	REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: QIraf,TIraf
	1534	REAL(kind=kind(1.D0)),dimension(nlraf) :: TEraf
	1535	REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: RIraf
	1536	REAL(kind=kind(1.D0)) :: masseU,kbolt,Konst,g
	1537	REAL(kind=kind(1.D0)) :: facP,rhoNloc
	1538	masseU=1.e-3/RNAVO
	1539	kbolt=RKBOL
	1540	Konst=kbolt/masseU
	1541	g=RG
	1542
	1543
	1544	! below layer iz_vertplasma no effect due to vertical diffusion
	1545	do l=1,nlayer
	1546	if (P1D(l) .ge. Praf(1)) then
	1547	qk1dnew(l,:)=qk1d(l,:)
	1548	Ti1dnew(l,:)=Ti1d(l,:)
	1549	Te1dnew(l)=TE1d(l)
	1550	endif
	1551
	1552	if (P1D(l) .lt. Praf(1)) then
	1553	iP=1
	1554	do while(Praf(iP) .gt. P1D(l))
	1555	iP=iP+1
	1556	enddo
	1557	IP=iP-1
	1558	facP=(P1D(l)-Praf(iP))/(Praf(iP+1)-Praf(iP))
	1559	rhoNloc=P1D(l)/kbolt/TN1D(l)M1d(l)masseU
	1560
	1561	do k=1,ndiffions
	1562	Rknew(l,k)=RIraf(iP,k)+(RIraf(iP+1,k)-RIraf(iP,k))*facP
	1563	TI1dnew(l,k)=TIraf(IP,k)+(TIraf(iP+1,k)-TIraf(iP,k))*facP
	1564	enddo
	1565	TE1Dnew(l)=TEraf(iP)+(TEraf(iP+1)-TEraf(iP))*facP
	1566
	1567	do k=1,ndiffions
	1568	Qk1dnew(l,k)=Rknew(l,k)/rhoNloc
	1569	enddo
	1570
	1571	endif
	1572	enddo
	1573
	1574	END
	1575
	1576	SUBROUTINE GCMGRID_P2K(Zraf,Praf,Traf,Mraf,QIraf,RIraf,
	1577	& Tiraf,Teraf,Wk,Ez,P1Di,P1D,Tn1d,M1d,qk1d,qk1dnew,
	1578	& TI1D,TI1Dnew,TE1D,TE1Dnew,Wk1D,Wk1D2,Wk1D3,Ez1d,ZZ1d,
	1579	& nlraf,ndiffions,nlayer,izv,itrans,etrans,FacQ,FacTI,FacTe)
	1580
	1581	IMPLICIT NONE
	1582	#include "YOMCST.h"
	1583	INTEGER :: nlraf,ndiffions,nlayer,l,iP,k,izv
	1584	INTEGER :: itrans(ndiffions)
	1585	INTEGER :: etrans(1)
	1586	REAL(kind=kind(1.D0)),dimension(nlayer) :: P1d,Tn1d,M1d,ZZ1d
	1587	REAL(kind=kind(1.D0)),dimension(nlayer+1) :: P1di
	1588	REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: Qk1d,Qk1dnew
	1589	REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: TI1d,TI1dnew
	1590	REAL(kind=kind(1.D0)),dimension(nlayer) :: TE1D,TE1Dnew
	1591	REAL(kind=kind(1.D0)),dimension(nlayer+1,ndiffions) :: Wk1D,WK1d2
	1592	REAL(kind=kind(1.D0)),dimension(nlayer+1,ndiffions) :: Wk1D3
	1593	REAL(kind=kind(1.D0)),dimension(nlayer+1) :: Ez1D
	1594	! REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: Rknew
	1595	REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: FacQ,FacTI
	1596	REAL(kind=kind(1.D0)),dimension(nlayer) :: FacTE
	1597	REAL(kind=kind(1.D0)),dimension(nlraf) :: Zraf,Praf,Traf
	1598	REAL(kind=kind(1.D0)),dimension(nlraf) :: Mraf
	1599	REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: QIraf
	1600	REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: TIraf
	1601	REAL(kind=kind(1.D0)),dimension(nlraf) :: TEraf
	1602	REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: RIraf
	1603	REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: Wk
	1604	REAL(kind=kind(1.D0)),dimension(nlraf) :: Ez
	1605	REAL(kind=kind(1.D0)) masseU,kbolt,Konst,g
	1606	REAL(kind=kind(1.D0)) facP,rhoNloc
	1607	masseU=1.e-3/RNAVO
	1608	kbolt=RKBOL
	1609	Konst=kbolt/masseU
	1610	g=RG
	1611
	1612
	1613	! below layer iz_vertplasma no effect due to vertical diffusion
	1614	do l=1,nlayer
	1615	if (P1D(l) .ge. Praf(1)) then
	1616	qk1dnew(l,1:ndiffions)=qk1d(l,1:ndiffions)
	1617	Ti1dnew(l,1:ndiffions)=Ti1d(l,1:ndiffions)
	1618	endif
	1619
	1620	if (P1di(l) .ge. Praf(1)) then
	1621	Ez1D(l)=0D0
	1622	Wk1D(l,1:ndiffions)=0D0
	1623	Wk1D2(l,1:ndiffions)=0D0
	1624	Wk1D3(l,1:ndiffions)=0D0
	1625	endif
	1626
	1627	if (P1D(l) .lt. Praf(1)) then
	1628	iP=1
	1629	do while(Praf(iP) .gt. P1D(l))
	1630	iP=iP+1
	1631	enddo
	1632	iP=iP-1
	1633
	1634	! indP=maxloc(Praf,mask=Praf < P1D(l))
	1635	! IP=indP(1)-1
	1636
	1637	facP=(P1D(l)-Praf(iP))/(Praf(iP+1)-Praf(iP))
	1638	rhoNloc=P1D(l)/kbolt/TN1D(l)M1d(l)masseU
	1639	do k=1,ndiffions
	1640	Qk1dnew(l,k)=(RIraf(iP,k)+(RIraf(iP+1,k)-RIraf(iP,k))*facP)
	1641	& *facQ(l,k)/rhoNloc
	1642	TI1dnew(l,k)=(TIraf(IP,k)+(TIraf(iP+1,k)-TIraf(iP,k))*facP)
	1643	& *facTI(l,k)
	1644	enddo
	1645	TE1Dnew(l)=(TEraf(iP)+(TEraf(iP+1)-TEraf(iP))facP)facTE(l)
	1646
	1647	! if (l .eq. nlayer) then
	1648	! print*,'iP',iP,Praf(iP),P1D(l),Praf(iP+1),nlraf,facP
	1649	! print*,RIraf(iP,1),RIraf(iP+1,1),facQ(l,1),rhoNloc,facTI(l,1)
	1650	! print*,Qk1Dnew(l,1),TI1Dnew(l,1)
	1651	! endif
	1652
	1653	endif
	1654
	1655	if (P1Di(l) .lt. Praf(1)) then
	1656	iP=1
	1657	do while(Praf(iP) .gt. P1Di(l))
	1658	iP=iP+1
	1659	enddo
	1660	iP=iP-1
	1661
	1662	facP=(P1Di(l)-Praf(iP))/(Praf(iP+1)-Praf(iP))
	1663	do k=1,ndiffions
	1664	wk1d(l,k)=(Wk(iP,k)+(Wk(iP+1,k)-Wk(iP,k))*facP)
	1665	& /(ZZ1D(l)-ZZ1D(l-1))
	1666	wk1d2(l,k)=(Wk(iP,k)+(Wk(iP+1,k)-Wk(iP,k))*facP)
	1667	& *(ZZ1D(l)-ZZ1D(l-1))
	1668	wk1d3(l,k)=(Wk(iP,k)+(Wk(iP+1,k)-Wk(iP,k))*facP)
	1669	enddo
	1670	Ez1D(l)=(EZ(iP)+(Ez(iP+1)-Ez(iP))*facP)
	1671	endif
	1672
	1673	enddo ! l
	1674	Wk1d(nlayer+1,1:ndiffions)=0D0
	1675	Wk1d2(nlayer+1,1:ndiffions)=0D0
	1676	Wk1d3(nlayer+1,1:ndiffions)=0D0
	1677	Ez1D(nlayer+1)=0D0
	1678
	1679
	1680	END
	1681
	1682
	1683	SUBROUTINE CORRFACK(qbef,qaft,Tbef,Taft,Tebef,Teaf,
	1684	& Fq,Ft,Fte,nl,nq)
	1685
	1686	IMPLICIT NONE
	1687
	1688	INTEGER :: nl,nq,l,k
	1689	REAL(kind=kind(1.D0)),dimension(nl,nq) :: qbef,qaft,Fq
	1690	REAL(kind=kind(1.D0)),dimension(nl,nq) :: Tbef,Taft,Ft
	1691	REAL(kind=kind(1.D0)),dimension(nl) :: Tebef,Teaf,Fte
	1692
	1693	do l=1,nl
	1694	do k=1,nq
	1695	Fq(l,k)=qbef(l,k)/qaft(l,k)
	1696	Ft(l,k)=Tbef(l,k)/Taft(l,k)
	1697	enddo
	1698	Fte(l)=Tebef(l)/Teaf(l)
	1699	enddo
	1700
	1701	END
	1702
	1703	SUBROUTINE DCOEFFK(TI,NUI,M,DI,nl,nqdyn,nqtot,ind)
	1704
	1705	IMPLICIT NONE
	1706
	1707	INTEGER :: nl,l,nqdyn,nqtot,k,kn
	1708	INTEGER,dimension(nqdyn) :: ind
	1709	REAL(kind=kind(1.D0)),dimension(nl,nqdyn) :: NUI,TI,DI
	1710	REAL,dimension(nqtot) :: M
	1711	REAL(kind=kind(1.D0)) :: kbolt,masseU
	1712	masseU=1.660538782d-27
	1713	kbolt=1.3806504d-23
	1714
	1715	do k=1,nqdyn
	1716	kn=ind(k)
	1717	do l=1,nl
	1718	DI(l,k)=kbolt/masseU*TI(l,k)/NUI(l,k)/M(kn)
	1719	enddo
	1720	enddo
	1721
	1722	END
	1723
	1724	SUBROUTINE DIFFPARAMK(T,Te,RE,D,Dz,
	1725	& alpha,beta,gama,delta,eps,dzeta,eta,nl,Wmax)
	1726
	1727	IMPLICIT NONE
	1728
	1729	INTEGER :: nl,l
	1730	REAL(kind=kind(1.D0)),DIMENSION(nl) :: T,D
	1731	REAL(kind=kind(1.D0)),DIMENSION(nl) :: Te,Re
	1732	REAL(kind=kind(1.D0)) :: DZ,Wmax
	1733	REAL(kind=kind(1.D0)),DIMENSION(nl) :: alpha,beta,gama,delta,eps
	1734	REAL(kind=kind(1.D0)),DIMENSION(nl) :: dzeta,eta
	1735
	1736	alpha(1)=1D0/T(1)(-3D0T(1)+4D0T(2)-T(3))/(2D0DZ)+
	1737	& 1D0/T(1)(-3D0Te(1)+4D0Te(2)-Te(3))/(2D0DZ)
	1738	delta(1)=(-3D0D(1)+4D0D(2)-D(3))/(2D0*Dz)
	1739	beta(1)=1D0/T(1)
	1740	dzeta(1)=Te(1)/T(1)*
	1741	& (-3D0log(Re(1))+4D0log(Re(2))-log(Re(3)))/(2D0*DZ)
	1742
	1743	! Add a limit to dzeta (test to avoid unstabilities)
	1744
	1745	! if (abs(dzeta(1)) .ge. 1d0 .and. Wmax .ge. 500d0) then
	1746	! dzeta(1)=1d0*dzeta(1)/abs(dzeta(1))
	1747	! endif
	1748
	1749	! dzeta(1)=Te(1)/T(1)/Re(1)*
	1750	! & (-3D0Re(1)+4D0Re(2)-Re(3))/(2D0*Dz)
	1751
	1752	do l=2,nl-1
	1753	alpha(l)=1D0/T(l)(T(l+1)-T(l-1))/(2D0DZ)+
	1754	& 1D0/T(l)(Te(l+1)-Te(l-1))/(2D0DZ)
	1755	delta(l)=(D(l+1)-D(l-1))/(2D0*Dz)
	1756	beta(l)=1D0/T(l)
	1757	dzeta(l)=Te(l)/T(l)(log(Re(l+1))-log(Re(l-1)))/(2D0DZ)
	1758	! dzeta(l)=Te(l)/T(l)/Re(l)(Re(l+1)-Re(l-1))/(2D0DZ)
	1759
	1760	! if (abs(dzeta(l)) .ge. 1d0 .and. Wmax .ge. 500d0) then
	1761	! dzeta(l)=1d0*dzeta(l)/abs(dzeta(l))
	1762	! endif
	1763
	1764	enddo
	1765
	1766	! Upper altitude values
	1767
	1768	alpha(nl)=1D0/T(nl)(3D0T(nl)-4D0T(nl-1)+T(nl-2))/(2D0DZ)+
	1769	& 1D0/T(nl)(3D0Te(nl)-4D0Te(nl-1)+Te(nl-2))/(2D0DZ)
	1770	delta(nl)=(3D0D(nl)-4D0D(nl-1)+D(nl-2))/(2D0*DZ)
	1771	beta(nl)=1D0/T(nl)
	1772	dzeta(nl)=Te(nl)/T(nl)*
	1773	& (3D0log(Re(nl))-4D0log(Re(nl-1))+log(Re(nl-2)))/(2D0*DZ)
	1774
	1775	! if (abs(dzeta(nl)) .ge. 1d0 .and. Wmax .ge. 500d0) then
	1776	! dzeta(nl)=1d0*dzeta(nl)/abs(dzeta(nl))
	1777	! endif
	1778
	1779	! dzeta(nl)=Te(nl)/T(nl)/Re(nl)*
	1780	! & (3D0Re(nl)-4D0Re(nl-1)+Re(nl-2))/(2D0*Dz)
	1781
	1782	! Compute the gama and eps coefficients
	1783	! Lower altitude values
	1784
	1785	gama(1)=(-3D0beta(1) +4D0beta(2) -beta(3) )/(2d0*dz)
	1786	eps(1) =(-3D0alpha(1)+4D0alpha(2)-alpha(3))/(2d0*dz)
	1787	eta(1) =(-3D0dzeta(1)+4D0dzeta(2)-dzeta(3))/(2D0*dz)
	1788
	1789	do l=2,nl-1
	1790	gama(l)=(beta(l+1) -beta(l-1) )/(2d0*Dz)
	1791	eps(l) =(alpha(l+1)-alpha(l-1))/(2d0*Dz)
	1792	eta(l) =(dzeta(l+1)-dzeta(l-1))/(2D0*Dz)
	1793	end do
	1794
	1795	gama(nl)=(3D0beta(nl) -4D0beta(nl-1) +beta(nl-2) )/(2D0*DZ)
	1796	eps(nl) =(3D0alpha(nl)-4D0alpha(nl-1)+alpha(nl-2))/(2D0*DZ)
	1797	eta(nl) =(3D0dzeta(nl)-4D0dzeta(nl-1)+dzeta(nl-2))/(2D0*dz)
	1798	! print*,'test',alpha(nl-1),beta(nl-1),dzeta(nl-1)
	1799
	1800
	1801	END
	1802
	1803
	1804	SUBROUTINE MATCOEFFK(alpha,beta,gama,delta,eps,dzeta,eta,
	1805	& Dad,Rhoad,Facesc,dz,dt,A,B,C,D,Ti,Te,nl)
	1806
	1807	IMPLICIT NONE
	1808
	1809	INTEGER :: nl,l
	1810	REAL*8,DIMENSION(nl) :: alpha,beta,gama,delta,eps,dzeta,eta
	1811	REAL*8,DIMENSION(nl) :: Dad,RHoad,Ti,Te
	1812	REAL*8 :: dz,dt,del1,del2,del3,FacEsc
	1813	REAL*8,DIMENSION(nl) :: A,B,C,D
	1814	del1=dt/2d0/dz
	1815	del2=dt/dz/dz
	1816	del3=dt
	1817
	1818	! lower boundary coefficients no change
	1819	A(1)=0d0
	1820	B(1)=1d0
	1821	C(1)=0d0
	1822	D(1)=rhoAd(1)
	1823
	1824	do l=2,nl-1
	1825	A(l)=(delta(l)+(alpha(l)+beta(l)+dzeta(l))Dad(l))del1
	1826	& -Dad(l)*del2
	1827	B(l)=-(delta(l)*(alpha(l)+beta(l)+dzeta(l))
	1828	& +Dad(l)(gama(l)+eps(l)+eta(l)))del3
	1829	B(l)=B(l)+1d0+2d0Dad(l)del2
	1830	C(l)=-(delta(l)+(alpha(l)+beta(l)+dzeta(l))Dad(l))del1
	1831	& -Dad(l)*del2
	1832	D(l)=rhoAD(l)
	1833	enddo
	1834
	1835	! Upper boundary coefficients Diffusion profile
	1836	C(nl)=0d0
	1837	B(nl)=-1d0
	1838	A(nl)=exp(-dZ(alpha(nl-1)+beta(nl-1)+dzeta(nl-1)))FacEsc
	1839	! A(nl)=exp(-dZ(alpha(nl)+beta(nl)+dzeta(nl)))FacEsc
	1840	! A(nl)=exp(-dZ(alpha(nl-1)+beta(nl-1)))FacEsc
	1841	! A(nl)=exp(-dZ(alpha(nl-1)+beta(nl-1))
	1842	! & (Ti(nl-1)/(Te(nl-1)+Ti(nl-1))))*Facesc
	1843	D(nl)=0D0
	1844
	1845	END
	1846
	1847	SUBROUTINE CheckmassK(X,Y,nl)
	1848
	1849	IMPLICIT NONE
	1850
	1851	INTEGER :: nl
	1852	REAL*8,DIMENSION(nl) :: X,Y
	1853	REAL*8 Xtot,Ytot
	1854
	1855	Xtot=sum(X)
	1856	Ytot=sum(Y)
	1857
	1858	if (abs((Xtot-Ytot)/Xtot) .gt. 1d-4) then
	1859	print*,'no conservation for mass',Xtot,Ytot
	1860	endif
	1861
	1862	END
	1863
	1864

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

Download in other formats:

Original Format