Index: /trunk/LMDZ.MARS/libf/phymars/callkeys_mod.F90 =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/callkeys_mod.F90 (revision 4062) +++ /trunk/LMDZ.MARS/libf/phymars/callkeys_mod.F90 (revision 4063) @@ -64,8 +64,6 @@ logical,save :: topflows ! entrainment by mountain top dust flows parametrization !$OMP THREADPRIVATE(topflows) -logical,save :: coagulation ! coagulation of dust particles -logical,save :: fullcoag ! coagulation of dust particles using full equations instead of lookup tables -logical,save :: kernel_b, kernel_g, kernel_de, kernel_ti ! kernels for different types of dust coagulation -!$OMP THREADPRIVATE(coagulation,fullcoag,kernel_b,kernel_g,kernel_de,kernel_ti) +logical,save :: dust_coagulation ! coagulation of dust particles +!$OMP THREADPRIVATE(coagulation) integer,save :: coal_kg ! mode for coalescence 0,1, or 2 !$OMP THREADPRIVATE(coal_kg) Index: unk/LMDZ.MARS/libf/phymars/coagulation_mod.F90 =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/coagulation_mod.F90 (revision 4062) +++ (revision ) @@ -1,719 +1,0 @@ -module coagulation_mod - -! ------------------------ -! Dust coagulation scheme -! Based on Bertrand et al., 2022 -! Author: T. Bertrand -! ------------------------ - -use comcstfi_h, only: pi,g,r,mugaz -use tracer_mod, only: igcm_dust_mass,igcm_dust_number, & - igcm_stormdust_mass,igcm_stormdust_number, & - igcm_topdust_mass,igcm_topdust_number -use microphys_h, only: kbz,m0co2 -use callkeys_mod, only: kernel_b,kernel_g,kernel_de,kernel_ti, & - fullcoag, coal_kg -implicit none - -public :: coagul_init - -! ------------------------------- -! Coagulation Input Parameters -! ------------------------------- -real :: effb=1. ! Sticking efficiency for B coagulation -real :: coal_fac=1. ! Overall coalescence factor -! Bin discretization of particle sizes -real*8, save :: r1_coag = 0.01e-6 ! min radius -real*8, save :: rn_coag = 40.e-6 ! max radius -integer,parameter :: nres_coag = 10 ! nb bins -real*8, save :: vrat_coag ! volume ratio -real*8, save, allocatable :: rads_coag(:),vols_coag(:),deltar_coag(:) - -!$OMP THREADPRIVATE(r1_coag,rn_coag,vrat_coag) -!$OMP THREADPRIVATE(rads_coag,vols_coag,deltar_coag) - -! parameters of the look up tables of coagulation coefficients -integer,parameter :: table_numt = 15 ! Nb of temperature bins -integer,parameter :: table_nump = 25 ! Nb of pressure bins -integer,parameter :: table_numm = 25 ! Nb of mass bins -real, save :: table_pt(table_numt) -real, save :: table_pres(table_nump) -real, save :: table_mfp(table_numm) -real, save :: table_b(table_numt,table_nump,nres_coag,nres_coag) ! Brownian -real, save :: table_g(table_numt,table_nump,nres_coag,nres_coag) ! Gravitational -real, save :: table_de(table_numt,table_nump,nres_coag,nres_coag) ! Brownian diffusion enhancement -real, save :: table_ti(table_numt,table_nump,nres_coag,nres_coag) ! Turbulent -real*8, public :: dev_dt = 0.63676 ! standard deviation of dust distribution -real, public ::rho_dust=2500. ! Mars dust density (kg.m-3) - -!$OMP THREADPRIVATE(table_pt,table_pres,table_mfp) -!$OMP THREADPRIVATE(table_b,table_g,table_de,table_ti) - -contains - -!####################################################################### - -subroutine coagul_init() - -! Initialisation of the particle bin discretization - -integer :: i -allocate(rads_coag(nres_coag)) -allocate(vols_coag(nres_coag)) -allocate(deltar_coag(nres_coag)) - -vrat_coag=(rn_coag/r1_coag)**(3./(nres_coag-1)) - -do i = 1, nres_coag - rads_coag(i) = r1_coag*vrat_coag**((i-1)/3.) - vols_coag(i) = 4./3.*pi*r1_coag**3*vrat_coag**(i-1) -enddo -! diameter width -deltar_coag(:)=2.*rads_coag(:)*2.**(1/3.)*(vrat_coag**(1./3.)-1)/(1+vrat_coag)**(1./3.) - -! Making Lookup Tables -if (.not.fullcoag) then - call make_tables() -endif -!deallocate(rads_coag) -!deallocate(vols_coag) -!deallocate(deltar_coag) - -end subroutine coagul_init - -!####################################################################### - -SUBROUTINE coagul_main( & - ngrid, nlayer, nq, ptime, ptimestep, pq, pdqfi, pt, pdtfi, & - pplay, pplev, pdqcoag) - -!-------------------------------------------------------- -! Input Variables -!-------------------------------------------------------- -INTEGER, INTENT(IN) :: ngrid ! number of horizontal grid points -INTEGER, INTENT(IN) :: nlayer ! number of vertical grid points -INTEGER, INTENT(IN) :: nq ! number of tracer species -REAL, INTENT(IN) :: ptime -REAL, INTENT(IN) :: ptimestep - -REAL, INTENT(IN) :: pq(ngrid,nlayer,nq) ! advected field nq -REAL, INTENT(IN) :: pdqfi(ngrid,nlayer,nq)! tendancy field pq -REAL, INTENT(IN) :: pt(ngrid,nlayer) ! temperature at mid-layer (K) -REAL, INTENT(IN) :: pdtfi(ngrid,nlayer) ! tendancy temperature at mid-layer (K) -REAL, INTENT(IN) :: pplay(ngrid,nlayer) ! pressure at middle of the layers -REAL, INTENT(IN) :: pplev(ngrid,nlayer+1) ! pressure at intermediate levels - -!-------------------------------------------------------- -! Output Variables -!-------------------------------------------------------- -REAL, INTENT(OUT) :: pdqcoag(ngrid,nlayer,nq) ! tendancy field for dust coagulation - -!----------------------------------------------------------------------- -! Local variables -!----------------------------------------------------------------------- -integer :: ie, je, id, jd, iq, i, j, k, l, ii, jj -real dens,dev,cst,sig0 -real mfp0,pt0,rad0,rad10,rad20,rho0,pres0 -real term1,term2,kernel,norm -real, dimension(ngrid, nlayer) :: r0, rho, rn, ntot, mtot, mtotnew -real, dimension(ngrid, nlayer, nq) :: numb_new, numb_ini, numb_ini_dis, mass_ini, mass_new -real, dimension(ngrid, nlayer, nres_coag) :: ndis, ndis_new -real, dimension(ngrid,nlayer,nres_coag,nq) :: ndis_tab,rat_tab,ndis_tab_new -real, dimension(ngrid, nlayer) :: zt -real, dimension(ngrid, nlayer, nq) :: zq - -logical, save :: firstcall = .true. -!$OMP THREADPRIVATE(firstcall) -!for boxinterp function: -integer :: t1,t2,p1,p2,m1,m2 - -!====================================================================== -! Main routine - -!----------------------------------------------------------------------- -!*** 1) Initialisations -!----------------------------------------------------------------------- - -!!! Particles properties -dens = rho_dust ! density dust -sig0 = dev_dt ! Standard deviation of the dust distribution -cst = .75 / (pi*dens) * exp( -4.5*sig0**2. ) ! fixed distribution - -!!! Radius and Volumes defined in coagul_init -! rads_coag, vols_coag, deltar_coag, vrat_coag - -!!! Getting tracers and temperature fields -zq(:,:,:)=pq(:,:,:)+pdqfi(:,:,:)*ptimestep -zt(:,:)=pt(:,:)+pdtfi(:,:)*ptimestep - -!!! Atmospheric state -DO l=1,nlayer - rho(:,l)=pplay(:,l)/(r*pt(:,l)) -ENDDO - -!!! Initializing tracers and other fields -pdqcoag(:,:,:)=0.d0 -ndis_tab(:,:,:,:)=0.d0 -ndis_new(:,:,:)=0.d0 -mtot(:,:)=0.d0 -rat_tab(:,:,:,:) = 0.d0 -numb_ini_dis(:,:,:) = 0.d0 -numb_ini(:,:,:) = 0.d0 -numb_new(:,:,:) = 0.d0 -mass_ini(:,:,:) = 0.d0 -mass_new(:,:,:) = 0.d0 - -!!! Preparing the initial log normal distribution -do iq=1,nq - if ((iq.eq.igcm_dust_mass).or.(iq.eq.igcm_stormdust_mass).or.(iq.eq.igcm_topdust_mass)) then - mass_ini(:,:,iq)=max(zq(:,:,iq),1e-15) - numb_ini(:,:,iq)=max(zq(:,:,iq-1),1e-15) - - !!! Get total mass to ensure mass conservation - mtot(:,:)=mtot(:,:)+mass_ini(:,:,iq) - - !!! Initial Lognormal distribution - ! sig0 remains constant - r0(:,:)=(3/4.*mass_ini(:,:,iq)/(pi*numb_ini(:,:,iq)*dens))**(1/3.)*exp(-1.5*sig0**2) - - do i=1,nres_coag ! Sum of the ndis for all dust modes - ndis_tab(:,:,i,iq)=numb_ini(:,:,iq)*rho(:,:)*deltar_coag(i)/(2.*rads_coag(i)*(2*pi)**0.5*sig0)*exp(-0.5*(log(rads_coag(i)/r0(:,:)))**2/sig0**2) - enddo - !!! Initial number mixing ratio from bin discretization nb/kg - numb_ini_dis(:,:,iq)=sum(ndis_tab(:,:,:,iq),3)/rho(:,:) - numb_new(:,:,iq)=numb_ini_dis(:,:,iq) - endif -enddo - - -!!! Defining initial number ratio for each mode -do iq=1,nq - if ((iq.eq.igcm_dust_mass).or.(iq.eq.igcm_stormdust_mass).or.(iq.eq.igcm_topdust_mass)) then - rat_tab(:,:,:,iq)=ndis_tab(:,:,:,iq)/sum(ndis_tab(:,:,:,:),4) - endif -enddo - -!!! Normalizing the ratio (make sure sum=1) -do iq=1,nq - if ((iq.eq.igcm_dust_mass).or.(iq.eq.igcm_stormdust_mass).or.(iq.eq.igcm_topdust_mass)) then - rat_tab(:,:,:,iq)=rat_tab(:,:,:,iq)/sum(rat_tab(:,:,:,:),4) - endif -enddo - -!!! Total initial distribution nb/m3 per bin -ndis(:,:,:)=sum(ndis_tab(:,:,:,:),4) -!!! Total initial distribution nb/kg -ntot(:,:)=sum(ndis(:,:,:),3)/rho(:,:) - -!----------------------------------------------------------------------- -!*** 2) Coagulation -!----------------------------------------------------------------------- -! Spatial loop => 1D - -if (fullcoag) then ! computing full coagulation equations - do i=1,ngrid - do j=1,nlayer - if (ntot(i,j).gt.1000.) then ! activating coagulation above this threshold - pt0=pt(i,j) - rho0=rho(i,j) - mfp0=mfp(pt0,rho0) - - ! Bin loop - do k=1,nres_coag - - ! Term 1 - term1=0. - do jj=1,k - do ii=1,k-1 - kernel=0. - if (kernel_b) kernel=kernel+betab(pt0,mfp0,rads_coag(ii),rads_coag(jj)) - if (kernel_g) kernel=kernel+betag(pt0,mfp0,rads_coag(ii),rads_coag(jj)) - if (kernel_de) kernel=kernel+betade(pt0,mfp0,rho0,rads_coag(ii),rads_coag(jj)) - if (kernel_ti) kernel=kernel+betati(pt0,mfp0,rho0,rads_coag(ii),rads_coag(jj)) - term1=term1+frac(ii,jj,k)*vols_coag(ii)*kernel*ndis_new(i,j,ii)*ndis(i,j,jj) - enddo - enddo - term1=term1*ptimestep/vols_coag(k) - - ! Term 2 - term2=0. - do jj=1,nres_coag - kernel=0. - if (kernel_b) kernel=kernel+betab(pt0,mfp0,rads_coag(k),rads_coag(jj)) - if (kernel_g) kernel=kernel+betag(pt0,mfp0,rads_coag(k),rads_coag(jj)) - if (kernel_de) kernel=kernel+betade(pt0,mfp0,rho0,rads_coag(k),rads_coag(jj)) - if (kernel_ti) kernel=kernel+betati(pt0,mfp0,rho0,rads_coag(k),rads_coag(jj)) - term2=term2+(1-frac(k,jj,k))*kernel*ndis(i,j,jj) - enddo - term2=term2*ptimestep - - ndis_new(i,j,k)=(ndis(i,j,k)+coal_fac*term1)/(1.+coal_fac*term2) - - enddo - - norm=sum(ndis(i,j,:)*vols_coag(:))/sum(ndis_new(i,j,:)*vols_coag(:)) - ndis_new(i,j,:)=ndis_new(i,j,:)*norm - - endif - enddo - enddo - -else ! fullcoag - - ! Using Lookup tables - do i=1,ngrid - do j=1,nlayer - if (ntot(i,j).gt.1000.) then ! activating coagulation above this threshold - pt0=pt(i,j) - rho0=rho(i,j) - mfp0=mfp(pt0,rho0) - pres0=pplay(i,j) - !find the corners for interpolation - t2=findval(table_pt,pt0) - if (t2==0) t2=table_numt - if (t2==1) t2=2 - t1=t2-1 - m2=findval(table_mfp,mfp0) - if (m2==0) m2=table_numm - if (m2==1) m2=2 - m1=m2-1 - p2=findval(table_pres,pres0) - if (p2==0) p2=table_numt - if (p2==1) p2=2 - p1=p2-1 - - ! Bin loop - do k=1,nres_coag - - ! Term 1 - term1=0. - do jj=1,k - do ii=1,k-1 - kernel=0. - if (kernel_b) kernel=kernel+boxinterp(table_b(t1,m1,ii,jj), & - table_b(t1,m2,ii,jj), & - table_b(t2,m1,ii,jj), & - table_b(t2,m2,ii,jj), & - table_pt(t1),table_pt(t2),log10(table_mfp(m1)),log10(table_mfp(m2)),pt0,log10(mfp0)) - if (kernel_g) kernel=kernel+boxinterp(table_g(t1,m1,ii,jj), & - table_g(t1,m2,ii,jj), & - table_g(t2,m1,ii,jj), & - table_g(t2,m2,ii,jj), & - table_pt(t1),table_pt(t2),log10(table_mfp(m1)),log10(table_mfp(m2)),pt0,log10(mfp0)) - if (kernel_de) kernel=kernel+boxinterp(table_de(t1,p1,ii,jj), & - table_de(t1,p2,ii,jj), & - table_de(t2,p1,ii,jj), & - table_de(t2,p2,ii,jj), & - table_pt(t1),table_pt(t2),log10(table_pres(p1)),log10(table_pres(p2)),pt0,log10(pres0)) - ! if (kernel_ti) kernel=kernel+betati(pt0,mfp0,rho(i,j,l),ii,jj) !not implemented yet - term1=term1+frac(ii,jj,k)*vols_coag(ii)*kernel*ndis_new(i,j,ii)*ndis(i,j,jj) - enddo - enddo - term1=term1*ptimestep/vols_coag(k) - - - ! Term 2 - term2=0. - do jj=1,nres_coag - kernel=0. - if (kernel_b) kernel=kernel+boxinterp(table_b(t1,m1,k,jj), & - table_b(t1,m2,k,jj), & - table_b(t2,m1,k,jj), & - table_b(t2,m2,k,jj), & - table_pt(t1),table_pt(t2),log10(table_mfp(m1)),log10(table_mfp(m2)),pt0,log10(mfp0)) - if (kernel_g) kernel=kernel+boxinterp(table_g(t1,m1,k,jj), & - table_g(t1,m2,k,jj), & - table_g(t2,m1,k,jj), & - table_g(t2,m2,k,jj), & - table_pt(t1),table_pt(t2),log10(table_mfp(m1)),log10(table_mfp(m2)),pt0,log10(mfp0)) - if (kernel_de) kernel=kernel+boxinterp(table_de(t1,p1,k,jj), & - table_de(t1,p2,k,jj), & - table_de(t2,p1,k,jj), & - table_de(t2,p2,k,jj), & - table_pt(t1),table_pt(t2),log10(table_pres(p1)),log10(table_pres(p2)),pt0,log10(pres0)) - ! if (kernel_ti) kernel=kernel+betati(pt0,mfp0,rho(i,j,l),k,jj) !not implemented yet - term2=term2+(1-frac(k,jj,k))*kernel*ndis(i,j,jj) - enddo - term2=term2*ptimestep - - ndis_new(i,j,k)=(ndis(i,j,k)+coal_fac*term1)/(1.+coal_fac*term2) - enddo - - !! Volume conservation - norm=sum(ndis(i,j,:)*vols_coag(:))/sum(ndis_new(i,j,:)*vols_coag(:)) - ndis_new(i,j,:)=ndis_new(i,j,:)*norm - - endif - - enddo - enddo -endif ! fullcoag - -!----------------------------------------------------------------------- -!*** 3) new moments mass and number -!----------------------------------------------------------------------- -mtotnew(:,:)=1.d-15 - -do iq=1,nq - if ((iq.eq.igcm_dust_mass).or.(iq.eq.igcm_stormdust_mass).or.(iq.eq.igcm_topdust_mass)) then - !! Estimate new distribution in each dust mode using the initial number ratios - ndis_tab_new(:,:,:,iq)=ndis_new(:,:,:)*rat_tab(:,:,:,iq) - - do i=1,ngrid - do j=1,nlayer - !! New number and mass mixing ratio in each mode - numb_new(i,j,iq)=sum(ndis_tab_new(i,j,:,iq))/rho(i,j) - mass_new(i,j,iq)=sum(ndis_tab_new(i,j,:,iq)*vols_coag(:))*dens/rho(i,j) - enddo - enddo - - !! New total mass - mtotnew(:,:)=mtotnew(:,:)+mass_new(:,:,iq) - endif -enddo - -do iq=1,nq - if ((iq.eq.igcm_dust_mass).or.(iq.eq.igcm_stormdust_mass).or.(iq.eq.igcm_topdust_mass)) then - !! mass conservation - mass_new(:,:,iq)=mass_new(:,:,iq)/mtotnew(:,:)*mtot(:,:) - !! New tendancy - pdqcoag(:,:,iq)=(mass_new(:,:,iq)-mass_ini(:,:,iq))/ptimestep - !pdqcoag(:,:,iq-1)=(numb_new(:,:,iq)-numb_ini(:,:,iq))/ptimestep - pdqcoag(:,:,iq-1)=(numb_new(:,:,iq)-numb_ini_dis(:,:,iq))/ptimestep - !! check - where (zq(:,:,iq)+pdqcoag(:,:,iq)*ptimestep.le.0.) - pdqcoag(:,:,iq)=-zq(:,:,iq)/ptimestep+1.e-15 - end where - where (zq(:,:,iq-1)+pdqcoag(:,:,iq-1)*ptimestep.le.0.) - pdqcoag(:,:,iq-1)=-zq(:,:,iq-1)/ptimestep+1.e-14 - end where - endif -enddo - -end subroutine coagul_main - - -!====================================================================== -!====================================================================== -real function dynvis(pt) - implicit none - ! dynamic viscosity following sutherland's formula for CO2 - ! Pa.s = kg m-1 s-1 - real, intent(in) :: pt - dynvis=1.37e-5*(273.15 + 222.)/(pt + 222.)*(pt/273.15)**(3./2.) -end function dynvis - -real function kinvis(pt,rho) - implicit none - real, intent(in) :: pt,rho - kinvis=dynvis(pt)/rho -end function kinvis - -real function thervel(pt,massm) - implicit none - ! thermal velocity of an air molecule or a dust particle - ! m s-1 - real, intent(in) :: pt, massm - thervel=(8*kbz*pt/(pi*massm))**0.5 -end function thervel - -real function mfp(pt,rho) - implicit none - ! thermal velocity of an air molecule - ! m s-1 - real, intent(in) :: pt, rho - mfp=2.*kinvis(pt,rho)/thervel(pt,m0co2) ! mean free path (m) -end function mfp - -real function cun(pt,mfp,rad) - implicit none - !! Cunningham slip flow correction - real, intent(in) :: pt, mfp, rad - cun=1.+mfp/rad*(1.246+0.42*exp(-0.87/(mfp/rad))) -end function cun - -real function diff(pt,mfp,rad) - implicit none - !! diffusion coefficient - real, intent(in) :: pt, mfp, rad - diff=kbz*pt/(6.*pi*rad*dynvis(pt))*cun(pt,mfp,rad) -end function diff - -real function lambda(pt, mfp, rad) - implicit none - !! particle mean free path - real, intent(in) :: pt,mfp,rad - real massm,dens - dens = rho_dust ! density dust - massm=dens*4./3.*pi*rad**3 - lambda=2.*diff(pt, mfp, rad)/(pi*thervel(pt,massm)) -end function lambda - -real function delta(pt,mfp,rad) - implicit none - !! mean distance - real, intent(in) :: pt,mfp,rad - real :: ltmp - - ltmp=lambda(pt,mfp,rad) - delta=( (2.*rad+ltmp)**3 - (4.*rad**2+ltmp**2)**1.5 ) /(6*rad*ltmp) - 2*rad -end function delta - -real function fallv(pt,mfp,rad) - implicit none - !! fall velocity - real, intent(in) :: pt,mfp,rad - real dens - dens = rho_dust ! density dust - fallv=2./9.*rad**2*dens*g/dynvis(pt)*cun(pt,mfp,rad) -end function fallv - -real function reyn(pt,mfp,rad,rho) - implicit none - !! Reynold number - real, intent(in) :: pt,rad,mfp,rho - reyn=2.*rad*fallv(pt,mfp,rad)/kinvis(pt,rho) -end function reyn - -real function stokes(pt,mfp,rad1,rad2) - implicit none - !! Stokes number - real, intent(in) :: pt,mfp,rad1,rad2 - real :: f1 - - f1=fallv(pt,mfp,rad1) - stokes=f1*abs(fallv(pt,mfp,rad2)-f1)/(rad2*g) -end function stokes - -real function schmidt(pt,mfp,rho,rad) - implicit none - !! Schmidt number - real, intent(in) :: pt,mfp,rho,rad - schmidt=kinvis(pt,rho)/diff(pt,mfp,rad) -end function schmidt - -real function coal(pt,mfp,rad1,rad2,rho) - implicit none - !! Coalescence efficiency - real, intent(in) :: pt,mfp,rad1,rad2,rho - real :: coalea,coalev - real :: stmp,rtmp - - stmp=stokes(pt,mfp,rad1,rad2) - rtmp=reyn(pt,mfp,rad2,rho) - coalea=stmp**2/(stmp+0.5)**2 - if (stmp>1.214) then - coalev=(1+0.75*log(2.*stmp)/(stmp-1.214))**(-2) - else - coalev=0. - endif - coal=(60.*coalev+coalea*rtmp)/(60.+rtmp) -end function coal - -real function betab(pt,mfp,rad1,rad2) - implicit none - real, intent(in) :: pt,mfp,rad1,rad2 - real :: m1,m2,num,den1,den2 - real :: dens - real :: d1,d2 - - dens = rho_dust ! density dust - d1=diff(pt,mfp,rad1) - d2=diff(pt,mfp,rad2) - !! kernel transition regime - ! Mass of 1 particle 1 and 2 - m1=dens*4./3.*pi*rad1**3 - m2=dens*4./3.*pi*rad2**3 - num=4.*pi*(rad1+rad2)*(d1+d2) - den1=(rad1+rad2)/(rad1+rad2+(delta(pt,mfp,rad1)**2+delta(pt,mfp,rad2)**2)**0.5) - den2=4./effb*(d1+d2)/((thervel(pt,m1)**2+thervel(pt,m2)**2)**0.5*(rad1+rad2)) - betab=num/(den1+den2) -end function betab - -real function betag(pt,mfp,rad1,rad2) - implicit none - real, intent(in) :: pt,mfp,rad1,rad2 - real w1,w2,ee - real :: radt - !! kernel gravitation - ! fall vel - w1=fallv(pt,mfp,rad1) - w2=fallv(pt,mfp,rad2) - radt=(rad1+rad2)**2 - if (coal_kg.eq.0) then ! E=1 - ee=1. - elseif (coal_kg.eq.1) then - ee=1.5*(min(rad1,rad2))**2/radt - elseif (coal_kg.eq.2) then - ee=0.25*(min(rad1,rad2))**2/radt - endif - !b=coal(pt,r1,r2)*pi*(r1+r2)**2*abs(w1-w2) - betag=ee*pi*radt*abs(w1-w2) -end function betag - -real function betade(pt,mfp,rho,rad1,rad2) - implicit none - real, intent(in) :: pt,mfp,rad1,rad2,rho - real rd,sc - !! kernel diffusion enhancement - rd=reyn(pt,mfp,max(rad1,rad2),rho) - sc=schmidt(pt,mfp,rho,min(rad1,rad2)) - if (rd.le.1.) then - betade=betab(pt,mfp,rad1,rad2)*0.45*rd**(1/3.)*sc**(1/3.) - else - betade=betab(pt,mfp,rad1,rad2)*0.45*rd**(1/2.)*sc**(1/3.) - endif -end function betade - -real function betade_tabfunc(pt,pres,rad1,rad2) - implicit none - real, intent(in) :: pt,pres,rad1,rad2 - real rd,sc,rho,mf - !! kernel diffusion enhancement - rho=pres/(pt*r) - mf=mfp(pt,rho) - rd=reyn(pt,mf,max(rad1,rad2),rho) - sc=schmidt(pt,mf,rho,min(rad1,rad2)) - if (rd.le.1.) then - betade_tabfunc=betab(pt,mf,rad1,rad2)*0.45*rd**(1/3.)*sc**(1/3.) - else - betade_tabfunc=betab(pt,mf,rad1,rad2)*0.45*rd**(1/2.)*sc**(1/3.) - endif -end function betade_tabfunc - -real function betati(pt,mfp,rho,rad1,rad2) - implicit none - real, intent(in) :: pt,mfp,rad1,rad2,rho - real w1,w2,eps - ! fall vel - w1=fallv(pt,mfp,rad1) - w2=fallv(pt,mfp,rad2) - betati=eps**(0.75)*pi/(g*kinvis(pt,rho)**(0.25))*(rad1+rad2)**2*abs(w1-w2) -end function betati - -real function betati_tabfunc(pt,pres,rad1,rad2) - implicit none - real, intent(in) :: pt,rad1,rad2 - real w1,w2,eps,mf,rho,pres - rho=pres/(pt*r) - mf=mfp(pt,rho) - ! fall vel - w1=fallv(pt,mf,rad1) - w2=fallv(pt,mf,rad2) - betati_tabfunc=eps**(0.75)*pi/(g*kinvis(pt,rho)**(0.25))*(rad1+rad2)**2*abs(w1-w2) -end function betati_tabfunc - -real function frac(i,j,k) - implicit none - integer, intent(in) :: i,j,k - real vint - ! intermediate volume - vint=vols_coag(i)+vols_coag(j) - if (k.lt.nres_coag) then - if ( (vint.ge.vols_coag(k)) .and. (vint.lt.vols_coag(k+1))) then - frac=(vols_coag(k+1)-vint)/(vols_coag(k+1)-vols_coag(k))*vols_coag(k)/vint - return - endif - endif - if ((k.eq.nres_coag) .and. (vint.ge.vols_coag(k))) then - frac=1. - return - endif - if (k.gt.1) then - if ( (vint.lt.vols_coag(k)) .and. (vint.gt.vols_coag(k-1))) then - frac=(-vols_coag(k-1)+vint)/(-vols_coag(k-1)+vols_coag(k))*vols_coag(k)/vint - return - endif - endif - - frac=0. - return -end function frac - -! ********************************************************* -! Building lookup tables : called in the firstcall -! ********************************************************* - -subroutine make_tables() - -implicit none -integer :: i,j,k,l -real :: t1,t2,p1,p2 - -do i = 1,table_numt - table_pt(i)=50.+25.*real(i) -end do - -do j = 1,table_nump - table_pres(j)=(1.e-8)*10.**(real(j)/2.) -end do - -do j = 1,table_numm - table_mfp(j)=(1.e-8)*10.**(real(j)/2.) -end do - -do i = 1,table_numt - do j = 1,table_numm - do k = 1,nres_coag - do l = 1,nres_coag - table_b(i,j,k,l)=betab(table_pt(i),table_mfp(j),rads_coag(k),rads_coag(l)) - table_g(i,j,k,l)=betag(table_pt(i),table_mfp(j),rads_coag(k),rads_coag(l)) - end do - end do - end do -end do - -do i = 1,table_numt - do j = 1,table_nump - do k = 1,nres_coag - do l = 1,nres_coag - table_de(i,j,k,l)=betade_tabfunc(table_pt(i),table_pres(j),rads_coag(k),rads_coag(l)) - end do - end do - end do -end do - -end subroutine make_tables - -!======================================================================C -!======================================================================C - -real function boxinterp(ft1p1,ft2p1,ft1p2,ft2p2,tref1,tref2,pref1, & - pref2,tmid,pmid) -! -! Calculate 2d interpolation from pressure and temperature -!----------------------------------------------------------------------! -! T2 .FT2P1 .FT2P2 -! -! -! T1 .FT1P1 .FT1P2 -! P1 P2 -!----------------------------------------------------------------------! - -implicit none - -real :: ft1p1, ft2p1, ft1p2, ft2p2, tref1, tref2 -real :: pref1, pref2, tmid, pmid, ans1, ans2, ans - -!======================================================================C - -ans1 = ft1p1 + (ft2p1 - ft1p1)*(tmid - tref1)/(tref2 - tref1) -ans2 = ft1p2 + (ft2p2 - ft1p2)*(tmid - tref1)/(tref2 - tref1) -boxinterp = ans1 + (ans2 - ans1)*(pmid - pref1)/(pref2 - pref1) - -end function boxinterp - -!======================================================================C -!======================================================================C - -integer function findval(array,value) -implicit none -real :: array(:) -real :: value - -findval=minloc(array,dim=1,mask = (array > value)) - -end function findval - -end module coagulation_mod - Index: /trunk/LMDZ.MARS/libf/phymars/conf_phys.F =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/conf_phys.F (revision 4062) +++ /trunk/LMDZ.MARS/libf/phymars/conf_phys.F (revision 4063) @@ -46,4 +46,7 @@ & ads_massive_ice use nonoro_gwd_mix_mod, only: calljliu_gwimix + use dust_coagulation_mod, only: dust_coag_kernel_b, + & dust_coag_kernel_g,dust_coag_kernel_de,dust_coag_kernel_ti, + & coal_coeff_mode,full_coag_equations use lwxn_mod, only: lwxn_linear, lwxn_alphan, lwxn_ncouche use callkeys_mod, only: startphy_file, activice, activeco2ice, @@ -72,7 +75,5 @@ & tf_injection, ti_injection, thermochem, & thermoswater, tituscap, tke_heat_flux, - & topflows, water, coagulation, kernel_b, - & kernel_g,kernel_de,kernel_ti,coal_kg, - & fullcoag + & topflows, water, dust_coagulation use write_output_mod, only: output_diagfi @@ -389,17 +390,19 @@ ! dust particle coagulation write(*,*)"call coagulation of dust" - coagulation=.false. ! default value - call getin_p("coagulation",coagulation) - write(*,*)" coagulation = ",coagulation - fullcoag=.false. ! default value - kernel_b=.true. ! default value (if coagulation=true) - kernel_g=.false. ! default value - kernel_de=.false. ! default value - kernel_ti=.false. ! default value - coal_kg=0 ! default value - if (coagulation) then - write(*,*)" coagulation fullcoag and mode= ",fullcoag,coal_kg - write(*,*)" coagulation kernels= ",kernel_b,kernel_g, - & kernel_de,kernel_ti + dust_coagulation=.false. ! default value + call getin_p("dust_coagulation",dust_coagulation) + write(*,*)" dust_coagulation = ",dust_coagulation + full_coag_equations=.false. ! default value + dust_coag_kernel_b=.true. ! default value (if coagulation=true) + dust_coag_kernel_g=.false. ! default value + dust_coag_kernel_de=.false. ! default value + dust_coag_kernel_ti=.false. ! default value + coal_coeff_mode=0 ! default value + if (dust_coagulation) then + write(*,*)" coagulation full_coag_equations: ", + & full_coag_equations + write(*,*)" coagulation coal_coeff_mode: ",coal_coeff_mode + write(*,*)" coagulation kernels= ",dust_coag_kernel_b, + & dust_coag_kernel_g,dust_coag_kernel_de,dust_coag_kernel_ti endif Index: /trunk/LMDZ.MARS/libf/phymars/dust_coagulation_mod.F90 =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/dust_coagulation_mod.F90 (revision 4063) +++ /trunk/LMDZ.MARS/libf/phymars/dust_coagulation_mod.F90 (revision 4063) @@ -0,0 +1,726 @@ +module dust_coagulation_mod + +! ------------------------ +! Dust coagulation scheme +! Based on Bertrand et al., 2022 +! Author: T. Bertrand +! ------------------------ + +use comcstfi_h, only: pi,g,r,mugaz +use tracer_mod, only: igcm_dust_mass,igcm_dust_number, & + igcm_stormdust_mass,igcm_stormdust_number, & + igcm_topdust_mass,igcm_topdust_number +use microphys_h, only: kbz,m0co2 +implicit none + +public :: dust_coagulation_init + +logical,save :: full_coag_equations ! coagulation of dust particles using full equations instead of lookup tables +logical,save :: dust_coag_kernel_b ! flag to activate Brownian coagulation +logical,save :: dust_coag_kernel_g ! gravitational +logical,save :: dust_coag_kernel_de ! diffusion enhancement +logical,save :: dust_coag_kernel_ti ! Not implemented yet. +!$OMP THREADPRIVATE(full_coag_equations,dust_coag_kernel_b,dust_coag_kernel_g,dust_coag_kernel_de,dust_coag_kernel_ti,coal_coeff_mode) +integer,save :: coal_coeff_mode ! Flag to choose how to calculate coalescence coefficient +!$OMP THREADPRIVATE(coal_coeff_mode) + +! ------------------------------- +! Coagulation Input Parameters +! ------------------------------- +real :: effb=1. ! Sticking efficiency for B coagulation +real :: coal_fac=1. ! Overall coalescence factor +! Bin discretization of particle sizes +real*8, save :: r1_coag = 0.01e-6 ! min radius +real*8, save :: rn_coag = 40.e-6 ! max radius +integer,parameter :: nres_coag = 10 ! nb bins +real*8, save :: vrat_coag ! volume ratio +real*8, save, allocatable :: rads_coag(:),vols_coag(:),deltar_coag(:) + +!$OMP THREADPRIVATE(r1_coag,rn_coag,vrat_coag) +!$OMP THREADPRIVATE(rads_coag,vols_coag,deltar_coag) + +! parameters of the look up tables of coagulation coefficients +integer,parameter :: table_numt = 15 ! Nb of temperature bins +integer,parameter :: table_nump = 25 ! Nb of pressure bins +integer,parameter :: table_numm = 25 ! Nb of mass bins +real, save :: table_pt(table_numt) +real, save :: table_pres(table_nump) +real, save :: table_mfp(table_numm) +real, save :: table_b(table_numt,table_nump,nres_coag,nres_coag) ! Brownian +real, save :: table_g(table_numt,table_nump,nres_coag,nres_coag) ! Gravitational +real, save :: table_de(table_numt,table_nump,nres_coag,nres_coag) ! Brownian diffusion enhancement +real, save :: table_ti(table_numt,table_nump,nres_coag,nres_coag) ! Turbulent +real*8, public :: dev_dt = 0.63676 ! standard deviation of dust distribution +real, public ::rho_dust=2500. ! Mars dust density (kg.m-3) + +!$OMP THREADPRIVATE(table_pt,table_pres,table_mfp) +!$OMP THREADPRIVATE(table_b,table_g,table_de,table_ti) + +contains + +!####################################################################### + +subroutine dust_coagulation_init() + +! Initialisation of the particle bin discretization + +integer :: i +allocate(rads_coag(nres_coag)) +allocate(vols_coag(nres_coag)) +allocate(deltar_coag(nres_coag)) + +vrat_coag=(rn_coag/r1_coag)**(3./(nres_coag-1)) + +do i = 1, nres_coag + rads_coag(i) = r1_coag*vrat_coag**((i-1)/3.) + vols_coag(i) = 4./3.*pi*r1_coag**3*vrat_coag**(i-1) +enddo +! diameter width +deltar_coag(:)=2.*rads_coag(:)*2.**(1/3.)*(vrat_coag**(1./3.)-1)/(1+vrat_coag)**(1./3.) + +! Making Lookup Tables +if (.not.full_coag_equations) then + call make_tables() +endif +!deallocate(rads_coag) +!deallocate(vols_coag) +!deallocate(deltar_coag) + +end subroutine dust_coagulation_init + +!####################################################################### + +SUBROUTINE dust_coagulation_main( & + ngrid, nlayer, nq, ptime, ptimestep, pq, pdqfi, pt, pdtfi, & + pplay, pplev, pdqcoag) + +!-------------------------------------------------------- +! Input Variables +!-------------------------------------------------------- +INTEGER, INTENT(IN) :: ngrid ! number of horizontal grid points +INTEGER, INTENT(IN) :: nlayer ! number of vertical grid points +INTEGER, INTENT(IN) :: nq ! number of tracer species +REAL, INTENT(IN) :: ptime +REAL, INTENT(IN) :: ptimestep + +REAL, INTENT(IN) :: pq(ngrid,nlayer,nq) ! advected field nq +REAL, INTENT(IN) :: pdqfi(ngrid,nlayer,nq)! tendancy field pq +REAL, INTENT(IN) :: pt(ngrid,nlayer) ! temperature at mid-layer (K) +REAL, INTENT(IN) :: pdtfi(ngrid,nlayer) ! tendancy temperature at mid-layer (K) +REAL, INTENT(IN) :: pplay(ngrid,nlayer) ! pressure at middle of the layers +REAL, INTENT(IN) :: pplev(ngrid,nlayer+1) ! pressure at intermediate levels + +!-------------------------------------------------------- +! Output Variables +!-------------------------------------------------------- +REAL, INTENT(OUT) :: pdqcoag(ngrid,nlayer,nq) ! tendancy field for dust coagulation + +!----------------------------------------------------------------------- +! Local variables +!----------------------------------------------------------------------- +integer :: ie, je, id, jd, iq, i, j, k, l, ii, jj +real dens,dev,cst,sig0 +real mfp0,pt0,rad0,rad10,rad20,rho0,pres0 +real term1,term2,kernel,norm +real, dimension(ngrid, nlayer) :: r0, rho, rn, ntot, mtot, mtotnew +real, dimension(ngrid, nlayer, nq) :: numb_new, numb_ini, numb_ini_dis, mass_ini, mass_new +real, dimension(ngrid, nlayer, nres_coag) :: ndis, ndis_new +real, dimension(ngrid,nlayer,nres_coag,nq) :: ndis_tab,rat_tab,ndis_tab_new +real, dimension(ngrid, nlayer) :: zt +real, dimension(ngrid, nlayer, nq) :: zq + +logical, save :: firstcall = .true. +!$OMP THREADPRIVATE(firstcall) +!for boxinterp function: +integer :: t1,t2,p1,p2,m1,m2 + +!====================================================================== +! Main routine + +!----------------------------------------------------------------------- +!*** 1) Initialisations +!----------------------------------------------------------------------- + +!!! Particles properties +dens = rho_dust ! density dust +sig0 = dev_dt ! Standard deviation of the dust distribution +cst = .75 / (pi*dens) * exp( -4.5*sig0**2. ) ! fixed distribution + +!!! Radius and Volumes defined in coagul_init +! rads_coag, vols_coag, deltar_coag, vrat_coag + +!!! Getting tracers and temperature fields +zq(:,:,:)=pq(:,:,:)+pdqfi(:,:,:)*ptimestep +zt(:,:)=pt(:,:)+pdtfi(:,:)*ptimestep + +!!! Atmospheric state +DO l=1,nlayer + rho(:,l)=pplay(:,l)/(r*pt(:,l)) +ENDDO + +!!! Initializing tracers and other fields +pdqcoag(:,:,:)=0.d0 +ndis_tab(:,:,:,:)=0.d0 +ndis_new(:,:,:)=0.d0 +mtot(:,:)=0.d0 +rat_tab(:,:,:,:) = 0.d0 +numb_ini_dis(:,:,:) = 0.d0 +numb_ini(:,:,:) = 0.d0 +numb_new(:,:,:) = 0.d0 +mass_ini(:,:,:) = 0.d0 +mass_new(:,:,:) = 0.d0 + +!!! Preparing the initial log normal distribution +do iq=1,nq + if ((iq.eq.igcm_dust_mass).or.(iq.eq.igcm_stormdust_mass).or.(iq.eq.igcm_topdust_mass)) then + mass_ini(:,:,iq)=max(zq(:,:,iq),1e-15) + numb_ini(:,:,iq)=max(zq(:,:,iq-1),1e-15) + + !!! Get total mass to ensure mass conservation + mtot(:,:)=mtot(:,:)+mass_ini(:,:,iq) + + !!! Initial Lognormal distribution + ! sig0 remains constant + r0(:,:)=(3/4.*mass_ini(:,:,iq)/(pi*numb_ini(:,:,iq)*dens))**(1/3.)*exp(-1.5*sig0**2) + + do i=1,nres_coag ! Sum of the ndis for all dust modes + ndis_tab(:,:,i,iq)=numb_ini(:,:,iq)*rho(:,:)*deltar_coag(i)/(2.*rads_coag(i)*(2*pi)**0.5*sig0)*exp(-0.5*(log(rads_coag(i)/r0(:,:)))**2/sig0**2) + enddo + !!! Initial number mixing ratio from bin discretization nb/kg + numb_ini_dis(:,:,iq)=sum(ndis_tab(:,:,:,iq),3)/rho(:,:) + numb_new(:,:,iq)=numb_ini_dis(:,:,iq) + endif +enddo + + +!!! Defining initial number ratio for each mode +do iq=1,nq + if ((iq.eq.igcm_dust_mass).or.(iq.eq.igcm_stormdust_mass).or.(iq.eq.igcm_topdust_mass)) then + rat_tab(:,:,:,iq)=ndis_tab(:,:,:,iq)/sum(ndis_tab(:,:,:,:),4) + endif +enddo + +!!! Normalizing the ratio (make sure sum=1) +do iq=1,nq + if ((iq.eq.igcm_dust_mass).or.(iq.eq.igcm_stormdust_mass).or.(iq.eq.igcm_topdust_mass)) then + rat_tab(:,:,:,iq)=rat_tab(:,:,:,iq)/sum(rat_tab(:,:,:,:),4) + endif +enddo + +!!! Total initial distribution nb/m3 per bin +ndis(:,:,:)=sum(ndis_tab(:,:,:,:),4) +!!! Total initial distribution nb/kg +ntot(:,:)=sum(ndis(:,:,:),3)/rho(:,:) + +!----------------------------------------------------------------------- +!*** 2) Coagulation +!----------------------------------------------------------------------- +! Spatial loop => 1D + +if (full_coag_equations) then ! computing full coagulation equations + do i=1,ngrid + do j=1,nlayer + if (ntot(i,j).gt.1000.) then ! activating coagulation above this threshold + pt0=pt(i,j) + rho0=rho(i,j) + mfp0=mfp(pt0,rho0) + + ! Bin loop + do k=1,nres_coag + + ! Term 1 + term1=0. + do jj=1,k + do ii=1,k-1 + kernel=0. + if (dust_coag_kernel_b) kernel=kernel+betab(pt0,mfp0,rads_coag(ii),rads_coag(jj)) + if (dust_coag_kernel_g) kernel=kernel+betag(pt0,mfp0,rads_coag(ii),rads_coag(jj)) + if (dust_coag_kernel_de) kernel=kernel+betade(pt0,mfp0,rho0,rads_coag(ii),rads_coag(jj)) + if (dust_coag_kernel_ti) kernel=kernel+betati(pt0,mfp0,rho0,rads_coag(ii),rads_coag(jj)) + term1=term1+frac(ii,jj,k)*vols_coag(ii)*kernel*ndis_new(i,j,ii)*ndis(i,j,jj) + enddo + enddo + term1=term1*ptimestep/vols_coag(k) + + ! Term 2 + term2=0. + do jj=1,nres_coag + kernel=0. + if (dust_coag_kernel_b) kernel=kernel+betab(pt0,mfp0,rads_coag(k),rads_coag(jj)) + if (dust_coag_kernel_g) kernel=kernel+betag(pt0,mfp0,rads_coag(k),rads_coag(jj)) + if (dust_coag_kernel_de) kernel=kernel+betade(pt0,mfp0,rho0,rads_coag(k),rads_coag(jj)) + if (dust_coag_kernel_ti) kernel=kernel+betati(pt0,mfp0,rho0,rads_coag(k),rads_coag(jj)) + term2=term2+(1-frac(k,jj,k))*kernel*ndis(i,j,jj) + enddo + term2=term2*ptimestep + + ndis_new(i,j,k)=(ndis(i,j,k)+coal_fac*term1)/(1.+coal_fac*term2) + + enddo + + norm=sum(ndis(i,j,:)*vols_coag(:))/sum(ndis_new(i,j,:)*vols_coag(:)) + ndis_new(i,j,:)=ndis_new(i,j,:)*norm + + endif + enddo + enddo + +else ! full_coag_equations + + ! Using Lookup tables + do i=1,ngrid + do j=1,nlayer + if (ntot(i,j).gt.1000.) then ! activating coagulation above this threshold + pt0=pt(i,j) + rho0=rho(i,j) + mfp0=mfp(pt0,rho0) + pres0=pplay(i,j) + !find the corners for interpolation + t2=findval(table_pt,pt0) + if (t2==0) t2=table_numt + if (t2==1) t2=2 + t1=t2-1 + m2=findval(table_mfp,mfp0) + if (m2==0) m2=table_numm + if (m2==1) m2=2 + m1=m2-1 + p2=findval(table_pres,pres0) + if (p2==0) p2=table_numt + if (p2==1) p2=2 + p1=p2-1 + + ! Bin loop + do k=1,nres_coag + + ! Term 1 + term1=0. + do jj=1,k + do ii=1,k-1 + kernel=0. + if (dust_coag_kernel_b) kernel=kernel+boxinterp(table_b(t1,m1,ii,jj), & + table_b(t1,m2,ii,jj), & + table_b(t2,m1,ii,jj), & + table_b(t2,m2,ii,jj), & + table_pt(t1),table_pt(t2),log10(table_mfp(m1)),log10(table_mfp(m2)),pt0,log10(mfp0)) + if (dust_coag_kernel_g) kernel=kernel+boxinterp(table_g(t1,m1,ii,jj), & + table_g(t1,m2,ii,jj), & + table_g(t2,m1,ii,jj), & + table_g(t2,m2,ii,jj), & + table_pt(t1),table_pt(t2),log10(table_mfp(m1)),log10(table_mfp(m2)),pt0,log10(mfp0)) + if (dust_coag_kernel_de) kernel=kernel+boxinterp(table_de(t1,p1,ii,jj), & + table_de(t1,p2,ii,jj), & + table_de(t2,p1,ii,jj), & + table_de(t2,p2,ii,jj), & + table_pt(t1),table_pt(t2),log10(table_pres(p1)),log10(table_pres(p2)),pt0,log10(pres0)) + ! if (dust_coag_kernel_ti) kernel=kernel+betati(pt0,mfp0,rho(i,j,l),ii,jj) !not implemented yet + term1=term1+frac(ii,jj,k)*vols_coag(ii)*kernel*ndis_new(i,j,ii)*ndis(i,j,jj) + enddo + enddo + term1=term1*ptimestep/vols_coag(k) + + + ! Term 2 + term2=0. + do jj=1,nres_coag + kernel=0. + if (dust_coag_kernel_b) kernel=kernel+boxinterp(table_b(t1,m1,k,jj), & + table_b(t1,m2,k,jj), & + table_b(t2,m1,k,jj), & + table_b(t2,m2,k,jj), & + table_pt(t1),table_pt(t2),log10(table_mfp(m1)),log10(table_mfp(m2)),pt0,log10(mfp0)) + if (dust_coag_kernel_g) kernel=kernel+boxinterp(table_g(t1,m1,k,jj), & + table_g(t1,m2,k,jj), & + table_g(t2,m1,k,jj), & + table_g(t2,m2,k,jj), & + table_pt(t1),table_pt(t2),log10(table_mfp(m1)),log10(table_mfp(m2)),pt0,log10(mfp0)) + if (dust_coag_kernel_de) kernel=kernel+boxinterp(table_de(t1,p1,k,jj), & + table_de(t1,p2,k,jj), & + table_de(t2,p1,k,jj), & + table_de(t2,p2,k,jj), & + table_pt(t1),table_pt(t2),log10(table_pres(p1)),log10(table_pres(p2)),pt0,log10(pres0)) + ! if (dust_coag_kernel_ti) kernel=kernel+betati(pt0,mfp0,rho(i,j,l),k,jj) !not implemented yet + term2=term2+(1-frac(k,jj,k))*kernel*ndis(i,j,jj) + enddo + term2=term2*ptimestep + + ndis_new(i,j,k)=(ndis(i,j,k)+coal_fac*term1)/(1.+coal_fac*term2) + enddo + + !! Volume conservation + norm=sum(ndis(i,j,:)*vols_coag(:))/sum(ndis_new(i,j,:)*vols_coag(:)) + ndis_new(i,j,:)=ndis_new(i,j,:)*norm + + endif + + enddo + enddo +endif ! full_coag_equations + +!----------------------------------------------------------------------- +!*** 3) new moments mass and number +!----------------------------------------------------------------------- +mtotnew(:,:)=1.d-15 + +do iq=1,nq + if ((iq.eq.igcm_dust_mass).or.(iq.eq.igcm_stormdust_mass).or.(iq.eq.igcm_topdust_mass)) then + !! Estimate new distribution in each dust mode using the initial number ratios + ndis_tab_new(:,:,:,iq)=ndis_new(:,:,:)*rat_tab(:,:,:,iq) + + do i=1,ngrid + do j=1,nlayer + !! New number and mass mixing ratio in each mode + numb_new(i,j,iq)=sum(ndis_tab_new(i,j,:,iq))/rho(i,j) + mass_new(i,j,iq)=sum(ndis_tab_new(i,j,:,iq)*vols_coag(:))*dens/rho(i,j) + enddo + enddo + + !! New total mass + mtotnew(:,:)=mtotnew(:,:)+mass_new(:,:,iq) + endif +enddo + +do iq=1,nq + if ((iq.eq.igcm_dust_mass).or.(iq.eq.igcm_stormdust_mass).or.(iq.eq.igcm_topdust_mass)) then + !! mass conservation + mass_new(:,:,iq)=mass_new(:,:,iq)/mtotnew(:,:)*mtot(:,:) + !! New tendancy + pdqcoag(:,:,iq)=(mass_new(:,:,iq)-mass_ini(:,:,iq))/ptimestep + !pdqcoag(:,:,iq-1)=(numb_new(:,:,iq)-numb_ini(:,:,iq))/ptimestep + pdqcoag(:,:,iq-1)=(numb_new(:,:,iq)-numb_ini_dis(:,:,iq))/ptimestep + !! check + where (zq(:,:,iq)+pdqcoag(:,:,iq)*ptimestep.le.0.) + pdqcoag(:,:,iq)=-zq(:,:,iq)/ptimestep+1.e-15 + end where + where (zq(:,:,iq-1)+pdqcoag(:,:,iq-1)*ptimestep.le.0.) + pdqcoag(:,:,iq-1)=-zq(:,:,iq-1)/ptimestep+1.e-14 + end where + endif +enddo + +end subroutine dust_coagulation_main + + +!====================================================================== +!====================================================================== +real function dynvis(pt) + implicit none + ! dynamic viscosity following sutherland's formula for CO2 + ! Pa.s = kg m-1 s-1 + real, intent(in) :: pt + dynvis=1.37e-5*(273.15 + 222.)/(pt + 222.)*(pt/273.15)**(3./2.) +end function dynvis + +real function kinvis(pt,rho) + implicit none + real, intent(in) :: pt,rho + kinvis=dynvis(pt)/rho +end function kinvis + +real function thervel(pt,massm) + implicit none + ! thermal velocity of an air molecule or a dust particle + ! m s-1 + real, intent(in) :: pt, massm + thervel=(8*kbz*pt/(pi*massm))**0.5 +end function thervel + +real function mfp(pt,rho) + implicit none + ! thermal velocity of an air molecule + ! m s-1 + real, intent(in) :: pt, rho + mfp=2.*kinvis(pt,rho)/thervel(pt,m0co2) ! mean free path (m) +end function mfp + +real function cun(pt,mfp,rad) + implicit none + !! Cunningham slip flow correction + real, intent(in) :: pt, mfp, rad + cun=1.+mfp/rad*(1.246+0.42*exp(-0.87/(mfp/rad))) +end function cun + +real function diff(pt,mfp,rad) + implicit none + !! diffusion coefficient + real, intent(in) :: pt, mfp, rad + diff=kbz*pt/(6.*pi*rad*dynvis(pt))*cun(pt,mfp,rad) +end function diff + +real function lambda(pt, mfp, rad) + implicit none + !! particle mean free path + real, intent(in) :: pt,mfp,rad + real massm,dens + dens = rho_dust ! density dust + massm=dens*4./3.*pi*rad**3 + lambda=2.*diff(pt, mfp, rad)/(pi*thervel(pt,massm)) +end function lambda + +real function delta(pt,mfp,rad) + implicit none + !! mean distance + real, intent(in) :: pt,mfp,rad + real :: ltmp + + ltmp=lambda(pt,mfp,rad) + delta=( (2.*rad+ltmp)**3 - (4.*rad**2+ltmp**2)**1.5 ) /(6*rad*ltmp) - 2*rad +end function delta + +real function fallv(pt,mfp,rad) + implicit none + !! fall velocity + real, intent(in) :: pt,mfp,rad + real dens + dens = rho_dust ! density dust + fallv=2./9.*rad**2*dens*g/dynvis(pt)*cun(pt,mfp,rad) +end function fallv + +real function reyn(pt,mfp,rad,rho) + implicit none + !! Reynold number + real, intent(in) :: pt,rad,mfp,rho + reyn=2.*rad*fallv(pt,mfp,rad)/kinvis(pt,rho) +end function reyn + +real function stokes(pt,mfp,rad1,rad2) + implicit none + !! Stokes number + real, intent(in) :: pt,mfp,rad1,rad2 + real :: f1 + + f1=fallv(pt,mfp,rad1) + stokes=f1*abs(fallv(pt,mfp,rad2)-f1)/(rad2*g) +end function stokes + +real function schmidt(pt,mfp,rho,rad) + implicit none + !! Schmidt number + real, intent(in) :: pt,mfp,rho,rad + schmidt=kinvis(pt,rho)/diff(pt,mfp,rad) +end function schmidt + +real function coal(pt,mfp,rad1,rad2,rho) + implicit none + !! Coalescence efficiency + real, intent(in) :: pt,mfp,rad1,rad2,rho + real :: coalea,coalev + real :: stmp,rtmp + + stmp=stokes(pt,mfp,rad1,rad2) + rtmp=reyn(pt,mfp,rad2,rho) + coalea=stmp**2/(stmp+0.5)**2 + if (stmp>1.214) then + coalev=(1+0.75*log(2.*stmp)/(stmp-1.214))**(-2) + else + coalev=0. + endif + coal=(60.*coalev+coalea*rtmp)/(60.+rtmp) +end function coal + +real function betab(pt,mfp,rad1,rad2) + implicit none + real, intent(in) :: pt,mfp,rad1,rad2 + real :: m1,m2,num,den1,den2 + real :: dens + real :: d1,d2 + + dens = rho_dust ! density dust + d1=diff(pt,mfp,rad1) + d2=diff(pt,mfp,rad2) + !! kernel transition regime + ! Mass of 1 particle 1 and 2 + m1=dens*4./3.*pi*rad1**3 + m2=dens*4./3.*pi*rad2**3 + num=4.*pi*(rad1+rad2)*(d1+d2) + den1=(rad1+rad2)/(rad1+rad2+(delta(pt,mfp,rad1)**2+delta(pt,mfp,rad2)**2)**0.5) + den2=4./effb*(d1+d2)/((thervel(pt,m1)**2+thervel(pt,m2)**2)**0.5*(rad1+rad2)) + betab=num/(den1+den2) +end function betab + +real function betag(pt,mfp,rad1,rad2) + implicit none + real, intent(in) :: pt,mfp,rad1,rad2 + real w1,w2,ee + real :: radt + !! kernel gravitation + ! fall vel + w1=fallv(pt,mfp,rad1) + w2=fallv(pt,mfp,rad2) + radt=(rad1+rad2)**2 + if (coal_coeff_mode.eq.0) then ! E=1 + ee=1. + elseif (coal_coeff_mode.eq.1) then + ee=1.5*(min(rad1,rad2))**2/radt + elseif (coal_coeff_mode.eq.2) then + ee=0.25*(min(rad1,rad2))**2/radt + endif + !b=coal(pt,r1,r2)*pi*(r1+r2)**2*abs(w1-w2) + betag=ee*pi*radt*abs(w1-w2) +end function betag + +real function betade(pt,mfp,rho,rad1,rad2) + implicit none + real, intent(in) :: pt,mfp,rad1,rad2,rho + real rd,sc + !! kernel diffusion enhancement + rd=reyn(pt,mfp,max(rad1,rad2),rho) + sc=schmidt(pt,mfp,rho,min(rad1,rad2)) + if (rd.le.1.) then + betade=betab(pt,mfp,rad1,rad2)*0.45*rd**(1/3.)*sc**(1/3.) + else + betade=betab(pt,mfp,rad1,rad2)*0.45*rd**(1/2.)*sc**(1/3.) + endif +end function betade + +real function betade_tabfunc(pt,pres,rad1,rad2) + implicit none + real, intent(in) :: pt,pres,rad1,rad2 + real rd,sc,rho,mf + !! kernel diffusion enhancement + rho=pres/(pt*r) + mf=mfp(pt,rho) + rd=reyn(pt,mf,max(rad1,rad2),rho) + sc=schmidt(pt,mf,rho,min(rad1,rad2)) + if (rd.le.1.) then + betade_tabfunc=betab(pt,mf,rad1,rad2)*0.45*rd**(1/3.)*sc**(1/3.) + else + betade_tabfunc=betab(pt,mf,rad1,rad2)*0.45*rd**(1/2.)*sc**(1/3.) + endif +end function betade_tabfunc + +real function betati(pt,mfp,rho,rad1,rad2) + implicit none + real, intent(in) :: pt,mfp,rad1,rad2,rho + real w1,w2,eps + ! fall vel + w1=fallv(pt,mfp,rad1) + w2=fallv(pt,mfp,rad2) + betati=eps**(0.75)*pi/(g*kinvis(pt,rho)**(0.25))*(rad1+rad2)**2*abs(w1-w2) +end function betati + +real function betati_tabfunc(pt,pres,rad1,rad2) + implicit none + real, intent(in) :: pt,rad1,rad2 + real w1,w2,eps,mf,rho,pres + rho=pres/(pt*r) + mf=mfp(pt,rho) + ! fall vel + w1=fallv(pt,mf,rad1) + w2=fallv(pt,mf,rad2) + betati_tabfunc=eps**(0.75)*pi/(g*kinvis(pt,rho)**(0.25))*(rad1+rad2)**2*abs(w1-w2) +end function betati_tabfunc + +real function frac(i,j,k) + implicit none + integer, intent(in) :: i,j,k + real vint + ! intermediate volume + vint=vols_coag(i)+vols_coag(j) + if (k.lt.nres_coag) then + if ( (vint.ge.vols_coag(k)) .and. (vint.lt.vols_coag(k+1))) then + frac=(vols_coag(k+1)-vint)/(vols_coag(k+1)-vols_coag(k))*vols_coag(k)/vint + return + endif + endif + if ((k.eq.nres_coag) .and. (vint.ge.vols_coag(k))) then + frac=1. + return + endif + if (k.gt.1) then + if ( (vint.lt.vols_coag(k)) .and. (vint.gt.vols_coag(k-1))) then + frac=(-vols_coag(k-1)+vint)/(-vols_coag(k-1)+vols_coag(k))*vols_coag(k)/vint + return + endif + endif + + frac=0. + return +end function frac + +! ********************************************************* +! Building lookup tables : called in the firstcall +! ********************************************************* + +subroutine make_tables() + +implicit none +integer :: i,j,k,l +real :: t1,t2,p1,p2 + +do i = 1,table_numt + table_pt(i)=50.+25.*real(i) +end do + +do j = 1,table_nump + table_pres(j)=(1.e-8)*10.**(real(j)/2.) +end do + +do j = 1,table_numm + table_mfp(j)=(1.e-8)*10.**(real(j)/2.) +end do + +do i = 1,table_numt + do j = 1,table_numm + do k = 1,nres_coag + do l = 1,nres_coag + table_b(i,j,k,l)=betab(table_pt(i),table_mfp(j),rads_coag(k),rads_coag(l)) + table_g(i,j,k,l)=betag(table_pt(i),table_mfp(j),rads_coag(k),rads_coag(l)) + end do + end do + end do +end do + +do i = 1,table_numt + do j = 1,table_nump + do k = 1,nres_coag + do l = 1,nres_coag + table_de(i,j,k,l)=betade_tabfunc(table_pt(i),table_pres(j),rads_coag(k),rads_coag(l)) + end do + end do + end do +end do + +end subroutine make_tables + +!======================================================================C +!======================================================================C + +real function boxinterp(ft1p1,ft2p1,ft1p2,ft2p2,tref1,tref2,pref1, & + pref2,tmid,pmid) +! +! Calculate 2d interpolation from pressure and temperature +!----------------------------------------------------------------------! +! T2 .FT2P1 .FT2P2 +! +! +! T1 .FT1P1 .FT1P2 +! P1 P2 +!----------------------------------------------------------------------! + +implicit none + +real :: ft1p1, ft2p1, ft1p2, ft2p2, tref1, tref2 +real :: pref1, pref2, tmid, pmid, ans1, ans2, ans + +!======================================================================C + +ans1 = ft1p1 + (ft2p1 - ft1p1)*(tmid - tref1)/(tref2 - tref1) +ans2 = ft1p2 + (ft2p2 - ft1p2)*(tmid - tref1)/(tref2 - tref1) +boxinterp = ans1 + (ans2 - ans1)*(pmid - pref1)/(pref2 - pref1) + +end function boxinterp + +!======================================================================C +!======================================================================C + +integer function findval(array,value) +implicit none +real :: array(:) +real :: value + +findval=minloc(array,dim=1,mask = (array > value)) + +end function findval + +end module dust_coagulation_mod + Index: /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F (revision 4062) +++ /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F (revision 4063) @@ -25,5 +25,6 @@ use calcstormfract_mod, only: calcstormfract use topmons_mod, only: topmons,topmons_setup - use coagulation_mod, only: coagul_main,coagul_init + use dust_coagulation_mod, only: dust_coagulation_main, + & dust_coagulation_init use nltecool_mod, only: nltecool use nlte_tcool_mod, only: nlte_tcool @@ -150,5 +151,5 @@ use callkeys_mod, only: photochem, callthermos use callkeys_mod, only: startphy_file - use callkeys_mod, only: coagulation + use callkeys_mod, only: dust_coagulation IMPLICIT NONE @@ -794,7 +795,7 @@ if (topflows) call topmons_setup(ngrid) -c Initialize coagulation parameters +c Initialize dust coagulation parameters c ~~~~~~~~~~~~~~~ - if (coagulation) call coagul_init() + if (dust_coagulation) call dust_coagulation_init() c Parameterization of the ATKE routine @@ -1455,7 +1456,7 @@ c 3.3 Dust coagulation c ------------------------------------------- - IF (coagulation) THEN + IF (dust_coagulation) THEN pdqcoag(:,:,:)=0. - CALL coagul_main(ngrid,nlayer,nq,ptime,ptimestep, + CALL dust_coagulation_main(ngrid,nlayer,nq,ptime,ptimestep, & pq,pdq,pt,pdt,pplay,pplev, & pdqcoag) @@ -1496,5 +1497,5 @@ ENDIF ! end of if (rdstorm) - ENDIF ! end of if (coagulation) + ENDIF ! end of if (dust coagulation) c 3.4 Dust injection from the surface @@ -3791,10 +3792,10 @@ endif ! (topflows) - if (coagulation) then + if (dust_coagulation) then call write_output('zdqcoag_dustm', - & 'coagulation tendency', + & 'dust coagulation tendency', & 'kg kg-1 s-1',pdqcoag(:,:,igcm_dust_mass)) call write_output('zdqcoag_dustn', - & 'coagulation tendency', + & 'dust coagulation tendency', & 'nbp kg-1 s-1',pdqcoag(:,:,igcm_dust_number)) endif