Skip to content

Commit

Permalink
ini push
Browse files Browse the repository at this point in the history
  • Loading branch information
Jingyi Chen committed Feb 18, 2022
0 parents commit 9ee7d41
Show file tree
Hide file tree
Showing 27 changed files with 22,120 additions and 0 deletions.
18,952 changes: 18,952 additions & 0 deletions DVODE_F90_M.f90

Large diffs are not rendered by default.

62 changes: 62 additions & 0 deletions aerospec.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,62 @@
!****************************************************************
! This module provides dry aerosol spectra
!****************************************************************

MODULE aerospecMD

use constant

implicit none

integer :: naero = 50

CONTAINS

SUBROUTINE aerospec(iaero,aerobin,na0,rhos,solfracv,rhois,aerostate0,&
&conc0,prad0,bprad0,pmass0)

implicit none
integer, intent(in) :: aerobin,iaero
real(RLK), intent(in) :: rhos,solfracv,rhois
real(RLK), intent(out) :: na0
integer, dimension(aerobin),intent(out) :: aerostate0
real(RLK), dimension(aerobin),intent(out) :: conc0,prad0,pmass0
real(RLK), dimension(aerobin+1),intent(out) :: bprad0

real(RLK) :: dev,geo,radmin,radmax,raddiff,coe1,coe2
integer :: ii

dev=1.5
geo=0.06e-4

if (iaero .eq. 1) then
na0 = 200
endif

radmin=0.01e-4;radmax=1.0e-4
raddiff=(log(radmax)-log(radmin))/real(aerobin,8)
bprad0(1)=radmin

do ii=1,aerobin
bprad0(ii+1)=exp(log(bprad0(ii))+raddiff)
prad0(ii)=(bprad0(ii)+bprad0(ii+1))/2.0

coe1=na0/(2.0*pi)**0.5/log(dev)
coe2=2.0*log(dev)**2.0
conc0(ii)=(coe1/prad0(ii))*exp(-(log(prad0(ii))-log(geo))**2.0/coe2)!#/CM3/CM
conc0(ii)=conc0(ii)*(bprad0(ii+1)-bprad0(ii)) ! #/CM3

pmass0(ii)=rhos*solfracv*4.0/3.0*pi*prad0(ii)**3.0
pmass0(ii)=pmass0(ii)+rhois*(1.0-solfracv)*4.0/3.0*pi*prad0(ii)**3.0

enddo

!aerosol states
do ii=1,aerobin
aerostate0(ii)=0
enddo


END SUBROUTINE aerospec

END MODULE aerospecMD
Binary file added cdconc.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
113 changes: 113 additions & 0 deletions chem.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,113 @@
!****************************************************************
! This module provides aerosol chemical compositions
!****************************************************************

MODULE chemMD

use constant

implicit none


CONTAINS

SUBROUTINE chem(ichem,rhos,rhosap,soluba,solubb,solubc,Ms,kap,rhois,&
& solfracm,solfracv)

implicit none

integer, intent(in) :: ichem
real(RLK), intent(out) :: rhos,rhosap,soluba,solubb,solubc,&
& Ms,kap,rhois,solfracm,solfracv
integer :: solid,inid

! soid : soluble part chemical composition:
! 1-sulfate ammonium; 2-sea salt; 3-organic carbon
! inid : insoluble part chemical composition:
! 1-dust;
! sfr : soluble fraction (mass fraction)
solid = ichem
inid = 1
solfracm = 1.0

call solcomp(solid,rhos,rhosap,soluba,solubb,solubc,Ms,kap)
call insolcomp(inid,rhois)

solfracv=solfracm*rhois/(rhos+solfracm*(rhois-rhos))

END SUBROUTINE chem


!*************************************************************
!This subroutine provides parameters of soluble coomposition
!**************************************************************
SUBROUTINE solcomp(idin,rhost,rhosapt,solubat,solubbt,solubct,Mst,kapt)

implicit none

integer, intent(in) :: idin
real(RLK), intent(out) :: rhost,rhosapt,solubat,solubbt,solubct,Mst,kapt

select case (idin)

!SULFATE AMMONIUM
case (1)
rhost=1.769
rhosapt=2.093
solubat=0.1149
solubbt=-4.489e-4
solubct=1.385e-6
Mst=132.14
kapt=0.61

