******************************************************* * * subroutine nucleaCO2(pco2,temp,sat,n_ccn,nucrate, & n_ccn_h2oice,rad_h2oice,nucrate_h2oice) USE comcstfi_h implicit none * * * This subroutine computes the nucleation rate * * as given in Pruppacher & Klett (1978) in the * * case of water ice forming on a solid substrate. * * Definition refined by Keese (jgr,1989) * * Authors: F. Montmessin * * Adapted for the LMD/GCM by J.-B. Madeleine * * (October 2011) * * Optimisation by A. Spiga (February 2012) * ******************************************************* ! nucrate = output ! nucrate_h2o en sortie aussi : !nucleation sur dust et h2o separement ici !#include "tracer.h" #include "microphys.h" c#include "microphysCO2.h" c Inputs DOUBLE PRECISION pco2,sat DOUBLE PRECISION n_ccn(nbinco2_cld), n_ccn_h2oice(nbinco2_cld) REAL temp c Output ! DOUBLE PRECISION nucrate(nbinco2_cld) DOUBLE PRECISION nucrate(nbinco2_cld) DOUBLE PRECISION nucrate_h2oice(nbinco2_cld) ! h2o as substrate double precision rad_h2oice(nbinco2_cld) ! h2o ice grid (as substrate) c Local variables DOUBLE PRECISION nco2 c DOUBLE PRECISION sigco2 ! Water-ice/air surface tension (N.m) c external sigco2 DOUBLE PRECISION rstar ! Radius of the critical germ (m) DOUBLE PRECISION gstar ! # of molecules forming a critical embryo DOUBLE PRECISION fistar ! Activation energy required to form a critical embryo (J) ! DOUBLE PRECISION zeldov ! Zeldovitch factor (no dim) DOUBLE PRECISION fshapeco2 ! function defined at the end of the file DOUBLE PRECISION deltaf c Ratio rstar/radius of the nucleating dust particle c double precision xratio double precision mtetalocal ! local mteta in double precision double precision fshapeco2simple,zefshapeco2 integer i LOGICAL firstcall DATA firstcall/.true./ SAVE firstcall c ************************************************* mtetalocal = dble(mtetaco2) !! use mtetalocal for better performance cccccccccccccccccccccccccccccccccccccccccccccccccc ccccccccccc ESSAIS TN MTETA = F (T) cccccccccccccc c if (temp .gt. 200) then c mtetalocal = mtetalocal c else if (temp .lt. 190) then c mtetalocal = mtetalocal-0.05 c else c mtetalocal = mtetalocal - (190-temp)*0.005 c endif c----------------exp law, see Trainer 2008, J. Phys. Chem. C 2009, 113, 2036\u20132040 !mtetalocal = max(mtetalocal - 6005*exp(-0.065*temp),0.1) !mtetalocal = max(mtetalocal - 6005*exp(-0.068*temp),0.1) !print*, mtetalocal, temp cccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccc c IF (firstcall) THEN c print*, ' ' c print*, 'dear user, please keep in mind that' c print*, 'contact parameter IS constant' !print*, 'contact parameter IS NOT constant:' !print*, 'max(mteta - 6005*exp(-0.065*temp),0.1)' !print*, 'max(mteta - 6005*exp(-0.068*temp),0.1)' c print*, ' ' c firstcall=.false. c END IF cccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccc c write(*,*) "IN nuc, SAT = ",sat c write(*,*) "IN nuc, mtetalocal = ",mtetalocal if (sat .gt. 1.) then ! minimum condition to activate nucleation nco2 = pco2 / kbz / temp rstar = 2. * sigco2 * vo1co2 / (kbz*temp*dlog(sat)) gstar = 4. * pi * (rstar * rstar * rstar) / (3.*vo1co2) fshapeco2simple = (2.+mtetalocal)*(1.-mtetalocal)*(1.-mtetalocal) & / 4. c Loop over size bins do 200 i=1,nbinco2_cld c write(*,*) "IN NUCLEA, i, RAD_CLDCO2(i) = ",i, rad_cldco2(i), c & n_ccn(i) if ( n_ccn(i) .lt. 1e-10 ) then c no dust, no need to compute nucleation! nucrate(i)=0. goto 210 endif if (rad_cldco2(i).gt.3000.*rstar) then zefshapeco2 = fshapeco2simple else zefshapeco2 = fshapeco2(mtetalocal,rad_cldco2(i)/rstar) endif fistar = (4./3.*pi) * sigco2 * (rstar * rstar) * & zefshapeco2 deltaf = (2.*desorpco2-surfdifco2-fistar)/ & (kbz*temp) deltaf = min( max(deltaf, -100.d0), 100.d0) if (deltaf.eq.-100.) then nucrate(i) = 0. else nucrate(i)= dble(sqrt ( fistar / & (3.*pi*kbz*temp*(gstar*gstar)) ) & * kbz * temp * rstar & * rstar * 4. * pi & * ( nco2*rad_cldco2(i) ) & * ( nco2*rad_cldco2(i) ) & / ( zefshapeco2 * nusco2 * m0co2 ) & * dexp (deltaf)) endif 210 continue if ( n_ccn_h2oice(i) .lt. 1e-10 ) then c no dust, no need to compute nucleation! nucrate_h2oice(i)=0. goto 200 endif if (rad_h2oice(i).gt.3000.*rstar) then zefshapeco2 = fshapeco2simple else zefshapeco2 = fshapeco2(mtetalocal,rad_h2oice(i)/rstar) ! same m for dust/h2o ice endif fistar = (4./3.*pi) * sigco2 * (rstar * rstar) * & zefshapeco2 deltaf = (2.*desorpco2-surfdifco2-fistar)/ & (kbz*temp) deltaf = min( max(deltaf, -100.d0), 100.d0) if (deltaf.eq.-100.) then nucrate_h2oice(i) = 0. else nucrate_h2oice(i)= dble(sqrt ( fistar / & (3.*pi*kbz*temp*(gstar*gstar)) ) & * kbz * temp * rstar & * rstar * 4. * pi & * ( nco2*rad_h2oice(i) ) & * ( nco2*rad_h2oice(i) ) & / ( zefshapeco2 * nusco2 * m0co2 ) & * dexp (deltaf)) endif 200 continue else do i=1,nbinco2_cld nucrate(i) = 0. nucrate_h2oice(i) = 0. enddo endif return end ********************************************************* double precision function fshapeco2(cost,rap) implicit none * function computing the f(m,x) factor * * related to energy required to form a critical embryo * ********************************************************* double precision cost,rap double precision yeah !! PHI yeah = sqrt( 1. - 2.*cost*rap + rap*rap ) !! FSHAPECO2 = TERM A fshapeco2 = (1.-cost*rap) / yeah fshapeco2 = fshapeco2 * fshapeco2 * fshapeco2 fshapeco2 = 1. + fshapeco2 !! ... + TERM B yeah = (rap-cost)/yeah fshapeco2 = fshapeco2 + & rap*rap*rap*(2.-3.*yeah+yeah*yeah*yeah) !! ... + TERM C fshapeco2 = fshapeco2 + 3. * cost * rap * rap * (yeah-1.) !! FACTOR 1/2 fshapeco2 = 0.5*fshapeco2 return end