! SUBROUTINE plasmaphys(q,rhoN,P,pk,teta,mmean,mudyn, ! & tempk,tetak,tempe,tetae,wcontk,wcovk,Efieldz,wphysk,alt) SUBROUTINE iondiff_red(ngrid,nlayer,nqmx,pplay,pplev, ! & ptneutre,ptelect,ption,pq, & ptneutre,pq,pphis, & gmtime,lat,lon, & ptimestep,pdteuv,pdtconduc,pdqdiff, & pdqiondiff) !**************************************************************** ! ! Ion ambipolar diffusion routine ! ! Authors: J-Y. Chaufray ! ------- ! ! J-Y. Chaufray ! ! Vitesse verticale et ses effets sur nk, Tk, uk, vk, etc... ! Seule les effets du terme (wk - wn) sont pris en compte ici ! Pour l instant seul les effets sur la densite nk sont pris en compte ! ! Version: 14/11/2020 ! ! ------- ! ! 2022/10/22: adapted and adding by Antoine Martinez ! !***************************************************************** USE chemparam_mod USE infotrac_phy, only: tname !USE dimphy use plasmaphys_venus_mod use iono_h, only: temp_elect, temp_ion IMPLICIT NONE ! Solve the vertical dynamics equation ! The equation is solved on a altitude refined grid after interpolation ! final results are reinterpolated on the atmospheric pressure grids ! Upper boundary ! The velocity at upper altitude should be parametrized (escape) ! Here I assume w2up (nlayer)=0 ! Lower boundaru ! W2(1) = O (The ion velocity is equal to neutral velocity) !#include "dimensions.h" !#include "paramet.h" !#include "plasmadyn.h" !#include "comgeom.h" ! ========================= ! ==== INPUT / OUTPUT ===== ! ========================= ! ====== INPUT ====== INTEGER,intent(in) :: ngrid ! number of atmospheric columns INTEGER,intent(in) :: nlayer ! number of atmospheric layers INTEGER,intent(in) :: nqmx ! number of advected tracers REAL,intent(in) :: ptimestep ! physical timestep (s) REAL,intent(in) :: pphis(ngrid) ! geopotential of the surface (m2/s) REAL,intent(in) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa) REAL,intent(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa) REAL,intent(in) :: pq(ngrid,nlayer,nqmx) ! mid-layer mass mixing ratio tracers (kg/kg) REAL,intent(in) :: ptneutre(ngrid,nlayer) ! mid-layer NEUTRAL temperature (K) ! REAL pdt(ngrid,nlayer) ! REAL,intent(in) :: ptelect(ngrid,nlayer) ! mid-layer ELECTRON temperature (K) ! REAL,intent(in) :: ption(ngrid,nlayer) ! mid-layer ION temperature (K) REAL,intent(in) :: pdteuv(ngrid,nlayer) ! Temperature EUV heating tendancy (K/s) REAL,intent(in) :: pdtconduc(ngrid,nlayer) ! Temperature conduction tendancy (K/s) REAL,intent(in) :: pdqdiff(ngrid,nlayer,nqmx) ! neutral mmr molecular diffusion tendancy (kg/kg/s) ! real pdq(ngrid,nlayer,nq) REAL,intent(in) :: gmtime ! day fraction ratio [from 0 to 1] REAL,intent(in) :: lat(ngrid) ! grid latitude [degree] REAL,intent(in) :: lon(ngrid) ! grid longitude [degree] ! ====== OUTPUT ====== REAL :: pdqiondiff(ngrid,nlayer,nqmx) ! ion mmr tendancy (kg/kg/s) ! ========================= ! ======== LOCAL ========== ! ========================= integer, parameter :: t_elec_origin= 2 integer, parameter :: t_ion_origin = 1 !Electronic temperature. Index 1 -> Theis et al. 1980 - model data ; Index 2-> Theis et al. 1984 - model data !Ion temperature. Index 1 -> from VIRA Model based on the model of Miller et al. 1980 & 1984. REAL :: lon_sun, sza_local,szad_local,mu_local REAL :: lon_local(ngrid), lat_local(ngrid) REAL :: ptelect(ngrid,nlayer) ! mid-layer ELECTRON temperature (K) REAL :: ption(ngrid,nlayer) ! mid-layer ION temperature (K) REAL :: qnew(ngrid,nlayer,nqmx), qold(ngrid,nlayer,nqmx) REAL :: mmean_new(ngrid,nlayer), vmr_new(ngrid,nlayer,nqmx) REAL :: vmr_ion(ngrid,nlayer) !REAL wcontk(ip1jmp1,llp+1,ndiffions),Efieldz(ip1jmp1,llp+1) !REAL wcovk(ip1jmp1,llp+1,ndiffions),wphysk(ip1jmp1,llp+1,ndiffions) ! 1D field I pass to double precision REAL(kind=kind(1.D0)) tdiff1,tdiff2 REAL(kind=kind(1.D0)),dimension(nlayer+1) :: P1D,Ez1d REAL(kind=kind(1.D0)),dimension(nlayer) :: Mmean1D,RhoN1D,qe1D REAL(kind=kind(1.D0)),dimension(nlayer) :: Pnlay1D,ZZ1D,TempN1d REAL(kind=kind(1.D0)),dimension(nlayer) :: TempE1D,TempE1dint REAL(kind=kind(1.D0)),dimension(nlayer) :: FacTe,TempE1dnew REAL(kind=kind(1.D0)),dimension(naltvert) :: Zraf,Praf,Tnraf,Mraf REAL(kind=kind(1.D0)),dimension(naltvert) :: Rraf,QEraf,REraf REAL(kind=kind(1.D0)),dimension(naltvert) :: RErafold,TEraf,Ez !REAL(kind=kind(1.D0)) X(lld),Y(lld),a(lld),b(lld),c(lld) !REAL(kind=kind(1.D0)) X2(lld),Y2(lld),a2(lld),b2(lld),c2(lld) INTEGER,dimension(1) :: indZ REAL(kind=kind(1.D0)) :: Rho0,Nu0,T0,H0,D0,dt,dtad,time0 REAL(kind=kind(1.D0)) :: tdiff3,tmin,Wmax,ttot,tsim,tdiff REAL(kind=kind(1.D0)) :: pphis1D ! geopotential of the surface (m2/s) REAL(kind=kind(1.D0)),save :: Wlim REAL(kind=kind(1.D0)),dimension(naltvert) :: RIad,Tnad,TIad,Tead REAL(kind=kind(1.D0)),dimension(naltvert) :: RElad,Dad,Zad REAL(kind=kind(1.D0)) :: Dzraf,Dz,FacEsc,Time REAL(kind=kind(1.D0)),dimension(naltvert) :: alpha,beta,delta REAL(kind=kind(1.D0)),dimension(naltvert) :: eps,gama,dzeta,eta REAL(kind=kind(1.D0)),dimension(naltvert) :: Atri,Btri,Ctri REAL(kind=kind(1.D0)),dimension(naltvert) :: Dtri,Xtri,Ytri INTEGER :: ig,k,l,ln,kn,ij0,l0,ln0,iter,ld,iq,istep,ke,niter,nn INTEGER :: iqelec INTEGER :: nsubreal LOGICAL,SAVE :: firstcall=.true. REAL(kind=kind(1.D0)) masse,invsgmu,Konst REAL(kind=kind(1.D0)) masseU,kBolt,g,RgazP,Rvenus,Grav,Mvenus,PI integer,save :: ntracers integer,dimension(:),allocatable,save :: gcmind ! array of GCM indexes real,dimension(:),allocatable,save :: mol_tr ! array of mol mass traceurs real,dimension(:),allocatable,save :: Charge_tr ! array of electric charge traceurs integer,dimension(nqmx) :: indic_tracers,indic_iondiff integer,parameter :: nb_espIon_diff=24 character(len=20),dimension(nb_espIon_diff) :: ListeIonDiff character(len=20) :: electrans(1) integer,dimension(:),allocatable,save :: itrans ! array of sub index gcmind of fluid ions integer,dimension(1),save :: etrans ! index of gcmind electron INTEGER,save :: iz_vertplasma ! vertical index above which ion diffusion is computed INTEGER,save :: ndiffions ! nombre d'ions decrit ! ---------- ALLOCATABLE ARRAY ------------------------------ ! (ngrid,nlayer,ntracers) !REAL(kind=kind(1.D0)),dimension(:,:),allocatable :: qnew !REAL(kind=kind(1.D0)),dimension(:,:),allocatable :: qold ! (nlayer,ntracers) REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: q1D REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: RhoK1D ! (nlayer,ndiffions) REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: qk1d REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: qk1dnew REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: qk1dint REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: TempI1D REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: TempI1din REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: TempI1dnw REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: FacQ REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: FacTI ! (nlayer+1,ndiffions) REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Wk1d REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Wk1d2 REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Wk1d3 ! (ngrid,ndiffions) REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: ueff ! (ndiffions) REAL(kind=kind(1.D0)),dimension(:),allocatable,save :: Mtot REAL(kind=kind(1.D0)),dimension(:),allocatable,save :: Mtot2 ! ---------- ALLOCATABLE ARRAY ------------------------------ ! (naltvert,ntracers) REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Qraf REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Rkraf ! (naltvert,ndiffions) REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: QIraf REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: RIraf REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: RIrafold REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: TIraf REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: FIraf REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: DIraf REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Rhock ! (ndiffions) REAL(kind=kind(1.D0)),dimension(:),allocatable,save :: MtotZ1 REAL(kind=kind(1.D0)),dimension(:),allocatable,save :: MtotZ2 ! ========================================================== ! ========================================================== ! ========== FIRST CALL AND MAIN PARAMETERS ================ ! ========================================================== ! ========================================================== ! ======= Initializations at first call ========== if (firstcall) then ! === On elimine la microphysique ========== ntracers=0 indic_tracers(1:nqmx)=0 ! type_tr ==> 1: neutral gas, 2: ion gas, 3: liquid, 10: others do iq=1,nqmx if ((type_tr(iq) .le. 2.) .and. (type_tr(iq) .gt. 0.)) then indic_tracers(iq)=1 ntracers=ntracers+1 endif enddo allocate(gcmind(ntracers)) allocate(mol_tr(ntracers)) allocate(Charge_tr(ntracers)) ! Store gcm indexes in gcmind and molar masses in mol_tr nn=0 do iq=1,nqmx if (indic_tracers(iq) .eq. 1) then nn=nn+1 gcmind(nn) = iq mol_tr(nn) = M_tr(iq) Charge_tr(nn) = 0.D0 endif enddo ! === On selectionne les especes ioniques ========== ndiffions = 0 indic_iondiff(1:nqmx) = 0 ! define used ions ListeIonDiff(1)="o2plus" ListeIonDiff(2)="oplus" ListeIonDiff(3)="cplus" ListeIonDiff(4)="nplus" ListeIonDiff(5)="co2plus" ListeIonDiff(6)="noplus" ListeIonDiff(7)="coplus" ListeIonDiff(8)="hplus" ListeIonDiff(9)="h2oplus" ListeIonDiff(10)="h3oplus" ListeIonDiff(11)="hcoplus" ListeIonDiff(12)="n2plus" ListeIonDiff(13)="ohplus" ! define electron electrans(1)="elec" ! seek the index and number of ion/electron do iq=1,ntracers do nn=1,nb_espIon_diff if (ListeIonDiff(nn) .eq. tname(gcmind(iq))) then indic_iondiff(iq)=iq ! 1 endif ! we don't diffuse electron but we readapt electron density at the end if (electrans(1) .eq. tname(gcmind(iq))) then etrans(1)=iq endif enddo if (indic_iondiff(iq) .ge. 1) then print*,'Ion specie ', tname(gcmind(iq)), & 'diffused in Iondiff_red' ndiffions=ndiffions+1 endif enddo print*,'number of diffused ion:',ndiffions allocate(itrans(ndiffions)) ! Store gcm ion indexes in itrans nn=0 do iq=1,ntracers if (indic_iondiff(iq) .ge. 1) then nn=nn+1 itrans(nn)=indic_iondiff(iq) ! iq Charge_tr(itrans(nn))=1.!qcharge endif enddo Charge_tr(etrans(1))=-1.!qcharge*(-1.) ! find vertical index above which ion diffusion is computed do l=1,nlayer if (pplay(1,l) .gt. Pdiffion) then iz_vertplasma = l endif enddo iz_vertplasma = iz_vertplasma+1 ! seuil de la diffusion verticale print*,'vertical index for ion diffusion', & iz_vertplasma,pplay(1,iz_vertplasma) ! iz_plasma = iz_vertplasma-3 ! seuil de la dynamique horizontale ! allocatation des tableaux dependants du nombre d especes diffusees allocate(ueff(ngrid,ndiffions)) ! ----------------------------------- ! !allocate-(qnew(ngrid,nlayer,ntracers)) !allocate-(qold(ngrid,nlayer,ntracers)) allocate(q1D(nlayer,ntracers)) allocate(RhoK1D(nlayer,ntracers)) allocate(qk1d(nlayer,ndiffions)) allocate(qk1dnew(nlayer,ndiffions)) allocate(qk1dint(nlayer,ndiffions)) allocate(TempI1D(nlayer,ndiffions)) allocate(TempI1din(nlayer,ndiffions)) allocate(TempI1dnw(nlayer,ndiffions)) allocate(FacQ(nlayer,ndiffions)) allocate(FacTI(nlayer,ndiffions)) allocate(Wk1d(nlayer+1,ndiffions)) allocate(Wk1d2(nlayer+1,ndiffions)) allocate(Wk1d3(nlayer+1,ndiffions)) allocate(Mtot(ndiffions)) allocate(Mtot2(ndiffions)) ! ----------------------------------- ! allocate(Qraf(naltvert,ntracers)) allocate(Rkraf(naltvert,ntracers)) allocate(QIraf(naltvert,ndiffions)) allocate(RIraf(naltvert,ndiffions)) allocate(RIrafold(naltvert,ndiffions)) allocate(TIraf(naltvert,ndiffions)) allocate(FIraf(naltvert,ndiffions)) allocate(DIraf(naltvert,ndiffions)) allocate(Rhock(naltvert,ndiffions)) allocate(MtotZ1(ndiffions)) allocate(MtotZ2(ndiffions)) ! ----------------------------------- ! ! Conditions aux limites Wlim=500. ! Maximal absolute value of vertical velocity (m/s) (to avoid unstabilities in wdu/dz term) ueff(1:ngrid,1:ndiffions)=0. ! Effusion velocity at top in m/s !ueff(1:ngrid,1:ndiffions)=200. firstcall=.FALSE. endif ! firstcall masseU = 1.660538782d-27 kBolt = 1.3806504d-23 Konst = kBolt/masseU RgazP = 8.314472 PI = 2.*ASIN(1.) ! 3.141592653 g = 8.87d0 Rvenus = 6051800d0 ! Used to compute escape parameter no need to be more accurate Grav = 6.67d-11 Mvenus = 4.86d24 ij0 = 6000 ! For test invsgmu= 1d0/g/masseU ke=etrans(1) alpha(1:naltvert)=0D0 beta(1:naltvert)=0D0 delta(1:naltvert)=0D0 eps(1:naltvert)=0D0 dzeta(1:naltvert)=0D0 eta(1:naltvert)=0D0 !! Open Vertical file to check it ! if (firstcall) then ! open(unit=301,file='Plasma_vertical.dat',status='unknown') ! firstcall=.FALSE. ! endif !! First we derive only parameters useful for diffusion !! Extrapolation of paramaters to all altitudes (we solve the dynamics equation on larger altitude range) !! !! ===> ON A PAS BESOIN DE CETTE PARTIE CAR ON EST DEJA DANS LA PHYSIQUE !! CALL DYN2DIFFK(Pk,teta,P,mmean,rhoN,Preff,cpp,kappa, !! & q,Mtraceur,Tempk,tetak,TempE,TetaE,Preffplasma,itrans,etrans, !! & llm,nqtot,llp,ndiffions,ip1jmp1,iz_plasma,kappak,kappae, !! & PNlay,TempN,Rhok,TempI,TetaI,rhotetaI,TempF,TetaF,rhotetaF) !! ! ========================================================== ! ========================================================== ! ========== Main Loop on horizontal grids ================= ! ========================================================== ! ========================================================== ! ========== Initialization ================= dt=dble(ptimestep) !dtplasma*iappvert nsubreal=nsubvert niter=2 qold(1:ngrid,1:nlayer,1:nqmx)=pq(1:ngrid,1:nlayer,1:nqmx) qnew(1:ngrid,1:nlayer,1:nqmx)=pq(1:ngrid,1:nlayer,1:nqmx) ! tdiff=dt/nsubreal ! print*,'dt',dt lon_sun = (0.5 - gmtime)*2.*PI lon_local(:) = lon(:)*PI/180. lat_local(:) = lat(:)*PI/180. ! ========== Horizontal loop start ========== DO ig=1,ngrid sza_local = cos(lat_local(ig))*cos(lon_local(ig)) & *cos(lon_sun) + cos(lat_local(ig)) & *sin(lon_local(ig))*sin(lon_sun) sza_local = min(sza_local,1.) sza_local = max(-1.,sza_local) sza_local = acos(sza_local) ! RAD szad_local = sza_local * 180. / PI ! Degree mu_local = cos(sza_local) ttot = dt tsim = 0. tdiff1 = dt/nsubvert ! tdiff2 = dt/nsubmin ! !if (nsubreal .ne. nsubvert) print*,'ig',ig,tdiff ZZ1D(1:nlayer) = 0D0 P1D(1:nlayer) = 0D0 Pnlay1D(1:nlayer) = 0D0 TempN1D(1:nlayer) = 0D0 Mmean1D(1:nlayer) = 0D0 RHON1D(1:nlayer) = 0D0 RHOK1D(1:nlayer,1:ntracers)=0D0 Q1D(1:nlayer,1:ntracers)=0D0 TEMPI1D(1:nlayer,1:ndiffions)=0D0 !TETAI1D(1:nlayer,1:ndiffions)=0D0 TEMPE1D(1:nlayer)=0D0 !TETAE1D(1:nlayer)=0D0 Zraf(1:naltvert)=0D0 Praf(1:naltvert)=0D0 Tnraf(1:naltvert)=0D0 Mraf(1:naltvert)=0D0 Qraf(1:naltvert,1:ntracers)=0D0 Rkraf(1:naltvert,1:ntracers)=0D0 REraf(1:naltvert)=0D0 QIraf(1:naltvert,1:ndiffions)=0D0 TIraf(1:naltvert,1:ndiffions)=0D0 FIraf(1:naltvert,1:ndiffions)=0D0 QEraf(1:naltvert)=0D0 TEraf(1:naltvert)=0D0 rhock(1:naltvert,1:ndiffions)=0D0 ! print*,'ig',ig,llm ! First we compute 1d fields in double precision ! The mixing ratios slightly change here (< 0.1%) pphis1D = pphis(ig) P1D(nlayer+1) = pplev(ig,nlayer+1) do l=1,nlayer P1D(l) = pplev(ig,l) Pnlay1D(l) = pplay(ig,l) TempN1D(l) = ptneutre(ig,l) do iq=1,ntracers q1d(l,iq)=qold(ig,l,gcmind(iq)) enddo ! compute the mean molecular mass Mmean1d(l) = (1D0/sum(q1d(l,:)/dble(mol_tr(:)))) ! Compute the total mass density (kg/m3) rhok1d(l,:) = 0D0 RhoN1D(l) = Pnlay1D(l)*Mmean1d(l)/TempN1D(l)/Konst do iq=1,ntracers rhok1d(l,iq)=q1d(l,iq)*RhoN1D(l) enddo do k=1,ndiffions kn=itrans(k) qk1d(l,k)=q1D(l,kn) !TempI1D(l,k)=ption(ig,l) !TetaI1D(l,k)=TetaI(ig,l,k) !rhotetaI1D(l,k)=rhotetaI(ig,l,k) enddo !ke=etrans(1) qe1D(l)=q1D(l,ke) !TempE1D(l)=ptelect(ig,l) !TetaE1D(l)=TetaF(ig,l) !rhoTetaE1D(l)=rhoTetaf(ig,l) enddo ! l (vertical levels) ! Compute colonne mass of each dynamic ion Mtot(1:ndiffions)=0D0 ! do k=1,ndiffions ! kn=itrans(k) ! do l=1,lld ! ln=l+iz_vertplasma ! print*,ln,llm ! Mtot(k)=Mtot(k)+qk1d(ln,k)*rhoN1d(ln) ! enddo ! enddo ! print*,'Mtot before',Mtot ! Compute altitude of the scalar grid pressure level !if (isnan(TempN1d(1))) then ! write (*,*) 'PROBLEME NAN: TITI' !endif CALL ZVERTK(Pnlay1D,TempN1D,mmean1d,ZZ1d,nlayer,pphis1D) ! Compute lectron and Ion temperature of the scalar grid pressure level do l=1,nlayer ptelect(ig,l) = temp_elect(ZZ1d(l)/1000.,TempN1d(l), & szad_local,t_elec_origin) ption(ig,l) = temp_ion(ZZ1d(l)/1000.,TempN1d(l), & szad_local,t_ion_origin) TempE1D(l) = ptelect(ig,l) !!! 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. do k=1,ndiffions !kn=itrans(k) TempI1D(l,k)=max(0.5*ptelect(ig,l),TempN1d(l)) enddo enddo ! l (vertical levels) !write(*,*) szad_local,ZZ1d(:)/1000.,ptelect(ig,:), ption(ig,:) ! do l=1,llm ! print*,'Pnlay1D, Plev',Pnlay1D(l),P1D(l),rhok1D(l,1:ntracers) ! enddo ! I use a better vertial resolution above iz_vertplasma and altitude grid interpolation indZ(1)=0 CALL UPPER_RESOLK(ZZ1D,Pnlay1D,TempN1d,mmean1d,rhok1d,TempI1d, & TempE1D,mol_tr,nlayer,ndiffions,ntracers,naltvert,iz_vertplasma, & itrans,etrans,Zraf,Tnraf,Praf,Mraf,Qraf,Rkraf, & QIraf,TIraf,FIraf,Qeraf,Teraf,indZ,gcmind) ! print*,'z last level, P last level',Zraf(:),Praf(:) DZraf=Zraf(2)-Zraf(1) do l=1,naltvert do k=1,ndiffions kn=itrans(k) RIraf(l,k)=Rkraf(l,kn) enddo !ke=etrans(1) REraf(l)=Rkraf(l,ke) enddo ! Test ! do l=1,naltvert ! print*,'l',l,Zraf(l),Praf(l),Tnraf(l),Qiraf(l,1),Tiraf(l,1), ! & Qeraf(l),Teraf(l),Riraf(l,1),Reraf(l) ! enddo ! Before solving the diffusion equation, we reddo interpolation on GCM grid ! to derive mass correction factor due to interpolation (for mass conservation purpose) CALL GCMGRID_PK(Zraf,Praf,TNraf,Mraf, & QIraf,RIraf,TIraf,Teraf, & P1D,Pnlay1D,TempN1d,Mmean1d, & qk1d,qk1dint,TempI1D,TempI1Din,TempE1d,TempE1dint, & naltvert,ndiffions,nlayer,iz_vertplasma,itrans,etrans) !DO ln=1,nlayer !IF (k.eq.1) THEN ! write(*,*) real(qk1Dint(ln,:)), qold(ig,ln,gcmind(ke)) !ENDIF !ENDDO FacQ(1:nlayer,1:ndiffions)=1D0 FacTI(1:nlayer,1:ndiffions)=1D0 FacTE(1:nlayer)=1D0 CALL CORRFACK(qk1d,qk1dint,TempI1d,TempI1din,TempE1d, & TempE1dint,FacQ,FacTI,FacTE,nlayer,ndiffions) CALL DCOEFFK(TIraf,FIraf,mol_tr,DIraf, & naltvert,ndiffions,ntracers,itrans) Rho0=1d-20 T0=Tnraf(naltvert)!TIraf(naltvert,1) do ln=1,naltvert Tnad(ln)=TNraf(ln)/T0 enddo ! if(abs(lon(ig)-(48.75)).le.0.1) then ! if(abs(lat(ig)-(-80.625)).le.0.1) then ! write(*,*) 'ZZ1d= ', ZZ1d(iz_vertplasma:90) ! write(*,*) 'Zraf= ', Zraf(1:naltvert) ! write(*,*) 'Pnlay1D= ', Pnlay1D(iz_vertplasma:90) ! write(*,*) 'Praf= ', Praf(1:naltvert) ! DO k=1,ndiffions ! kn=itrans(k) ! if((tname(gcmind(kn)).eq.'hplus').or. ! & (tname(gcmind(kn)).eq.'1oplus').or. ! & (tname(gcmind(kn)).eq.'1noplus')) then ! write(*,*) tname(gcmind(kn))//' DEBUT' ! write(*,*) 'avant= ', ! & qk1d(iz_vertplasma:90,k) ! write(*,*) 'apres= ', ! & qk1dint(iz_vertplasma:90,k) ! write(*,*) 'interm= ', ! & Rkraf(1:naltvert,kn) ! write(*,*) 'corrmass= ', ! & real(FacQ(iz_vertplasma:90,k)) ! endif ! ENDDO ! ! endif ! endif ! do l=1,llm ! print*,'effect interp',l,qk1d(l,2),qk1dint(l,2)*facQ(l,2), ! & TempI1d(l,2),TempI1din(l,2)*facTI(l,2),TempE1D(l), ! & TempE1Dint(l)*facTE(l),facQ(l,:),FacTI(l,:),FacTe(l) ! enddo ! stop ! Compute vertical velocity to derive subtime step Tdiff=1.E-9 DO k=1,ndiffions kn=itrans(k) H0=kbolt*T0/mol_tr(kn)/g/masseU D0=DIraf(naltvert,k) Time0=H0*H0/D0 Time=Tdiff/Time0 do ln=1,naltvert RIad(ln)=RIraf(ln,k)/Rho0 RElad(ln)=Reraf(ln)/rho0 TIad(ln)=TIraf(ln,k)/T0 TEad(ln)=TEraf(ln)/T0 Zad(ln)=Zraf(ln)/H0 Dad(ln)=DIraf(ln,k)/D0 enddo DZ=Zad(2)-Zad(1) CALL DIFFPARAMK(Tiad,Tead,RElad,Dad,DZ, & alpha,beta,gama,delta,eps,dzeta,eta,naltvert,Wmax) ! if (ig .eq. ij0) then ! print*,'zeta',dzeta(:) ! print*,'eta',eta(:) ! endif FacEsc=exp(-ueff(ig,k)*Time0/H0*DZ/Dad(naltvert-1)) CALL MATCOEFFK(alpha,beta,gama,delta,eps,dzeta,eta,Dad,Riad, & FacEsc,dZ,time,Atri,Btri,Ctri,Dtri,Tiad,Tead,naltvert) rhock(1,k)=0D0 rhock(2,k)=0D0 do ln=3,naltvert-1 if (ln .lt. naltvert-1) then rhock(ln,k)=-Dad(ln)*(-RIraf(ln+2,k)+8.*RIraf(ln+1,k) & -8.*RIraf(ln-1,k)+RIraf(ln-2,k))/RIraf(ln,k)/12D0/DZ & -Dad(ln)*(alpha(ln)+beta(ln)+dzeta(ln)) endif if (ln .eq. naltvert-1) then rhock(ln,k)=-Dad(ln)* & (log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz & -Dad(ln)*(alpha(ln)+beta(ln)+dzeta(ln)) endif rhock(ln,k)=D0/H0*rhock(ln,k) enddo rhock(naltvert,k)=ueff(ig,k) enddo Wmax=maxval(abs(rhock)) if (mu_local .ge. 0.3) tdiff3=h0*Dz/Wmax/300d0 !20d0 !/100d0 if (mu_local .le. -0.3) tdiff3=h0*dz/Wmax/300d0 !300d0 if (abs(mu_local) .lt. 0.3) tdiff3=h0*dz/Wmax/400d0 !/1000d0 tmin=max(tdiff1,tdiff3) tmin=min(tmin,tdiff2) nsubreal=anint(dt/tmin) tdiff=tmin !! This is a weird factor to optimize time calculation !! Here, it is to be sure that the first dtime is very low tmin = 2*0.45*(DZ*H0)**2. tmin = tmin/(maxval(abs(DIraf(:,:)))+(DZ*H0)*Wmax) if (tdiff .ge. tmin) then tdiff = tmin endif ! if (ig .eq. ij0) then ! print*,'nreal',nsubreal,nsubvert,nsubmin, ! & Wmax,tmin,tdiff,tdiff1,tdiff2,tdiff3,h0*dZ,tsim,ttot,h0,dZ ! write(301,*) nsubreal ! endif !tdiff=8.75D-5 !tmin tdiff=tmin ! Fix a minimal value for tdiff ! ========== TEMPORAL Loop ================= ! Loop over dynamic ions to solve the diffusion equation ! Compute total mass on the Z grid before diffusion ! MtotZ1(:)=0D0 ! do k=1,ndiffions ! do l=1,naltvert ! MtotZ1(k)=MtotZ1(k)+RIraf(l,k)*Dzraf ! enddo ! enddo ! print*,'MtotZ',MtotZ1 !if (szad_local.ge.93.) tsim = ttot ! DO istep=1,nsubreal DO WHILE (tsim .lt. ttot) ! if (ig .eq. ij0) then ! print*,'tsim',tsim ! do l=1,naltvert ! print*,'density1',l,REraf(l)/Mtraceur(31)/masseU/1d6, ! & Riraf(l,2)/Mtraceur(22)/masseU/1d6 ! enddo ! endif RErafold=REraf RIrafold=RIraf DO iter=1,niter ! needed to use approximatively electron density (t+dt) DO ln=1,naltvert RElad(ln)=REraf(ln)/Rho0 REraf(ln)=RErafold(ln) ! if (ig .eq. ij0) then ! print*,'density',RElad(ln)/Mtraceur(31)/masseU/1d6*rho0 ! & ,Riraf(ln,2)/Mtraceur(22)/masseU/1d6 ! endif ENDDO ! ln DO k=1,ndiffions kn=itrans(k) H0=kbolt*T0/mol_tr(kn)/g/masseU D0=DIraf(naltvert,k) Time0=H0*H0/D0 Time=Tdiff/Time0 ! Compute the diffusion coefficient using the collion frequency ! print*,'now adimension' ! Adimension all parameters by value at iz_vertplasma+1 do ln=1,naltvert RIad(ln)=RIraf(ln,k)/Rho0 ! RElad(ln)=REraf(ln)/Rho0 TIad(ln)=TIraf(ln,k)/T0 TEad(ln)=TEraf(ln)/T0 Zad(ln)=Zraf(ln)/H0 Dad(ln)=DIraf(ln,k)/D0 enddo ! ln DZ=Zad(2)-Zad(1) ! Tead(:)=0D0 CALL DIFFPARAMK(Tiad,Tead,RElad,Dad,DZ, & alpha,beta,gama,delta,eps,dzeta,eta,naltvert,Wmax) ! if (ig .eq. ij0) then ! print*,'zeta',dzeta(:) ! print*,'eta',eta(:) ! endif FacEsc=exp(-ueff(ig,k)*Time0/H0*DZ/Dad(naltvert-1)) ! & *Tiad(ln-1)/(Tiad(ln-1)+Tead(ln-1))) ! eta(:)=0D0 ! dzeta(:)=0D0 ! do l=1,naltvert ! print*,l,alpha(l),dzeta(l),eta(l),RElad(l),Riad(l) ! enddo ! stop CALL MATCOEFFK(alpha,beta,gama,delta,eps,dzeta,eta,Dad, & Riad,FacEsc,dZ,time,Atri,Btri,Ctri,Dtri,Tiad,Tead,naltvert) Xtri(:)=0D0 Xtri=Riad Ytri=Xtri ! Solve system !write(*,*) tsim, iter CALL TridagDP(Atri,Btri,Ctri,Dtri,Xtri,naltvert) !STOP ! Check mass conservation ! CALL CheckmassK(RIad,Xtri,naltvert) ! if (ig .eq. 1653) then ! do l=1,naltvert ! print*,'lphysnew',l,rho0*Xtri(l),k ! enddo ! endif do l=1,naltvert if (iter .eq. niter) RIraf(l,k)=Rho0*Xtri(l) if (RIraf(l,k) .lt. 0. .or. isnan(RIraf(l,k))) then ! print*,'phys q <0',ig,l,k,iter,istep,Xtri(l),Ytri(l),Riraf(l,k), ! & REraf(l),Tead(l)*T0,Tiad(l)*T0 RIraf(l,k)=1D-24 ! RIraf(l,k)=RIraf(l-1,k)*exp(-DZ*(alpha(l-1)+beta(l-1)+dzeta(l-1))) ! nsubreal=nsubreal+50 ! if (nsubreal .gt. 2000) then ! print*,'oula nsubreal',ig,nsubreal ! endif ! goto 1 endif enddo ! l ! Update electron mass density for next iteration DO l=1,naltvert REraf(l)=REraf(l)+0.5*Charge_tr(kn)*mol_tr(ke)/mol_tr(kn) & * rho0*(Xtri(l)-Ytri(l)) ENDDO ! l ! if (ig .eq. ij0 .and. k .eq. 2) then ! do l=1,naltvert ! print*,'l1',l,Reraf(l)/mol_tr(ke)/masseU/1d6, ! & Xtri(l)*rho0/mol_tr(kn)/masseU/1d6, ! & Ytri(l)*rho0/mol_tr(kn)/masseU/1d6 ! enddo ! l ! endif ! Compute vertical velocity in atmospheric frame IF (iter .eq. niter) THEN rhock(1,k)=0D0 rhock(2,k)=0D0 do ln=3,naltvert-1 ! rhock(ln,k)=-Dad(ln)*((RIraf(ln+1,k)-RIraf(ln-1,k))/2D0/Dz ! & +RIraf(ln,k)*(alpha(ln)+beta(ln)+dzeta(ln))) ! rhock(ln,k)=D0/H0*rhock(ln,k)/RIraf(ln,k) ! rhock(ln,k)=-Dad(ln)* ! & ((log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz ! & +(alpha(ln)+beta(ln))*(Tiad(ln)/(Tead(ln)+Tiad(ln)))) ! rhock(ln,k)=D0/H0*rhock(ln,k) if (ln .lt. naltvert-1) then ! rhock(ln,k)=-Dad(ln)*(-log(RIraf(ln+2,k))+8.*log(RIraf(ln+1,k)) ! & -8.*log(RIraf(ln-1,k))+log(RIraf(ln-2,k))) ! & /12D0/Dz rhock(ln,k)=-Dad(ln)*(-RIraf(ln+2,k)+8.*RIraf(ln+1,k) & -8.*RIraf(ln-1,k)+RIraf(ln-2,k))/RIraf(ln,k)/12D0/DZ & -Dad(ln)*(alpha(ln)+beta(ln)+dzeta(ln)) endif if (ln .eq. naltvert-1) then ! rhock(ln,k)=-Dad(ln)*((Tiad(ln)+Tead(ln))/Tiad(ln)* ! & (log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz) ! & -Dad(ln)*(alpha(ln)+beta(ln)) rhock(ln,k)=-Dad(ln)* & (log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz & -Dad(ln)*(alpha(ln)+beta(ln)+dzeta(ln)) endif rhock(ln,k)=D0/H0*rhock(ln,k) ! if (RIraf(ln,k) .eq. 1D-24 .or. RIraf(ln+1,k) .eq. 1D-24) ! & rhock(ln,k)=0. ! I limit the vertical velocity to 1km/s to avoid instabilities in velocity advection (TO BE IMPROVED) ! print*,'rhock',istep,ln,k,Zraf(ln),ln,rhock(ln,k),RIraf(ln,k), ! & Facesc ! & RIraf(ln,k)*exp(-DZ*(alpha(ln)+beta(ln)+dzeta(ln))), ! & (alpha(ln)+beta(ln)+dzeta(ln)), ! & (log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz,Atri(ln+1), ! & exp(-DZ*(alpha(ln)+beta(ln)+dzeta(ln))) ! if (abs(rhock(ln,k)) .gt. 1000. .or. isnan(rhock(ln,k))) then ! print*,ig,istep ! print*,'rhock',Zraf(ln),ln,rhock(ln,k),RIraf(ln+1,k) ! & RIraf(ln,k)*exp(-DZ*(alpha(ln)+beta(ln)+dzeta(ln))), ! & (alpha(ln)+beta(ln)+dzeta(ln)), ! & (log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz,Atri(ln+1), ! & exp(-DZ*(alpha(ln)+beta(ln)+dzeta(ln))) ! if (abs(rhock(ln,k)) .gt. 10000.) stop ! rhock(ln,k) =1000.*rhock(ln,k)/abs(rhock(ln,k)) ! endif enddo ! ln rhock(naltvert,k)=ueff(ig,k) ENDIF ! (iter .eq. niter) ENDDO ! dynamic ions ! Upadte the electron density to assume neutral equilibrium IF (iter .eq. niter) THEN RERaf=Rerafold !ke=etrans(1) DO l=1,naltvert DO k=1,ndiffions kn=itrans(k) REraf(l)=REraf(l)+Charge_tr(kn)*mol_tr(ke)/mol_tr(kn) & * (RIraf(l,k)-Rirafold(l,k)) ENDDO ! if (ig .eq. ij0) then ! DO k=1,ndiffions ! kn=itrans(k) ! print*,'fin cycle',l,Rerafold(l)/Mtraceur(ke)/masseU/1d6, ! & Reraf(l)/Mtraceur(ke)/masseU/1d6, ! & Rirafold(l,k)/Mtraceur(kn)/masseU/1d6, ! & Riraf(l,k)/Mtraceur(kn)/masseU/1d6 ! enddo ! endif ENDDO ENDIF ! (iter .eq. niter) ! Compute vertical component of electric field ! DZ=Zraf(2)-Zraf(1) Ez(1)=0D0 DO ln=2,naltvert-1 Ez(ln)=(log(Teraf(ln+1))-log(Teraf(ln-1)))+ & (log(Reraf(ln+1))-log(Reraf(ln-1))) Ez(ln)=-kbolt/qcharge*Teraf(ln)/2D0/Dz/H0*Ez(ln) ENDDO Ez(naltvert)=3D0*log(Reraf(naltvert))-4D0*log(Reraf(naltvert-1)) & +log(Reraf(naltvert-2))+ & 3D0*log(Teraf(naltvert))-4D0*log(Teraf(naltvert-1)) & +log(Teraf(naltvert-2)) Ez(naltvert)=-kbolt/qcharge*Teraf(naltvert)/2D0/Dz/H0*Ez(naltvert) ! Ez(naltvert)=0D0 ! if (ig .eq. ij0) then ! do ln=1,naltvert ! ! print*,'reraf write',ln,Reraf(ln)/Mtraceur(31)/masseU/1d6 ! write(301,*) ig,istep,ln,Zraf(ln),RIraf(ln,:),rhock(ln,:), ! & REraf(ln),TIraf(ln,:),Teraf(ln),Tnraf(ln) ! enddo ! endif ENDDO ! iteration tsim=tsim+tdiff ! compute new time step Wmax=maxval(abs(rhock)) ! if (mu_local .ge. 0.) tdiff3=Dz/Wmax/300d0 ! if (mu_local .lt. 0.) tdiff3=dz/Wmax/500d0 !if (mu_local .ge. 0.5) tdiff3=h0*Dz/Wmax/20d0!20d0 !20d0 !/100d0 !if (mu_local .le. -0.5) tdiff3=2d0*16d0*h0*dz/Wmax/20d0!20d0 !300d0 !if (abs(mu_local) .lt. 0.5) tdiff3=2d0*16d0*h0*dz/Wmax/20d0!20d0 !400d0 !/1000d0 !! This is a weird factor to optimize time calculation if (mu_local .ge. 0.3) tdiff3=2d0*16d0*h0*Dz/Wmax/20d0!20d0 !20d0 !/100d0 if (mu_local .le. -0.3) tdiff3=2d0*16d0*h0*dz/Wmax/20d0!20d0 !300d0 if (abs(mu_local) .lt. 0.3) tdiff3=2d0*16d0*h0*dz/Wmax/20d0!20d0 !400d0 !/1000d0 tmin=max(tdiff1,tdiff3) tmin=min(tmin,tdiff2) tdiff=tmin !! This is a weird factor to optimize time calculation !!tmin = 2*0.45*(DZ*H0)**2./maxval(abs(DIraf(:,:))) tmin = 3.*8.*4.*200*2*0.45*(DZ*H0)**2. tmin = tmin/(maxval(abs(DIraf(:,:)))+(DZ*H0)*Wmax) if (tdiff .ge. tmin) then tdiff = tmin endif tdiff=max(tdiff1,tdiff) !if (szad_local .le. 20.) tdiff = tdiff/20. !!if (tsim .le. 8.75D-2) tdiff = 8.75D-5 if (tsim+tdiff .gt. ttot) tdiff=ttot-tsim ! if (ig .eq. ij0) then ! print*,'dt',Wmax,tdiff1,tdiff2,tdiff3,tdiff,tsim,ttot, ! & alpha(naltvert-1),beta(naltvert-1),dzeta(naltvert-1), ! & Atri(naltvert),Xtri(naltvert),Xtri(naltvert-1)*Atri(naltvert), ! & h0,dz,h0*dz,mu_local ! endif ENDDO ! time step ou while ! nsubreal=nsubvert ! end time step here ! if (ig .eq. ij0) close(301) MtotZ2(:)=0D0 ! Limites pour la vitesse verticale pour eviter problemes numeriques DO k=1,ndiffions DO ln=1,naltvert IF (abs(rhock(ln,k)).gt.Wlim .or. isnan(rhock(ln,k))) THEN rhock(ln,k) = Wlim*rhock(ln,k)/abs(rhock(ln,k)) ENDIF !IF (k.eq.1) THEN ! write(*,*) real(qk1Dnew(ln,:)), qold(ig,ln,gcmind(ke)) !ENDIF ENDDO ! ln ENDDO ! k ! print*,'rhock',l,k,Zraf(l),l,rhock(l,k),RIraf(l,k),Facesc ! MtotZ2(k)=MtotZ2(k)+RIraf(l,k)*Dzraf ! enddo ! enddo ! We reinterpolate results on the Pressure grid and correct with FacQ CALL GCMGRID_P2K(Zraf,Praf,TNraf,Mraf,QIraf,RIraf,TIraf,TEraf, & rhock,Ez,P1D,Pnlay1D,TempN1d,Mmean1d,qk1d,qk1dnew, & TempI1D,TempI1Dnw,TempE1D,TempE1dnew,Wk1D,Wk1D2,Wk1d3, & Ez1d,ZZ1d,naltvert,ndiffions,nlayer, & iz_vertplasma,itrans,etrans,FacQ,FacTI,FacTe) !iq = itrans(1) !write(*,*) tname(gcmind(iq)),ptimestep,qk1dnew(:,1) ! do l=1,llm+1 ! print*,'Ez,W',P1D(l),Ez1D(l),Wk1D(l,:) ! enddo ! Compute the mixing ratio and temperature & Update potential temperature ! Here I neglect advection of temperature (later) ! do ln=1,llm ! l=ln-iz_plasma ! Alt(ig,ln)=ZZ1d(ln) ! do k=1,ndiffions ! kn=itrans(k) !! if (ig .eq. ij0) print*,'qold',ln,q(ig,ln,kn),qk1d(ln,k) ! q(ig,ln,kn)=real(qk1Dnew(ln,k)) !! if (ig .eq. ij0) print*,'qnew',q(ig,ln,kn),qk1dnew(ln,k) ! if (l .gt. 0) then ! Tetak(ig,l,k)=real(tempI1Dnw(ln,k))**(1D0-kappak(k))* ! & (masseU*Mtraceur(kn)/kbolt*Preffplasma/RhoN(ig,ln)/q(ig,ln,kn)) ! & **(kappak(k)) ! wcontk(ig,l,k)=real(Wk1D(ln,k)) ! wcovk(ig,l,k)=real(Wk1D2(ln,k)) ! wphysk(ig,l,k)=real(Wk1D3(ln,k)) ! EfieldZ(ig,l)=real(Ez1D(ln)) ! endif ! ! if (q(ig,ln,kn) .lt. 0. .or. isnan(q(ig,ln,kn))) then ! print*,'aie q < 0',q(ig,ln,kn),ln,l,ig,k,kn ! print*,qk1dnew(:,k),qk1d(:,k) ! stop ! endif ! ! if (l .ge. 1) then ! if (isnan(Wcontk(ig,l,k))) then ! Wcontk(ig,l,k)=0. ! Wcovk(ig,l,k)=0. ! ! print*,'aie Wcontk',l,ln,ig,Wcontk(ig,l,k),Wk1D(ln,k) ! ! do l=1,naltvert ! ! print*,'rhock',l,rhock(l,:),RIraf(l,:),Dad(l),alpha(l),beta(l), ! ! & dzeta(l),Praf(l) ! ! enddo ! endif ! if (isnan(Efieldz(ig,l))) then ! Efieldz(ig,l)=0. ! endif ! endif ! ! enddo ! enddo ! do k=1,ndiffions ! Wcontk(ig,llp+1,k)=real(Wk1D(llm+1,k)) ! Wcovk(ig,llp+1,k)=real(Wk1D2(llm+1,k)) ! wphysk(ig,llp+1,k)=real(Wk1D3(llm+1,k)) ! enddo ! EfieldZ(ig,llp+1)=Ez1D(llm+1) DO ln=1,nlayer l=ln-iz_vertplasma DO k=1,ndiffions iq=gcmind(itrans(k)) qnew(ig,ln,iq)=real(qk1Dnew(ln,k)) ENDDO ENDDO ENDDO ! END ig - Main Loop on horizontal grids - plan horizontal ! ========================================================== ! ===== Compute the diffusion trends du to diffusion ======= ! ================== for the ions ========================== ! ========================================================== !ke = etrans(1) DO l=1,nlayer DO k=1,ndiffions iq = gcmind(itrans(k)) pdqiondiff(:,l,iq) = (qnew(:,l,iq)-qold(:,l,iq))/ptimestep ! CE CALCUL EST FAUX, CE N'EST PAS UNE CONSERVATION MMR QU'IL FAUT MAIS UNE CONSERVATION VMR !!!!!!!!!! ! pdqiondiff(:,l,iqelec)= pdqiondiff(:,l,iqelec) ! & + pdqiondiff(:,l,iq) !pdqiondiff(:,l,iq) =(qnew(:,l,kn)-qold(:,l,kn))/ptimestep !pdqiondiff(:,l,iqelec)= pdqiondiff(:,l,ke)+pdqiondiff(:,l,iq) ENDDO ENDDO !iq = gcmind(itrans(1)) !write(*,*) tname(iq),iq,ptimestep,pdqiondiff(1,:,iq),qold(1,:,iq) !write(*,*) maxval(pdqiondiff(:,:,iq)) ! ========================================================== ! ===== Compute the diffusion trends du to diffusion ======= ! ================== for the electron ====================== ! ========================================================== ! ------ update mmean ------ ! mmean_new(:,:) = 0. do iq = 1,ntracers mmean_new(:,:) = mmean_new(:,:) + & qnew(:,:,gcmind(iq))/M_tr(gcmind(iq)) enddo mmean_new(:,:) = 1./mmean_new(:,:) ! ------ conversion mmr into vmr ------ ! do iq = 1,ntracers vmr_new(:,:,gcmind(iq)) = qnew(:,:,gcmind(iq)) * & mmean_new(:,:)/M_tr(gcmind(iq)) enddo iqelec = gcmind(ke) ! ------ vmr ion density ------ ! vmr_ion(:,:) = 0. do iq = 1,ntracers if ((type_tr(gcmind(iq)) .eq. 2) .and. & ( tname(gcmind(iq)) .ne. tname(iqelec))) then vmr_ion(:,:) = vmr_ion(:,:) + vmr_new(:,:,gcmind(iq)) endif enddo ! ------ force charge neutrality ------ ! ! ------ vmr(ion) = vmr(elec) ------ ! DO ig=1,ngrid DO l=1,nlayer if (vmr_new(ig,l,iqelec) .ne. vmr_ion(ig,l)) then vmr_new(ig,l,iqelec) = vmr_ion(ig,l) ! DO k=1,ndiffions ! vmr_new(ig,l,iqelec) = vmr_new(ig,l,iqelec) + ! & vmr_new(ig,l,gcmind(itrans(k))) ! ENDDO endif ENDDO ENDDO ! ------ conversion vmr into mmr for electron ------ ! qnew(:,:,iqelec) = vmr_new(:,:,iqelec)*m_tr(iqelec)/mmean_new(:,:) ! ------ Compute the diffusion trends du to diffusion ------ !* ! ------ for the electron ------ ! DO l=1,nlayer pdqiondiff(:,l,iqelec) = & (qnew(:,l,iqelec)-qold(:,l,iqelec))/ptimestep ENDDO ! write(*,*) 'TATA EST PLUS LA' ! do ig=1,ip1jm+1,iip1 ! do l=1,llp ! ln=l+iz_plasma ! do k=1,ndiffions ! kn=itrans(k) ! q(ig+iim,ln,kn)=q(ig,ln,kn) ! tetak(ig+iim,l,k)=tetak(ig,l,k) ! wcontk(ig+iim,l,k)=wcontk(ig,l,k) ! if (ig .eq. ij0 .and. l .eq. 2 .and. k .eq.1) ! & print*,ij0+iim,wcontk(ij0,2,1),wcontk(ij0+iim,2,1) ! wcovk(ig+iim,l,k)=wcovk(ig,l,k) ! wphysk(ig+iim,l,k)=wphysk(ig,l,k) ! Efieldz(ig+iim,l)=Efieldz(ig,l) ! tetae(ig+iim,l)=tetae(ig,l) ! enddo ! enddo ! enddo ! print*,'wphys3',ij0,ij0+iim,wcontk(ij0,2,1),wcontk(ij0+iim,2,1) RETURN END ! ******************************************************************** ! ******************************************************************** ! ******************************************************************** SUBROUTINE TridagDP(a,b,c,r,u,n) ! USE mod_phys_lmdz_para, only: mpi_rank ! parameter (nmax=5000) ! dimension gam(nmax),a(n),b(n),c(n),r(n),u(n) real(kind=kind(1.D0)) gam(n),a(n),b(n),c(n),r(n),u(n) if(b(1).eq.0.)then stop 'tridag: error: b(1)=0 !!! ' endif bet=b(1) u(1)=r(1)/bet do 11 j=2,n gam(j)=c(j-1)/bet bet=b(j)-a(j)*gam(j) if(bet.eq.0.) then stop 'tridag: error: bet=0 !!! ' endif u(j)=(r(j)-a(j)*u(j-1))/bet 11 continue do 12 j=n-1,1,-1 u(j)=u(j)-gam(j+1)*u(j+1) 12 continue return END ! SUBROUTINE DYN2DIFFK(Pk,teta,P,mmean,rhoN,Preff,cpp,kappa, ! & q,Mtraceur,Tempk,tetak,TempE,tetaE,Preffplasma,itrans,etrans, ! & llm,nqtot,llp,ndiffions,ip1jmp1,iz_plasma,kappak,kappae, ! & Pnlay,TempN,Rhok,TempI,TetaI,rhotetaI,TempF,TetaF,rhotetaF) ! ! USE infotrac, only: tname ! ! IMPLICIT NONE ! ! INTEGER :: llm,llp,nqtot,ndiffions,ip1jmp1,l,ig,k,ln,kn,iz_plasma ! REAL :: Pk(ip1jmp1,llm),teta(ip1jmp1,llm) ! REAL :: mmean(ip1jmp1,llm),rhoN(ip1jmp1,llm) ! REAL :: P(ip1jmp1,llm+1),Pnlay(ip1jmp1,llm) ! REAL :: q(ip1jmp1,llm,nqtot) ! REAL :: Tempk(ip1jmp1,llp,ndiffions),tetak(ip1jmp1,llp,ndiffions) ! REAL :: TempE(ip1jmp1,llp),TetaE(ip1jmp1,llp) ! REAL :: Mtraceur(nqtot) ! INTEGER :: itrans(ndiffions),etrans(1) ! REAL :: Preff,Cpp,Preffplasma,kappa,unskappa ! REAL :: kappak(ndiffions) ! REAL :: kappae ! ! ! Output ! REAL :: TempN(ip1jmp1,llm) ! REAL :: Rhok(ip1jmp1,llm,nqtot) ! REAL :: TempI(ip1jmp1,llm,ndiffions) ! REAL :: TetaI(ip1jmp1,llm,ndiffions) ! REAL :: rhotetaI(ip1jmp1,llm,ndiffions) ! REAL :: TempF(ip1jmp1,llm) ! REAL :: TetaF(ip1jmp1,llm) ! REAL :: rhotetaF(ip1jmp1,llm) ! ! REAL :: masseU,kbolt ! masseU=1.660538782d-27 ! kbolt=1.3806504d-23 ! ! unskappa=1./kappa ! ! Rhok(:,:,:)=0. ! DO ig=1,ip1jmp1 ! DO ln=1,llm ! l=ln-iz_plasma ! TempN(ig,ln)=Teta(ig,ln)*pk(ig,ln)/Cpp ! Pnlay(ig,ln)=preff*(pk(ig,ln)/cpp)**unskappa ! do k=1,nqtot ! if (tname(k) .ne. "dust_number" .and. tname(k) .ne. "dust_mass" ! & .and. tname(k).ne."ccn_number".and.tname(k).ne."ccn_mass") then ! if (q(ig,ln,k) .le. 0.) q(ig,ln,k)=1d-30 ! Rhok(ig,ln,k)=q(ig,ln,k)*RhoN(ig,ln) !! if (q(ig,ln,k) .le. 0. .and. l .le. 0) Rhok(ig,ln,k)=1d-30*rhoN(ig,ln) ! endif ! enddo ! !! Ions parameters ! do k=1,ndiffions ! kn=itrans(k) ! if (ln .le. iz_plasma) TempI(ig,ln,k)=TempN(ig,ln) ! if (ln .gt. iz_plasma) TempI(ig,ln,k)=Tempk(ig,l,k) ! ! if (ln .le. iz_plasma) TetaI(ig,ln,k)=teta(ig,ln)* ! & (Preffplasma**kappak(k))/(Preff**kappa)* !! & (Mtraceur(kn)/mmean(ig,ln)/q(ig,ln,kn))**kappa ! & (rhoN(ig,ln)/mmean(ig,ln))**kappa* ! & 1D0/(rhok(ig,ln,kn)/Mtraceur(kn))**kappak(k)* ! & (kbolt*TempN(ig,ln)/masseu)**(kappa-kappak(k)) ! if (ln .gt. iz_plasma) TetaI(ig,ln,k)=Tetak(ig,l,k) ! ! rhotetaI(ig,ln,k)=rhok(ig,ln,kn)*TetaI(ig,ln,k) ! enddo ! !! Electron parameters ! ! kn=etrans(1) ! if (ln .le. iz_plasma) TempF(ig,ln)=TempN(ig,ln) ! if (ln .gt. iz_plasma) TempF(ig,ln)=TempE(ig,l) ! ! if (ln .le. iz_plasma) TetaF(ig,ln)=teta(ig,ln)* ! & (Preffplasma**kappae)/(Preff**kappa)* !! & (Mtraceur(kn)/mmean(ig,l)/q(ig,ln,kn))**kappa ! & (rhoN(ig,ln)/mmean(ig,ln))**kappa* ! & 1D0/(rhok(ig,ln,kn)/Mtraceur(kn))**kappae* ! & (kbolt*TempN(ig,ln)/masseU)**(kappa-kappae) ! if (ln .gt. iz_plasma) TetaF(ig,ln)=TetaE(ig,l) ! ! rhotetaF(ig,ln)=rhok(ig,ln,kn)*TetaF(ig,ln) ! ! ! ENDDO ! l or l=ln-iz_plasma ! ENDDO ! ig ! END SUBROUTINE ZVERTK(P,T,M,Z,nl,gsurf) IMPLICIT NONE #include "YOMCST.h" integer :: nl,l REAL(kind=kind(1.D0)):: P(nl),T(nl),M(nl),Z(nl) REAL(kind=kind(1.D0)):: masseU,kbolt,Konst,g,z0,Hpm REAL(kind=kind(1.D0)):: Tmean,Mmean,gsurf masseU=1.e-3/RNAVO kbolt=RKBOL Konst=kbolt/masseU g=RG Z0=gsurf/g!0d0 Z(1)=z0 Hpm=Konst*T(1)/g/M(1) do l=2,nl Tmean=T(l-1) Mmean=M(l-1) Hpm=Konst*Tmean/g/Mmean Z(l)=Z(l-1)-Hpm*log(P(l)/P(l-1)) ! print*,'z',l,Z(l) enddo END SUBROUTINE UPPER_RESOLK(ZZ1D,P1D,TN1d,M1d,rhok1d,TI1d,Te1d, & mol_tr,nlayer,ndiffions,ntracers,nlraf,idiffZ, & itrans,etrans,Zraf,TNraf,Praf,Mraf,Qraf,rhokraf, & QIraf,TIraf,FIraf,QEraf,Teraf,IndZ,gcmind) USE infotrac_phy, only: tname IMPLICIT NONE #include "YOMCST.h" INTEGER :: nlayer,ndiffions,ntracers,nlraf,idiffZ,iq,k,kn,ke,iz,l INTEGER :: indZ(1) INTEGER :: itrans(ndiffions),etrans(1),gcmind(ntracers) REAL :: mol_tr(ntracers) REAL(kind=kind(1.D0)),dimension(nlayer) :: P1D,TN1D,M1D,ZZ1D,Te1D REAL(kind=kind(1.D0)),dimension(nlayer,ntracers) :: rhok1d REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: TI1d REAL(kind=kind(1.D0)),dimension(nlraf) :: Zraf,Praf,TNraf REAL(kind=kind(1.D0)),dimension(nlraf) :: Mraf REAL(kind=kind(1.D0)),dimension(nlraf,ntracers) :: Qraf,Rhokraf REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: QIraf,FIraf REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: TIraf REAL(kind=kind(1.D0)),dimension(nlraf) :: Qeraf,TEraf REAL(kind=kind(1.D0)) :: kbolt,masseU,Konst,g,a0,apol REAL(kind=kind(1.D0)) :: mu,nlay,facZ,dZ masseU=1.e-3/RNAVO kbolt=RKBOL Konst=kbolt/masseU g=RG a0=2.9D-15 ! absolute constant (when densities expressed in /m3) !apol=2.911D0 ! in 10^-24 cm3 (Withers 2009 for CO2) Zraf(1) = ZZ1D(idiffZ) Praf(1) = P1D(idiffZ) TNraf(1) = TN1D(idiffZ) Qraf(:,:) = 0d0 Rhokraf(1:nlraf,1:ntracers)=0d0 do iq=1,ntracers ! if (tname(iq) .ne. "dust_number" .and. tname(iq) .ne. "dust_mass" ! & .and. tname(iq).ne."ccn_number".and.tname(iq).ne."ccn_mass") then Rhokraf(1,iq)=rhok1d(idiffZ,iq) Qraf(1,iq)=rhok1d(idiffZ,iq)/sum(rhok1d(idiffZ,1:ntracers)) ! endif enddo Mraf(1)=1D0/sum(Qraf(1,1:ntracers)/mol_tr(1:ntracers)) nlay=Praf(1)/kbolt/TNraf(1) do k=1,ndiffions kn=itrans(k) TIraf(1,k)=TI1D(idiffZ,k) QIraf(1,k)=Qraf(1,kn) ! mu=mol_tr(kn)*Mraf(1)/(mol_tr(kn)+Mraf(1)) ! FIraf(1,k)=nlay*a0*Mraf(1)/(Mraf(1)+mol_tr(kn))*sqrt(apol/mu) enddo ke=etrans(1) Teraf(1)=TE1D(idiffZ) Qeraf(1)=Qraf(1,ke) Zraf(nlraf)=ZZ1d(nlayer) do l=2,nlraf-1 Zraf(l)=Zraf(1)+(Zraf(nlraf)-Zraf(1))/DBLE(nlraf-1)*DBLE(l-1) iZ=1 do while (ZZ1D(iz) .le. Zraf(l)) iZ=iZ+1 enddo Iz=iz-1 ! indZ=maxloc(ZZ1D,MASK=ZZ1d < Zraf(l)) ! print*,'indZ',indZ ! iZ=indZ(1) dZ=Zraf(l)-Zraf(l-1) facZ=(Zraf(l)-zz1d(iz))/(zz1d(iZ+1)-ZZ1d(iz)) TNraf(l)=TN1D(iz)*(TN1D(iz+1)/TN1d(iz))**facZ do iq=1,ntracers if (Rhok1d(iz,iq) .gt. 0.) then Rhokraf(l,iq)=Rhok1d(iz,iq) * & (Rhok1d(iz+1,iq)/Rhok1d(iz,iq))**facZ endif enddo do iq=1,ntracers Qraf(l,iq)=Rhokraf(l,iq)/sum(Rhokraf(l,1:ntracers)) enddo Mraf(l)=1D0/sum(Qraf(l,1:ntracers)/mol_tr(1:ntracers)) nlay=sum(Rhokraf(l,1:ntracers))/Mraf(l)/masseU Praf(l)=nlay*kbolt*TNraf(l) do k=1,ndiffions kn=itrans(k) TIraf(l,k)=TI1D(iz,k)*(TI1D(iz+1,k)/TI1D(iz,k))**facZ Qiraf(l,k)=Qraf(l,kn) ! mu=mol_tr(kn)*Mraf(l)/(mol_tr(kn)+Mraf(l)) ! FIraf(l,k)=nlay*a0*Mraf(l)/(Mraf(l)+mol_tr(kn))*sqrt(apol/mu) enddo Teraf(l)=TE1D(iz)*(TE1D(iz+1)/TE1D(iz))**facZ Qeraf(l)=Qraf(l,ke) enddo Zraf(nlraf)=ZZ1d(nlayer) TNraf(nlraf)=TN1D(nlayer) ! Qraf(:,:)=0d0 do iq=1,ntracers Rhokraf(nlraf,iq)=rhok1d(nlayer,iq) Qraf(nlraf,iq)=rhok1d(nlayer,iq)/sum(rhok1d(nlayer,1:ntracers)) enddo Mraf(nlraf)=1D0/sum(Qraf(nlraf,1:ntracers)/mol_tr(1:ntracers)) nlay=sum(rhokraf(nlraf,1:ntracers))/Mraf(nlraf)/masseU Praf(nlraf)=P1D(nlayer) do k=1,ndiffions kn=itrans(k) TIraf(nlraf,k)=TI1D(nlayer,k) QIraf(nlraf,k)=Qraf(nlraf,kn) ! mu=mol_tr(kn)*Mraf(nlraf)/(mol_tr(kn)+Mraf(nlraf)) ! nlay=Praf(nlraf)/kbolt/TNraf(nlraf) ! FIraf(nlraf,k)=nlay*a0*Mraf(nlraf)/(Mraf(nlraf)+mol_tr(kn)) ! & *sqrt(apol/mu) enddo Teraf(nlraf)=TE1D(nlayer) QEraf(nlraf)=Qraf(nlraf,ke) apol=0D0 do iq=1,ntracers ! REF: CRC Handbook CHEMISTRY and Physics - 95th EDITION 2014-2015 if (tname(gcmind(iq)).eq.'co2') apol=2.911d0 ! in 10^-24 cm3 (Withers 2009 for CO2) if (tname(gcmind(iq)).eq.'co') apol=1.95d0 ! in 10^-24 cm3 if (tname(gcmind(iq)).eq.'o2') apol=1.5689d0 ! in 10^-24 cm3 if (tname(gcmind(iq)).eq.'n2') apol=1.7403d0 ! in 10^-24 cm3 if (tname(gcmind(iq)).eq.'h2') apol=0.8042d0 ! in 10^-24 cm3 if (tname(gcmind(iq)).eq.'he') apol=0.20550522d0 ! in 10^-24 cm3 if (tname(gcmind(iq)).eq.'h') apol=0.666793d0 ! in 10^-24 cm3 if (tname(gcmind(iq)).eq.'n') apol=1.10d0 ! in 10^-24 cm3 if (tname(gcmind(iq)).eq.'o') apol=0.802d0 ! in 10^-24 cm3 if(apol.ne.0d0) then do l=1,nlraf nlay=Rhokraf(l,iq)/Mraf(l)/masseU do k=1,ndiffions kn=itrans(k) mu=mol_tr(kn)*mol_tr(iq)/(mol_tr(kn)+mol_tr(iq)) FIraf(l,k) = FIraf(l,k) + nlay * a0 & * mol_tr(iq)/(mol_tr(iq)+mol_tr(kn)) * sqrt(apol/mu) enddo enddo endif apol=0d0 enddo END SUBROUTINE GCMGRID_PK(Zraf,Praf,Traf,Mraf,QIraf,RIraf, & TIraf,TEraf,P1di,P1D,Tn1d,M1d,qk1d,qk1dnew,TI1D,TI1Dnew, & TE1D,TE1Dnew,nlraf,ndiffions,nlayer,izv,itrans,etrans) IMPLICIT NONE #include "YOMCST.h" INTEGER :: nlraf,ndiffions,nlayer,l,iP,k,izv INTEGER,dimension(ndiffions) :: itrans INTEGER,dimension(1) :: etrans REAL(kind=kind(1.D0)),dimension(nlayer) :: P1d,Tn1d,M1d REAL(kind=kind(1.D0)),dimension(nlayer+1) :: P1di REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: Qk1d,Qk1dnew REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: TI1d,TI1dnew REAL(kind=kind(1.D0)),dimension(nlayer) :: TE1d,TE1Dnew REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: Rknew REAL(kind=kind(1.D0)),dimension(nlraf) :: Zraf,Praf,Traf,Mraf REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: QIraf,TIraf REAL(kind=kind(1.D0)),dimension(nlraf) :: TEraf REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: RIraf REAL(kind=kind(1.D0)) :: masseU,kbolt,Konst,g REAL(kind=kind(1.D0)) :: facP,rhoNloc masseU=1.e-3/RNAVO kbolt=RKBOL Konst=kbolt/masseU g=RG ! below layer iz_vertplasma no effect due to vertical diffusion do l=1,nlayer if (P1D(l) .ge. Praf(1)) then qk1dnew(l,:)=qk1d(l,:) Ti1dnew(l,:)=Ti1d(l,:) Te1dnew(l)=TE1d(l) endif if (P1D(l) .lt. Praf(1)) then iP=1 do while(Praf(iP) .gt. P1D(l)) iP=iP+1 enddo IP=iP-1 facP=(P1D(l)-Praf(iP))/(Praf(iP+1)-Praf(iP)) rhoNloc=P1D(l)/kbolt/TN1D(l)*M1d(l)*masseU do k=1,ndiffions Rknew(l,k)=RIraf(iP,k)+(RIraf(iP+1,k)-RIraf(iP,k))*facP TI1dnew(l,k)=TIraf(IP,k)+(TIraf(iP+1,k)-TIraf(iP,k))*facP enddo TE1Dnew(l)=TEraf(iP)+(TEraf(iP+1)-TEraf(iP))*facP do k=1,ndiffions Qk1dnew(l,k)=Rknew(l,k)/rhoNloc enddo endif enddo END SUBROUTINE GCMGRID_P2K(Zraf,Praf,Traf,Mraf,QIraf,RIraf, & Tiraf,Teraf,Wk,Ez,P1Di,P1D,Tn1d,M1d,qk1d,qk1dnew, & TI1D,TI1Dnew,TE1D,TE1Dnew,Wk1D,Wk1D2,Wk1D3,Ez1d,ZZ1d, & nlraf,ndiffions,nlayer,izv,itrans,etrans,FacQ,FacTI,FacTe) IMPLICIT NONE #include "YOMCST.h" INTEGER :: nlraf,ndiffions,nlayer,l,iP,k,izv INTEGER :: itrans(ndiffions) INTEGER :: etrans(1) REAL(kind=kind(1.D0)),dimension(nlayer) :: P1d,Tn1d,M1d,ZZ1d REAL(kind=kind(1.D0)),dimension(nlayer+1) :: P1di REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: Qk1d,Qk1dnew REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: TI1d,TI1dnew REAL(kind=kind(1.D0)),dimension(nlayer) :: TE1D,TE1Dnew REAL(kind=kind(1.D0)),dimension(nlayer+1,ndiffions) :: Wk1D,WK1d2 REAL(kind=kind(1.D0)),dimension(nlayer+1,ndiffions) :: Wk1D3 REAL(kind=kind(1.D0)),dimension(nlayer+1) :: Ez1D ! REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: Rknew REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: FacQ,FacTI REAL(kind=kind(1.D0)),dimension(nlayer) :: FacTE REAL(kind=kind(1.D0)),dimension(nlraf) :: Zraf,Praf,Traf REAL(kind=kind(1.D0)),dimension(nlraf) :: Mraf REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: QIraf REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: TIraf REAL(kind=kind(1.D0)),dimension(nlraf) :: TEraf REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: RIraf REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: Wk REAL(kind=kind(1.D0)),dimension(nlraf) :: Ez REAL(kind=kind(1.D0)) masseU,kbolt,Konst,g REAL(kind=kind(1.D0)) facP,rhoNloc masseU=1.e-3/RNAVO kbolt=RKBOL Konst=kbolt/masseU g=RG ! below layer iz_vertplasma no effect due to vertical diffusion do l=1,nlayer if (P1D(l) .ge. Praf(1)) then qk1dnew(l,1:ndiffions)=qk1d(l,1:ndiffions) Ti1dnew(l,1:ndiffions)=Ti1d(l,1:ndiffions) endif if (P1di(l) .ge. Praf(1)) then Ez1D(l)=0D0 Wk1D(l,1:ndiffions)=0D0 Wk1D2(l,1:ndiffions)=0D0 Wk1D3(l,1:ndiffions)=0D0 endif if (P1D(l) .lt. Praf(1)) then iP=1 do while(Praf(iP) .gt. P1D(l)) iP=iP+1 enddo iP=iP-1 ! indP=maxloc(Praf,mask=Praf < P1D(l)) ! IP=indP(1)-1 facP=(P1D(l)-Praf(iP))/(Praf(iP+1)-Praf(iP)) rhoNloc=P1D(l)/kbolt/TN1D(l)*M1d(l)*masseU do k=1,ndiffions Qk1dnew(l,k)=(RIraf(iP,k)+(RIraf(iP+1,k)-RIraf(iP,k))*facP) & *facQ(l,k)/rhoNloc TI1dnew(l,k)=(TIraf(IP,k)+(TIraf(iP+1,k)-TIraf(iP,k))*facP) & *facTI(l,k) enddo TE1Dnew(l)=(TEraf(iP)+(TEraf(iP+1)-TEraf(iP))*facP)*facTE(l) ! if (l .eq. nlayer) then ! print*,'iP',iP,Praf(iP),P1D(l),Praf(iP+1),nlraf,facP ! print*,RIraf(iP,1),RIraf(iP+1,1),facQ(l,1),rhoNloc,facTI(l,1) ! print*,Qk1Dnew(l,1),TI1Dnew(l,1) ! endif endif if (P1Di(l) .lt. Praf(1)) then iP=1 do while(Praf(iP) .gt. P1Di(l)) iP=iP+1 enddo iP=iP-1 facP=(P1Di(l)-Praf(iP))/(Praf(iP+1)-Praf(iP)) do k=1,ndiffions wk1d(l,k)=(Wk(iP,k)+(Wk(iP+1,k)-Wk(iP,k))*facP) & /(ZZ1D(l)-ZZ1D(l-1)) wk1d2(l,k)=(Wk(iP,k)+(Wk(iP+1,k)-Wk(iP,k))*facP) & *(ZZ1D(l)-ZZ1D(l-1)) wk1d3(l,k)=(Wk(iP,k)+(Wk(iP+1,k)-Wk(iP,k))*facP) enddo Ez1D(l)=(EZ(iP)+(Ez(iP+1)-Ez(iP))*facP) endif enddo ! l Wk1d(nlayer+1,1:ndiffions)=0D0 Wk1d2(nlayer+1,1:ndiffions)=0D0 Wk1d3(nlayer+1,1:ndiffions)=0D0 Ez1D(nlayer+1)=0D0 END SUBROUTINE CORRFACK(qbef,qaft,Tbef,Taft,Tebef,Teaf, & Fq,Ft,Fte,nl,nq) IMPLICIT NONE INTEGER :: nl,nq,l,k REAL(kind=kind(1.D0)),dimension(nl,nq) :: qbef,qaft,Fq REAL(kind=kind(1.D0)),dimension(nl,nq) :: Tbef,Taft,Ft REAL(kind=kind(1.D0)),dimension(nl) :: Tebef,Teaf,Fte do l=1,nl do k=1,nq Fq(l,k)=qbef(l,k)/qaft(l,k) Ft(l,k)=Tbef(l,k)/Taft(l,k) enddo Fte(l)=Tebef(l)/Teaf(l) enddo END SUBROUTINE DCOEFFK(TI,NUI,M,DI,nl,nqdyn,nqtot,ind) IMPLICIT NONE INTEGER :: nl,l,nqdyn,nqtot,k,kn INTEGER,dimension(nqdyn) :: ind REAL(kind=kind(1.D0)),dimension(nl,nqdyn) :: NUI,TI,DI REAL,dimension(nqtot) :: M REAL(kind=kind(1.D0)) :: kbolt,masseU masseU=1.660538782d-27 kbolt=1.3806504d-23 do k=1,nqdyn kn=ind(k) do l=1,nl DI(l,k)=kbolt/masseU*TI(l,k)/NUI(l,k)/M(kn) enddo enddo END SUBROUTINE DIFFPARAMK(T,Te,RE,D,Dz, & alpha,beta,gama,delta,eps,dzeta,eta,nl,Wmax) IMPLICIT NONE INTEGER :: nl,l REAL(kind=kind(1.D0)),DIMENSION(nl) :: T,D REAL(kind=kind(1.D0)),DIMENSION(nl) :: Te,Re REAL(kind=kind(1.D0)) :: DZ,Wmax REAL(kind=kind(1.D0)),DIMENSION(nl) :: alpha,beta,gama,delta,eps REAL(kind=kind(1.D0)),DIMENSION(nl) :: dzeta,eta alpha(1)=1D0/T(1)*(-3D0*T(1)+4D0*T(2)-T(3))/(2D0*DZ)+ & 1D0/T(1)*(-3D0*Te(1)+4D0*Te(2)-Te(3))/(2D0*DZ) delta(1)=(-3D0*D(1)+4D0*D(2)-D(3))/(2D0*Dz) beta(1)=1D0/T(1) dzeta(1)=Te(1)/T(1)* & (-3D0*log(Re(1))+4D0*log(Re(2))-log(Re(3)))/(2D0*DZ) ! Add a limit to dzeta (test to avoid unstabilities) ! if (abs(dzeta(1)) .ge. 1d0 .and. Wmax .ge. 500d0) then ! dzeta(1)=1d0*dzeta(1)/abs(dzeta(1)) ! endif ! dzeta(1)=Te(1)/T(1)/Re(1)* ! & (-3D0*Re(1)+4D0*Re(2)-Re(3))/(2D0*Dz) do l=2,nl-1 alpha(l)=1D0/T(l)*(T(l+1)-T(l-1))/(2D0*DZ)+ & 1D0/T(l)*(Te(l+1)-Te(l-1))/(2D0*DZ) delta(l)=(D(l+1)-D(l-1))/(2D0*Dz) beta(l)=1D0/T(l) dzeta(l)=Te(l)/T(l)*(log(Re(l+1))-log(Re(l-1)))/(2D0*DZ) ! dzeta(l)=Te(l)/T(l)/Re(l)*(Re(l+1)-Re(l-1))/(2D0*DZ) ! if (abs(dzeta(l)) .ge. 1d0 .and. Wmax .ge. 500d0) then ! dzeta(l)=1d0*dzeta(l)/abs(dzeta(l)) ! endif enddo ! Upper altitude values alpha(nl)=1D0/T(nl)*(3D0*T(nl)-4D0*T(nl-1)+T(nl-2))/(2D0*DZ)+ & 1D0/T(nl)*(3D0*Te(nl)-4D0*Te(nl-1)+Te(nl-2))/(2D0*DZ) delta(nl)=(3D0*D(nl)-4D0*D(nl-1)+D(nl-2))/(2D0*DZ) beta(nl)=1D0/T(nl) dzeta(nl)=Te(nl)/T(nl)* & (3D0*log(Re(nl))-4D0*log(Re(nl-1))+log(Re(nl-2)))/(2D0*DZ) ! if (abs(dzeta(nl)) .ge. 1d0 .and. Wmax .ge. 500d0) then ! dzeta(nl)=1d0*dzeta(nl)/abs(dzeta(nl)) ! endif ! dzeta(nl)=Te(nl)/T(nl)/Re(nl)* ! & (3D0*Re(nl)-4D0*Re(nl-1)+Re(nl-2))/(2D0*Dz) ! Compute the gama and eps coefficients ! Lower altitude values gama(1)=(-3D0*beta(1) +4D0*beta(2) -beta(3) )/(2d0*dz) eps(1) =(-3D0*alpha(1)+4D0*alpha(2)-alpha(3))/(2d0*dz) eta(1) =(-3D0*dzeta(1)+4D0*dzeta(2)-dzeta(3))/(2D0*dz) do l=2,nl-1 gama(l)=(beta(l+1) -beta(l-1) )/(2d0*Dz) eps(l) =(alpha(l+1)-alpha(l-1))/(2d0*Dz) eta(l) =(dzeta(l+1)-dzeta(l-1))/(2D0*Dz) end do gama(nl)=(3D0*beta(nl) -4D0*beta(nl-1) +beta(nl-2) )/(2D0*DZ) eps(nl) =(3D0*alpha(nl)-4D0*alpha(nl-1)+alpha(nl-2))/(2D0*DZ) eta(nl) =(3D0*dzeta(nl)-4D0*dzeta(nl-1)+dzeta(nl-2))/(2D0*dz) ! print*,'test',alpha(nl-1),beta(nl-1),dzeta(nl-1) END SUBROUTINE MATCOEFFK(alpha,beta,gama,delta,eps,dzeta,eta, & Dad,Rhoad,Facesc,dz,dt,A,B,C,D,Ti,Te,nl) IMPLICIT NONE INTEGER :: nl,l REAL*8,DIMENSION(nl) :: alpha,beta,gama,delta,eps,dzeta,eta REAL*8,DIMENSION(nl) :: Dad,RHoad,Ti,Te REAL*8 :: dz,dt,del1,del2,del3,FacEsc REAL*8,DIMENSION(nl) :: A,B,C,D del1=dt/2d0/dz del2=dt/dz/dz del3=dt ! lower boundary coefficients no change A(1)=0d0 B(1)=1d0 C(1)=0d0 D(1)=rhoAd(1) do l=2,nl-1 A(l)=(delta(l)+(alpha(l)+beta(l)+dzeta(l))*Dad(l))*del1 & -Dad(l)*del2 B(l)=-(delta(l)*(alpha(l)+beta(l)+dzeta(l)) & +Dad(l)*(gama(l)+eps(l)+eta(l)))*del3 B(l)=B(l)+1d0+2d0*Dad(l)*del2 C(l)=-(delta(l)+(alpha(l)+beta(l)+dzeta(l))*Dad(l))*del1 & -Dad(l)*del2 D(l)=rhoAD(l) enddo ! Upper boundary coefficients Diffusion profile C(nl)=0d0 B(nl)=-1d0 A(nl)=exp(-dZ*(alpha(nl-1)+beta(nl-1)+dzeta(nl-1)))*FacEsc ! A(nl)=exp(-dZ*(alpha(nl)+beta(nl)+dzeta(nl)))*FacEsc ! A(nl)=exp(-dZ*(alpha(nl-1)+beta(nl-1)))*FacEsc ! A(nl)=exp(-dZ*(alpha(nl-1)+beta(nl-1))* ! & (Ti(nl-1)/(Te(nl-1)+Ti(nl-1))))*Facesc D(nl)=0D0 END SUBROUTINE CheckmassK(X,Y,nl) IMPLICIT NONE INTEGER :: nl REAL*8,DIMENSION(nl) :: X,Y REAL*8 Xtot,Ytot Xtot=sum(X) Ytot=sum(Y) if (abs((Xtot-Ytot)/Xtot) .gt. 1d-4) then print*,'no conservation for mass',Xtot,Ytot endif END