!Sodium Chloride
case (2)
rhost=2.165
rhosapt=2.726
solubat=0.1805
solubbt=-5.310e-4
solubct=9.965e-7
Mst=58.44
kapt=1.28

!biomass burning aerosols
case (3)
rhost=1.662
rhosapt=2.000
solubat=0.1149
solubbt=-4.489e-4
solubct=1.385e-6
Mst=111.00
kapt=0.1
case default
print *,'THIS SALT DOES NOT EXIST IN THIS MODEL.'
end select

END SUBROUTINE solcomp

!***************************************************************
!This subroutine provides parameters of insoluble composition.
!**************************************************************
SUBROUTINE insolcomp(id,rhois1)

implicit none
integer, intent(in) :: id
real(RLK), intent(out) :: rhois1

select case (id)
case (1) !dust
rhois1=2.6


case (2) !soot BC
rhois1=2.0


case default
print *,'INPUT INSOLUBLE COMPOSITION DOES NOT EXIT.'

end select

END SUBROUTINE insolcomp

END MODULE chemMD
204 changes: 204 additions & 0 deletions cloudspec.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,204 @@
!*********************************************************************************
! This Module provides cloud droplet spectra statistical measures
!********************************************************************************

MODULE cloudspecMD

use constant

implicit none

Contains

SUBROUTINE cloudspec(i0bin,aerobin,aerostate,concv,&
&prad,pmass,cdconc,cdbeta,cdrmean,cdrsd,cddisp,acbin,&
&lwc,re,rv,skew,kurt,rrainLN,ThrsFuncLN)

implicit none

integer,parameter :: neye=3
integer, intent(in) :: i0bin,aerobin
integer, dimension(aerobin), intent(in) :: aerostate
real(RLK), dimension(aerobin), intent(in) :: prad,pmass,concv
integer, dimension(neye),intent(out) :: acbin
real(RLK),dimension(neye),intent(out) :: cdconc,cdbeta,cdrmean,cdrsd,cddisp,lwc
real(RLK),dimension(neye),intent(out) :: skew,kurt,rrainLN,ThrsFuncLN,re,rv

real(RLK),dimension(neye) :: cd_cm2, cd_cm3, cd_cm4, cd_m2, cd_m3, cd_m6
real(RLK), dimension(neye) :: rainLN_m6, rainLN_m3, vc
integer, dimension(neye) :: eye
integer :: ii, jj

!eye: 1-Reutter et al. 2009; 2-radius>1um; 3-radius>Kohler critical radius
cdconc=0.0 !number concentration per volume
cdrmean=0.0 !mean radius
cdrsd=0.0 !standard deviation
cddisp=0.0 !relative dispersion
cdbeta=0.0

!central moment of cloud droplet spectrum
cd_cm2=0.0;cd_cm3=0.0;cd_cm4=0.0

!moment of cloud droplet spectrum
cd_m2=0.0; cd_m3=0.0; cd_m6=0.0

acbin=0 !activated particle bin
lwc=0.0 !liqudi water content
vc=0.0 !total cloud water volume

!ThrsFuncLR6=0
ThrsFuncLN=0.0

eye=0
do ii=i0bin,aerobin

CALL valueeye(eye,neye,aerostate(ii),prad(ii))

do jj=1,neye
if (eye(jj)==1) then

acbin(jj)=acbin(jj)+1
cdconc(jj)=cdconc(jj)+concv(ii)
lwc(jj)=lwc(jj)+concv(ii)*pmass(ii)
vc(jj)=vc(jj)+concv(ii)*4.0/3.0*pi*prad(ii)**3.0 !cm3/cm3

cdrmean(jj)=cdrmean(jj)+prad(ii)*concv(ii)
cd_m2(jj)=cd_m2(jj)+(prad(ii)**2.0)*concv(ii)
cd_m3(jj)=cd_m3(jj)+(prad(ii)**3.0)*concv(ii)
cd_m6(jj)=cd_m6(jj)+(prad(ii)**6.0)*concv(ii)
else
acbin(jj)=0
cdconc(jj)=0.0
lwc(jj)=0.0
vc(jj)=0.0

cdrmean(jj)=0.0
cd_m2(jj)=0.0
cd_m3(jj)=0.0
cd_m6(jj)=0.0
endif
enddo

enddo

do jj=1,neye
if(cdconc(jj)/=0.0) then
cdrmean(jj)=cdrmean(jj)/cdconc(jj)
rv(jj)=(vc(jj)/cdconc(jj)/4.0*3.0/pi)**(1.0/3.0) ! volume mean radius
re(jj)=cd_m3(jj)/cd_m2(jj)
else
cdrmean(jj)=0.0
rv(jj)=0.0
re(jj)=0.0
endif
enddo

cdbeta=re/rv

eye=0
do ii=i0bin,aerobin

CALL valueeye(eye,neye,aerostate(ii),prad(ii))

do jj=1,neye
if(eye(jj)==1) then
cd_cm2(jj)=cd_cm2(jj)+((prad(ii)-cdrmean(jj))**2.0)*concv(ii)
cd_cm3(jj)=cd_cm3(jj)+((prad(ii)-cdrmean(jj))**3.0)*concv(ii)
cd_cm4(jj)=cd_cm4(jj)+((prad(ii)-cdrmean(jj))**4.0)*concv(ii)
else
cd_cm2(jj)=0.0
cd_cm3(jj)=0.0
cd_cm4(jj)=0.0
endif
enddo

enddo

do jj=1,neye
if(cdconc(jj)/=0.0) then
cdrsd(jj)=sqrt(cd_cm2(jj)/cdconc(jj))
cddisp(jj)=cdrsd(jj)/cdrmean(jj)
skew(jj)=cd_cm3(jj)/cdconc(jj)/(cdrsd(jj)**3.0)
kurt(jj)=cd_cm4(jj)/cdconc(jj)/(cdrsd(jj)**4.0)
else
cdrsd(jj)=0.0
cddisp(jj)=0.0
skew(jj)=0.0
kurt(jj)=0.0
endif
enddo

! Liu et al. (2004) eq.5 and eq.7
do jj=1,neye
!rrainLR6(jj)=23.72*1.0e-7/lwc(jj)**(1.0/6.0)/(rv(jj)*1.12)**(1.0/2.0)
rrainLN(jj)=2.8522e-6*cdconc(jj)**(1.0/6.0)/lwc(jj)**(1.0/3.0)
enddo

!rainLR6_m6=0; rainLR6_m3=0
rainLN_m6=0.0; rainLN_m3=0.0

do ii=i0bin,aerobin

do jj=1,neye

! if (prad(ii)>=rrainLR6(jj)) then
! rainLR6_m6(jj)=rainLR6_m6(jj)+(prad(ii)**6.0)*concv(ii)
! rainLR6_m3(jj)=rainLR6_m3(jj)+(prad(ii)**3.0)*concv(ii)
! else
! rainLR6_m6(jj)=0
! rainLR6_m3(jj)=0
! endif

if (prad(ii)>=rrainLN(jj)) then
rainLN_m6(jj)=rainLN_m6(jj)+(prad(ii)**6.0)*concv(ii)
rainLN_m3(jj)=rainLN_m3(jj)+(prad(ii)**3.0)*concv(ii)
else
rainLN_m6(jj)=0.0
rainLN_m3(jj)=0.0
endif

enddo

enddo


do jj=1,neye

! ThrsFuncLR6(jj)=rainLR6_m6(jj)*rainLR6_m3(jj)/cd_m6(jj)/cd_m3(jj)

ThrsFuncLN(jj)=rainLN_m6(jj)*rainLN_m3(jj)/cd_m6(jj)/cd_m3(jj)

enddo

END SUBROUTINE cloudspec

SUBROUTINE valueeye(eye,neye,aerostatei,pradi)

implicit none

integer, intent(in) :: neye
integer, dimension(neye), intent(out) :: eye
real(RLK), intent(in) :: pradi
integer, intent(in) :: aerostatei

if(aerostatei==2.or.aerostatei==4) then
eye(1)=1
else
eye(1)=0
endif

if (pradi>1.0E-04) then
eye(2)=1
else
eye(2)=0
endif

if (aerostatei==3.or.aerostatei==4) then
eye(3)=1
else
eye(3)=0
endif

END SUBROUTINE valueeye

END MODULE cloudspecMD
Loading

0 comments on commit 9ee7d41

Please sign in to comment.