subroutine improvedCO2clouds(ngrid,nlay,ptimestep, & pplay,pzlev,pt,pdt, & pq,pdq,pdqcloudco2,pdtcloudco2, & nq,tauscaling, & memdMMccn,memdMMh2o,memdNNccn) ! to use 'getin' USE comcstfi_h USE ioipsl_getincom USE updaterad use tracer_mod !, only: rho_ice_co2, nuiceco2_sed, igcm_co2, ! & rho_ice,igcm_h2o_ice, igcm_ccn_number, ! & igcm_co2_ice, igcm_dust_mass, ! & igcm_dust_number, igcm_ccnco2_mass, ! & igcm_ccnco2_number use conc_mod, only: mmean implicit none c------------------------------------------------------------------ c This routine is used to form CO2 clouds when a parcel of the GCM is c saturated. It includes the ability to have supersaturation, a c computation of the nucleation rates, growthrates and the c scavenging of dust particles by clouds. c It is worth noting that the amount of dust is computed using the c dust optical depth computed in aeropacity.F. That's why c the variable called "tauscaling" is used to convert c pq(dust_mass) and pq(dust_number), which are relative c quantities, to absolute and realistic quantities stored in zq. c This has to be done to convert the inputs into absolute c values, but also to convert the outputs back into relative c values which are then used by the sedimentation and advection c schemes. c CO2 ice particles can nucleate on both dust and on water ice particles c When CO2 ice is deposited onto a water ice particles, the particle is c removed from the water tracers. cWARNING: no sedimentation of the water ice origin is performed c in the microphysical timestep in co2cloud.F. c Authors of the water ice clouds microphysics c J.-B. Madeleine, based on the work by Franck Montmessin c (October 2011) c T. Navarro, debug,correction, new scheme (October-April 2011) c A. Spiga, optimization (February 2012) c Adaptation for CO2 clouds by Joachim Audouard (09/16), based on the work c of Constantino Listowski c------------------------------------------------------------------ !#include "dimensions.h" !#include "dimphys.h" #include "callkeys.h" !#include "tracer.h" !#include "comgeomfi.h" !#include "dimradmars.h" #include "microphys.h" #include "datafile.h" !#include "microphysCO2.h" !#include "conc.h" c------------------------------------------------------------------ c Inputs: INTEGER ngrid,nlay integer nq ! nombre de traceurs REAL ptimestep ! pas de temps physique (s) REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa) REAL pzlev(ngrid,nlay) ! altitude au milieu des couches () REAL pt(ngrid,nlay) ! temperature at the middle of the ! layers (K) REAL pdt(ngrid,nlay) ! tendance temperature des autres ! param. REAL pq(ngrid,nlay,nq) ! traceur (kg/kg) REAL pdq(ngrid,nlay,nq) ! tendance avant condensation ! (kg/kg.s-1) REAL tauscaling(ngrid) ! Convertion factor for qdust and Ndust REAL rice(ngrid,nlay) ! Water Ice mass mean radius (m) ! used for nucleation of CO2 on ice-coated ccns REAL rccnh2o(ngrid,nlay) ! Water Ice mass mean radius (m) c Outputs: REAL pdqcloudco2(ngrid,nlay,nq) ! tendance de la condensation ! CO2 (kg/kg.s-1) ! condensation si igcm_co2_ice REAL pdtcloudco2(ngrid,nlay) ! tendance temperature due ! a la chaleur latente c------------------------------------------------------------------ c Local variables: LOGICAL firstcall DATA firstcall/.true./ SAVE firstcall REAL*8 derf ! Error function !external derf !REAL*8 massflowrateCO2 !external massflowrateCO2 INTEGER ig,l,i,flag_pourri REAL zq(ngrid,nlay,nq) ! local value of tracers REAL zq0(ngrid,nlay,nq) ! local initial value of tracers REAL zt(ngrid,nlay) ! local value of temperature REAL zqsat(ngrid,nlay) ! saturation vapor pressure for CO2 REAL lw !Latent heat of sublimation (J.kg-1) REAL,save :: l0,l1,l2,l3,l4 REAL cste DOUBLE PRECISION dMice ! mass of condensed ice DOUBLE PRECISION sumcheck DOUBLE PRECISION pco2 ! Co2 vapor partial pressure (Pa) DOUBLE PRECISION satu ! Co2 vapor saturation ratio over ice DOUBLE PRECISION Mo,No DOUBLE PRECISION Rn, Rm, dev2,dev3, n_derf, m_derf DOUBLE PRECISION memdMMccn(ngrid,nlay) DOUBLE PRECISION memdMMh2o(ngrid,nlay) DOUBLE PRECISION memdNNccn(ngrid,nlay) ! Radius used by the microphysical scheme (m) DOUBLE PRECISION n_aer(nbinco2_cld) ! number concentration volume-1 of particle/each size bin DOUBLE PRECISION m_aer(nbinco2_cld) ! mass mixing ratio of particle/each size bin DOUBLE PRECISION m_aer2(nbinco2_cld) ! mass mixing ratio of particle/each size bin DOUBLE PRECISION m_aer_h2oice2(nbinco2_cld) ! mass mixing ratio of particle/each size bin DOUBLE PRECISION n_aer_h2oice(nbinco2_cld) ! Same - for CO2 nucleation DOUBLE PRECISION m_aer_h2oice(nbinco2_cld) ! Same - for CO2 nucleation DOUBLE PRECISION rad_h2oice(nbinco2_cld) c REAL*8 sigco2 ! Co2-ice/air surface tension (N.m) c EXTERNAL sigco2 DOUBLE PRECISION dN,dM, dNh2o, dMh2o, dNN,dMM,dNNh2o,dMMh2o DOUBLE PRECISION dMh2o_ice,dMh2o_ccn DOUBLE PRECISION rate(nbinco2_cld) ! nucleation rate DOUBLE PRECISION rateh2o(nbinco2_cld) ! nucleation rate REAL seq DOUBLE PRECISION rho_ice_co2T(ngrid,nlay) DOUBLE PRECISION riceco2(ngrid,nlay) ! CO2Ice mean radius (m) REAL rhocloud(ngrid,nlay) ! Cloud density (kg.m-3) REAL rhocloudco2(ngrid,nlay) ! Cloud density (kg.m-3) REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m) c REAL res ! Resistance growth DOUBLE PRECISION Ic_rice ! Mass transfer rate CO2 ice crystal DOUBLE PRECISION ratioh2o_ccn DOUBLE PRECISION vo2co2 c Parameters of the size discretization c used by the microphysical scheme DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-11 ! Minimum radius (m) DOUBLE PRECISION, PARAMETER :: rmax_cld = 3.e-6 ! Maximum radius (m) DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-12 ! Minimum bounary radius (m) DOUBLE PRECISION, PARAMETER :: rbmax_cld = 5.e-6 ! Maximum boundary radius (m) DOUBLE PRECISION vrat_cld ! Volume ratio DOUBLE PRECISION rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m) SAVE rb_cldco2 DOUBLE PRECISION dr_cld(nbinco2_cld) ! width of each rad_cldco2 bin (m) DOUBLE PRECISION vol_cld(nbinco2_cld) ! particle volume for each bin (m3) DOUBLE PRECISION Proba,Masse_atm,drsurdt,reff,Probah2o REAL sigma_iceco2 ! Variance of the ice and CCN distributions SAVE sigma_iceco2 REAL sigma_ice ! Variance of the ice and CCN distributions SAVE sigma_ice DOUBLE PRECISION Niceco2,Qccnco2,Nccnco2 integer coeffh2o !Variables for the meteoritic flux double precision meteo_ccn(ngrid,nlay,100) !100=nbinco2_cld !!! double precision,save :: meteo(130,100) double precision mtemp(100) logical file_ok integer altitudes_meteo(130),nelem,lebon1,lebon2 double precision :: ltemp1(130),ltemp2(130) integer ibin,uMeteo,j c---------------------------------- c TESTS REAL satubf(ngrid,nlay),satuaf(ngrid,nlay) REAL res_out(ngrid,nlay) c------------------------------------------------------------------ IF (firstcall) THEN !============================================================= ! 0. Definition of the size grid !============================================================= c rad_cldco2 is the primary radius grid used for microphysics computation. c The grid spacing is computed assuming a constant volume ratio c between two consecutive bins; i.e. vrat_cld. c vrat_cld is determined from the boundary values of the size grid: c rmin_cld and rmax_cld. c The rb_cldco2 array contains the boundary values of each rad_cldco2 bin. c dr_cld is the width of each rad_cldco2 bin. c Volume ratio between two adjacent bins ! vrat_cld = log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3. ! vrat_cld = exp(vrat_cld) vrat_cld = log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3. vrat_cld = exp(vrat_cld) c write(*,*) "vrat_cld", vrat_cld rb_cldco2(1) = rbmin_cld rad_cldco2(1) = rmin_cld vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld ! vol_cld(1) = 4./3. * pi * rmin_cld*rmin_cld*rmin_cld do i=1,nbinco2_cld-1 rad_cldco2(i+1) = rad_cldco2(i) * vrat_cld**(1./3.) vol_cld(i+1) = vol_cld(i) * vrat_cld enddo do i=1,nbinco2_cld rb_cldco2(i+1)= ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) * & rad_cldco2(i) dr_cld(i) = rb_cldco2(i+1) - rb_cldco2(i) enddo rb_cldco2(nbinco2_cld+1) = rbmax_cld dr_cld(nbinco2_cld) = rb_cldco2(nbinco2_cld+1) - & rb_cldco2(nbinco2_cld) print*, ' ' print*,'Microphysics co2: size bin information:' print*,'i,rb_cldco2(i), rad_cldco2(i),dr_cld(i)' print*,'-----------------------------------' do i=1,nbinco2_cld write(*,'(i3,3x,3(e12.6,4x))') i,rb_cldco2(i), rad_cldco2(i), & dr_cld(i) enddo write(*,'(i3,3x,e12.6)') nbinco2_cld+1,rb_cldco2(nbinco2_cld+1) print*,'-----------------------------------' do i=1,nbinco2_cld+1 rb_cldco2(i) = log(rb_cldco2(i)) !! we save that so that it is not computed !! at each timestep and gridpoint enddo c Contact parameter of co2 ice on dst ( m=cos(theta) ) c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c mteta = 0.952 c mtetaco2 = 0.952 write(*,*) 'co2_param contact parameter:', mtetaco2 c Volume of a co2 molecule (m3) vo1 = m0co2 / dble(rho_ice_co2) ! m0co2 et non mco2 vo1co2=vo1 ! AJOUT JA c Variance of the ice and CCN distributions sigma_iceco2 = sqrt(log(1.+nuiceco2_sed)) sigma_ice = sqrt(log(1.+nuice_sed)) c write(*,*) 'Variance of ice & CCN distribs :', sigma_iceco2 c write(*,*) 'nuice for sedimentation:', nuiceco2_sed c write(*,*) 'Volume of a co2 molecule:', vo1co2 write(*,*) 'Variance of ice & CCN CO2 distribs :', sigma_iceco2 write(*,*) 'nuice for co2 ice sedimentation:', nuiceco2_sed write(*,*) 'Volume of a co2 molecule:', vo1co2 coeffh2o=0 if (co2useh2o) then write(*,*) write(*,*) "co2useh2o=.true. in callphys.def" write(*,*) "This means water ice particles can " write(*,*) "serve as CCN for CO2 microphysics" coeffh2o=1 endif meteo_ccn(:,:,:)=0. if (meteo_flux) then write(*,*) write(*,*) "meteo_flux=.true. in callphys.def" write(*,*) "meteoritic dust particles are available" write(*,*) "for co2 ice nucleation! " write(*,*) "Flux given by J. Plane (altitude,size bins)" ! Initialisation of the flux: it is constant and is it saved !We must interpolate the table to the GCM altitudes INQUIRE(FILE=datafile(1:LEN_TRIM(datafile))// & '/Meteo_flux_Plane.dat', EXIST=file_ok) IF (.not. file_ok) THEN write(*,*) 'file Meteo_flux_Plane.dat should be in ' & ,datafile STOP endif !used Variables open(newunit=uMeteo,file=trim(datafile)// & '/Meteo_flux_Plane.dat' & ,FORM='formatted') !13000 records (130 altitudes x 100 bin sizes) do i=1,130 do j=1,100 ! les mêmes 100 bins size que la distri nuclea: on touche pas read(uMeteo,'(F12.6)') meteo(i,j) enddo altitudes_meteo(i)=i ! On doit maintenant réinterpoler le tableau(130,100) sur les altitudes du GCM (nlay,100) enddo close(uMeteo) endif !of if meteo_flux l0=595594d0 l1=903.111d0 l2=-11.5959d0 l3=0.0528288d0 l4=-0.000103183d0 firstcall=.false. END IF c write(*,*) "max memdNN =",maxval(memdNNccn) !============================================================= ! 1. Initialisation !============================================================= !cste = 4*pi*rho_ice*ptimestep !not used for co2 flag_pourri=0 res_out(:,:) = 0 rice(:,:) = 1.e-8 riceco2(:,:) = 1.e-11 c Initialize the tendencies pdqcloudco2(1:ngrid,1:nlay,1:nq)=0. pdtcloudco2(1:ngrid,1:nlay)=0. c pt temperature layer; pdt dT.s-1 c pq traceur kg/kg; pdq tendance idem .s-1 zt(1:ngrid,1:nlay) = & pt(1:ngrid,1:nlay) + & pdt(1:ngrid,1:nlay) * ptimestep c call WRITEDIAGFI(ngrid,"Ztclouds", c & "Ztclouds",'K',3,zt) c call WRITEDIAGFI(ngrid,"pdtclouds", c & "pdtclouds",'K',3,pdt*ptimestep) zq(1:ngrid,1:nlay,1:nq) = & pq(1:ngrid,1:nlay,1:nq) + & pdq(1:ngrid,1:nlay,1:nq) * ptimestep WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 ) & zq(1:ngrid,1:nlay,1:nq) = 1.e-30 zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq) !============================================================= ! 2. Compute saturation !============================================================= dev2 = 1. / ( sqrt(2.) * sigma_iceco2 ) dev3 = 1. / ( sqrt(2.) * sigma_ice ) call co2sat(ngrid*nlay,zt,pplay,zqsat) !zqsat is psat(co2) !============================================================= ! Bonus: additional meteoritic particles for nucleation !============================================================= if (meteo_flux) then !altitude_meteo(130) !pzlev(ngrid,nlay) !meteo(130,100) !resultat: meteo_ccn(ngrid,nlay,100) do l=1,nlay do ig=1,ngrid ltemp1=abs(altitudes_meteo(:)-pzlev(ig,l)) ltemp2=abs(altitudes_meteo(:)-pzlev(ig,l+1)) lebon1=minloc(ltemp1,DIM=1) lebon2=minloc(ltemp2,DIM=1) nelem=lebon2-lebon1+1. mtemp(:)=0d0 !mtemp(100) : valeurs pour les 100bins do ibin=1,100 mtemp(ibin)=sum(meteo(lebon1:lebon2,ibin)) enddo meteo_ccn(ig,l,:)=mtemp(:)/nelem enddo enddo endif c Main loop over the GCM's grid DO l=1,nlay DO ig=1,ngrid c Get the partial pressure of co2 vapor and its saturation ratio pco2 = zq(ig,l,igcm_co2) * (mmean(ig,l)/44.01) * pplay(ig,l) c satu = zq(ig,l,igcm_co2) / zqsat(ig,l) satu = pco2 / zqsat(ig,l) !============================================================= ! 3. Nucleation !============================================================= rho_ice_co2T(ig,l)=1000.*(1.72391-2.53e-4*zt(ig,l) & -2.87e-6*zt(ig,l)*zt(ig,l)) vo2co2 = m0co2 / dble(rho_ice_co2T(ig,l)) rho_ice_co2=rho_ice_co2T(ig,l) IF ( satu .ge. 1d0 ) THEN ! if there is condensation rdust(ig,l)= zq(ig,l,igcm_dust_mass) & *0.75/pi/rho_dust & / zq(ig,l,igcm_dust_number) rdust(ig,l)= rdust(ig,l)**(1./3.) c write(*,*) "Improved2, l,Rdust = ",l,rdust(ig,l) rdust(ig,l)=max(1.e-9,rdust(ig,l)) c rdust(ig,l)=min(1.e-5,rdust(ig,l)) ! write(*,*) "Improved3,Rdust = ",rdust(ig,l) c Expand the dust moments into a binned distribution Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)+1.e-30 No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+1.e-30 Rn = rdust(ig,l) Rn = -log(Rn) Rm = Rn - 3. * sigma_iceco2*sigma_iceco2 n_derf = erf( (rb_cldco2(1)+Rn) *dev2) m_derf = erf( (rb_cldco2(1)+Rm) *dev2) do i = 1, nbinco2_cld n_aer(i) = -0.5 * No * n_derf !! this ith previously computed m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed n_derf = derf((rb_cldco2(i+1)+Rn) *dev2) m_derf = derf((rb_cldco2(i+1)+Rm) *dev2) n_aer(i) = n_aer(i) + 0.5 * No * n_derf + & meteo_ccn(ig,l,i) !Ajout meteo_ccn m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf m_aer2(i)=4./3.*pi*rho_dust & *n_aer(i)*rad_cldco2(i)*rad_cldco2(i) & *rad_cldco2(i) c write(*,*) "diff =",rad_cldco2(i),m_aer(i),m_aer2(i) enddo c write(*,*) "Bilan =",sum(m_aer),sum(m_aer2) c sumcheck = 0 c do i = 1, nbinco2_cld c sumcheck = sumcheck + n_aer(i) c enddo c sumcheck = abs(sumcheck/No - 1) c if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then c print*, "WARNING, No sumcheck PROBLEM" c print*, "sumcheck, No, rdust",sumcheck, No,rdust(ig,l) c print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l c print*, "Dust binned distribution", n_aer c STOP c endif c sumcheck = 0 c do i = 1, nbinco2_cld c sumcheck = sumcheck + m_aer(i) c enddo c sumcheck = abs(sumcheck/Mo - 1) c if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) c & then c print*, "WARNING, Mo sumcheck PROBLEM" c print*, "sumcheck, Mo",sumcheck, Mo c print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig,l c print*, "Dust binned distribution", m_aer c STOP c endif m_aer(:)=m_aer2(:) rccnh2o(ig,l)= zq(ig,l,igcm_ccn_mass) & *0.75/pi/rho_dust & / zq(ig,l,igcm_ccn_number) rice(ig,l)=( zq(ig,l,igcm_h2o_ice)*3.0/ & (4.0*rho_ice_co2*zq(ig,l,igcm_ccn_number) & *pi*tauscaling(ig)) +rccnh2o(ig,l)*rccnh2o(ig,l) & *rccnh2o(ig,l) )**(1.0/3.0) rhocloud(ig,l)=( zq(ig,l,igcm_h2o_ice)*rho_ice & +zq(ig,l,igcm_ccn_mass) *tauscaling(ig)*rho_dust) & / (zq(ig,l,igcm_h2o_ice)+zq(ig,l,igcm_ccn_mass) & *tauscaling(ig)) c call updaterice_micro( c & zq(ig,l,igcm_h2o_ice), ! ice mass c & zq(ig,l,igcm_ccn_mass), ! ccn mass c & zq(ig,l,igcm_ccn_number), ! ccn number c & tauscaling(ig),rice(ig,l),rhocloud(ig,l)) ! rice radius of CCN + H20 crystal c write(*,*) "Improved1 Rice=",rice(ig,l) rice(ig,l)=max(1.e-10,rice(ig,l)) c rice(ig,l)=min(1.e-5,rice(ig,l)) !write(*,*) "Improved2 Rice=",rice(ig,l) Mo = zq(ig,l,igcm_h2o_ice) + & zq(ig,l,igcm_ccn_mass)*tauscaling(ig)+1.e-30!*tauscaling(ig) ! & + 1.e-30 !Total mass of H20 crystals,CCN included No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30 Rn = rice(ig,l) Rn = -log(Rn) Rm = Rn - 3. * sigma_ice*sigma_ice n_derf = erf( (rb_cldco2(1)+Rn) *dev3) m_derf = erf( (rb_cldco2(1)+Rm) *dev3) do i = 1, nbinco2_cld n_aer_h2oice(i) = -0.5 * No * n_derf m_aer_h2oice(i) = -0.5 * Mo * m_derf n_derf = derf( (rb_cldco2(i+1)+Rn) *dev3) m_derf = derf( (rb_cldco2(i+1)+Rm) *dev3) n_aer_h2oice(i) = n_aer_h2oice(i) + 0.5 * No * n_derf m_aer_h2oice(i) = m_aer_h2oice(i) + 0.5 * Mo * m_derf rad_h2oice(i) = rad_cldco2(i) m_aer_h2oice2(i)=4./3.*pi*rhocloud(ig,l) & *n_aer_h2oice(i)*rad_h2oice(i)*rad_h2oice(i) & *rad_h2oice(i) c write(*,*) "before nuc, i,rad_h2o(i)= ",i,rad_cldco2(i) c & ,m_aer_h2oice(i),n_aer_h2oice(i) enddo c sumcheck = 0 c do i = 1, nbinco2_cld c sumcheck = sumcheck + n_aer_h2oice(i) c enddo c sumcheck = abs(sumcheck/No - 1) c if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then c print*, "WARNING, Noh2o sumcheck PROBLEM" c print*, "sumcheck, No, rice",sumcheck, No,rice(ig,l) c print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l c print*, "Dust binned distribution", n_aer_h2oice c STOP c endif c sumcheck = 0 c do i = 1, nbinco2_cld c sumcheck = sumcheck + m_aer_h2oice(i) c enddo c sumcheck = abs(sumcheck/Mo - 1) c if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) c & then c print*, "WARNING, Moh2o sumcheck PROBLEM" c print*, "sumcheck, Mo",sumcheck, Mo c print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig,l c print*, "Dust binned distribution", m_aer_h2oice c STOP c endif c Get the rates of nucleation call nucleaCO2(dble(pco2),zt(ig,l),dble(satu) & ,n_aer,rate,n_aer_h2oice & ,rad_h2oice,rateh2o,vo2co2) m_aer_h2oice(:)=m_aer_h2oice2(:) dN = 0. dM = 0. dNh2o = 0. dMh2o = 0. do i = 1, nbinco2_cld Proba=1.0-dexp(-1.*ptimestep*rate(i)) Probah2o=coeffh2o*(1.0-dexp(-1.*ptimestep*rateh2o(i))) dNh2o = dNh2o + n_aer_h2oice(i) * Probah2o dMh2o = dMh2o + m_aer_h2oice(i) * Probah2o dN = dN + n_aer(i) * Proba dM = dM + m_aer(i) * Proba enddo ! dM masse activée (kg) et dN nb particules par kg d'air dNN= dN/tauscaling(ig) dMM= dM/tauscaling(ig) dNN=min(dNN,zq(ig,l,igcm_dust_number)) dMM=min(dMM,zq(ig,l,igcm_dust_mass)) ! if (dNN .gt. 0) then ! write(*,*) "Nuclea dNN crees=",dNN ! write(*,*) "Nuclea dMM crees=",dMM ! endif dNNh2o=dNh2o/tauscaling(ig) dNNh2o=min(dNNh2o,zq(ig,l,igcm_ccn_number)) ratioh2o_ccn=1./(zq(ig,l,igcm_h2o_ice) & +zq(ig,l,igcm_ccn_mass)*tauscaling(ig)) dMh2o_ice=dMh2o*zq(ig,l,igcm_h2o_ice)*ratioh2o_ccn dMh2o_ccn=dMh2o*zq(ig,l,igcm_ccn_mass)* & tauscaling(ig)*ratioh2o_ccn dMh2o_ccn=dMh2o_ccn/tauscaling(ig) dMh2o_ccn=min(dMh2o_ccn,zq(ig,l,igcm_ccn_mass)) dMh2o_ice=min(dMh2o_ice,zq(ig,l,igcm_h2o_ice)) zq(ig,l,igcm_ccnco2_mass) = & zq(ig,l,igcm_ccnco2_mass) + dMM zq(ig,l,igcm_ccnco2_number) = & zq(ig,l,igcm_ccnco2_number) + dNN zq(ig,l,igcm_dust_mass)= zq(ig,l,igcm_dust_mass)-dMM zq(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number)-dNN c Update CCN for CO2 nucleating on H2O CCN : ! Warning: must keep memory of it zq(ig,l,igcm_ccnco2_mass) = & zq(ig,l,igcm_ccnco2_mass) + dMh2o_ice+dMh2o_ccn zq(ig,l,igcm_ccnco2_number) = & zq(ig,l,igcm_ccnco2_number) + dNNh2o zq(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)-dNNh2o zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)-dMh2o_ice zq(ig,l,igcm_ccn_mass)= zq(ig,l,igcm_ccn_mass)-dMh2o_ccn memdMMh2o(ig,l)=memdMMh2o(ig,l)+dMh2o_ice memdMMccn(ig,l)=memdMMccn(ig,l)+dMh2o_ccn memdNNccn(ig,l)=memdNNccn(ig,l)+dNNh2o ENDIF ! of is satu >1 !============================================================= ! 4. Ice growth: scheme for radius evolution !============================================================= c We trigger crystal growth if and only if there is at least one nuclei (N>1). c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible. IF (zq(ig,l,igcm_ccnco2_number)*tauscaling(ig) .ge. 1)THEN ! we trigger crystal growth c c write(*,*) "ccn number mass=",zq(ig,l,igcm_ccnco2_number), c & zq(ig,l,igcm_ccnco2_mass) c Niceco2 = zq(ig,l,igcm_co2_ice) c Qccnco2 = zq(ig,l,igcm_ccnco2_mass) c Nccnco2 = zq(ig,l,igcm_ccnco2_number) c call updaterice_microco2(Niceco2,Qccnco2,Nccnco2, c & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l)) c write(*,*) "updater rice=",riceco2(ig,l) rdust(ig,l)= zq(ig,l,igcm_ccnco2_mass) & *0.75/pi/rho_dust & / zq(ig,l,igcm_ccnco2_number) rdust(ig,l)= rdust(ig,l)**(1./3.) rdust(ig,l)=max(1.e-10,rdust(ig,l)) c rdust(ig,l)=min(5.e-6,rdust(ig,l)) riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/ & (4.0*rho_ice_co2*zq(ig,l,igcm_ccnco2_number) & *pi*tauscaling(ig)) +rdust(ig,l)*rdust(ig,l) & *rdust(ig,l) )**(1.0/3.0) c riceco2(ig,l)=max(1.e-10,riceco2(ig,l)) c riceco2(ig,l)=min(1.e-5,riceco2(ig,l)) ! WATCH OUT: CO2 nuclei is supposed to be dust ! only when deriving rhocloud (otherwise would need to keep info on water embedded in co2) - listo c write(*,*) "Rdust before growth = ",rdust(ig,l) c write(*,*) "Riceco2 before growth = ",riceco2(ig,l) !! Niceco2,Qccnco2,Nccnco2 c Niceco2 = zq(ig,l,igcm_co2_ice) c Qccnco2 = zq(ig,l,igcm_ccnco2_mass) c Nccnco2 = zq(ig,l,igcm_ccnco2_number) c call updaterice_microco2(Niceco2,Qccnco2,Nccnco2, c & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l)) c write(*,*) "Riceco2 before growth = ",riceco2(ig,l) c write(*,*) "rdust before growth = ",rdust(ig,l) c write(*,*) "co2ice before growth=",zq(ig,l,igcm_co2_ice) c write(*,*) "co2 before growth=",zq(ig,l,igcm_co2) c write(*,*) "pplay before growth=",pplay(ig,l) c write(*,*) "zt before growth =",zt(ig,l) c write(*,*) "Satu before growth=",satu c riceco2(ig,l)=max(riceco2(ig,l),rdust(ig,l)) No = zq(ig,l,igcm_ccnco2_number)* tauscaling(ig)+1.e-30 ! No nb de particules de poussieres mis à l'échelle pour donner une opacité optique c saturation at equilibrium c rice should not be too small, otherwise seq value is not valid c seq = exp(2.*sigco2*mco2 / (rho_ice_co2*rgp*zt(ig,l)* c & max(riceco2(ig,l),1.e-7))) !Exponant sans unité OK ccccccc Scheme of microphys. mass growth for CO2 call massflowrateCO2(pplay(ig,l),zt(ig,l), & satu,riceco2(ig,l),mmean(ig,l),Ic_rice) ! Mass transfer rate (kg/s) for a rice particle ! Ic_rice mass flux kg.s-1 <0 si croissance ! if (isnan(Ic_rice)) then Ic_rice=0. flag_pourri=1 endif c drsurdt=-1.0/(4.0*pi*riceco2(ig,l)* c & riceco2(ig,l)*rho_ice_co2)*Ic_rice dMice = No * Ic_rice*ptimestep ! Kg par kg d'air, <0 si croissance ! c write(*,*) "dMicev0 in improved = " , dMice if (dMice .ge. 0d0) then dMice = min(dMice,abs(zq(ig,l,igcm_co2_ice))) else dMice =-1.* min(abs(dMice),abs(zq(ig,l,igcm_co2))) endif c riceco2(ig,l)=riceco2(ig,l)+drsurdt*ptimestep c write(*,*) "riceco2+dr/dt = ", riceco2(ig,l) !write(*,*) "dMice in improved = " , dMice zq(ig,l,igcm_co2_ice) = zq(ig,l,igcm_co2_ice)-dMice zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)+dMice ! latent heat release >0 if growth i.e. if dMice <0 lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)**2 + & l3 * zt(ig,l)**3 + l4 * zt(ig,l)**4 !J.kg-1 c write(*,*) "CPP= ",cpp ! = 744.5 !Peut etre un probleme de signe ici! pdtcloudco2(ig,l)= -1.*dMice*lw/cpp/ptimestep ! kg par couche * J par kg /J par K / s = K par seconde !write(*,*) "pdtcloudco2 after growth = ",pdtcloudco2(ig,l) !write(*,*) "co2 after growth = ",zq(ig,l,igcm_co2) !write(*,*) "co2_ice after growth = ",zq(ig,l,igcm_co2_ice) !deltaT par condens/subli. qui remplace le dT du CO2 du newcondens pré-constantino !PDT should be in K/s. !============================================================= ! 5. Dust cores released, tendancies, latent heat, etc ... !============================================================= c if (abs(dMice) .ge. 1.e-10) then c write(*,*) "DIAG zq pdt",(zq(ig,l,igcm_co2_ice)- c & zq0(ig,l,igcm_co2_ice))/ptimestep,pdtcloudco2(ig,l) c endif ENDIF ! of if NCCN > 1 rdust(ig,l)= zq(ig,l,igcm_ccnco2_mass) & *0.75/pi/rho_dust & / zq(ig,l,igcm_ccnco2_number) rdust(ig,l)= rdust(ig,l)**(1./3.) rdust(ig,l)=max(1.e-9,rdust(ig,l)) c rdust(ig,l)=min(5.e-6,rdust(ig,l)) riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/ & (4.0*rho_ice_co2*zq(ig,l,igcm_ccnco2_number) & *pi*tauscaling(ig)) +rdust(ig,l)*rdust(ig,l) & *rdust(ig,l) )**(1.0/3.0) !Niceco2 = zq(ig,l,igcm_co2_ice) !Qccnco2 = zq(ig,l,igcm_ccnco2_mass) !Nccnco2 = zq(ig,l,igcm_ccnco2_number) c c call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2, c & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l)) c If there is no more co2 ice, all the ice particles sublimate, all the condensation nuclei are released: if ((zq(ig,l,igcm_co2_ice).le. 1.e-25)) then c & .or. flag_pourri .eq. 1 c & .or.(riceco2(ig,l) .le. rdust(ig,l)) c & .or. (l .ge.1 .and. l .le. 5) c & .or. (zq(ig,l,igcm_co2_ice) .ge. 0.1) lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)**2 + & l3 * zt(ig,l)**3 + l4 * zt(ig,l)**4 !J.kg-1 c write(*,*) "CPP= ",cpp ! = 744.5 pdtcloudco2(ig,l)=pdtcloudco2(ig,l)-1. & *zq(ig,l,igcm_co2_ice)*lw/cpp/ptimestep ! !On sublime tout ! c write(*,*) "Riceco2 improved before reset=",riceco2(ig,l) c write(*,*) "Niceco2 improved before reset=", c & zq(ig,l,igcm_co2_ice) c write(*,*) "Rdust improved before reset=",rdust(ig,l) if (memdMMccn(ig,l) .gt. 0) then zq(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass) & +memdMMccn(ig,l) endif if (memdMMh2o(ig,l) .gt. 0) then zq(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice) & +memdMMh2o(ig,l) endif if (memdNNccn(ig,l) .gt. 0) then zq(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number) & +memdNNccn(ig,l) endif c if (zq(ig,l,igcm_ccnco2_mass) .gt. 1.e-30) then zq(ig,l,igcm_dust_mass) = & zq(ig,l,igcm_dust_mass) & + zq(ig,l,igcm_ccnco2_mass)- & (memdMMh2o(ig,l)+memdMMccn(ig,l)) c endif c if (zq(ig,l,igcm_ccnco2_number) .gt. 1.e-30) then zq(ig,l,igcm_dust_number) = & zq(ig,l,igcm_dust_number) & + zq(ig,l,igcm_ccnco2_number)-memdNNccn(ig,l) c endif c if (zq(ig,l,igcm_co2_ice) .gt. 1.e-30) then zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2) & + zq(ig,l,igcm_co2_ice) c endif zq(ig,l,igcm_ccnco2_mass)=0. zq(ig,l,igcm_co2_ice)=0. zq(ig,l,igcm_ccnco2_number)=0. memdNNccn(ig,l)=0. memdMMh2o(ig,l)=0. memdMMccn(ig,l)=0. riceco2(ig,l)=0. endif ENDDO ! of ig loop ENDDO ! of nlayer loop ! Get cloud tendencies pdqcloudco2(1:ngrid,1:nlay,igcm_co2) = & (zq(1:ngrid,1:nlay,igcm_co2) - & zq0(1:ngrid,1:nlay,igcm_co2))/ptimestep pdqcloudco2(1:ngrid,1:nlay,igcm_co2_ice) = & (zq(1:ngrid,1:nlay,igcm_co2_ice) - & zq0(1:ngrid,1:nlay,igcm_co2_ice))/ptimestep pdqcloudco2(1:ngrid,1:nlay,igcm_h2o_ice) = & (zq(1:ngrid,1:nlay,igcm_h2o_ice) - & zq0(1:ngrid,1:nlay,igcm_h2o_ice))/ptimestep pdqcloudco2(1:ngrid,1:nlay,igcm_ccn_mass) = & (zq(1:ngrid,1:nlay,igcm_ccn_mass) - & zq0(1:ngrid,1:nlay,igcm_ccn_mass))/ptimestep pdqcloudco2(1:ngrid,1:nlay,igcm_ccn_number) = & (zq(1:ngrid,1:nlay,igcm_ccn_number) - & zq0(1:ngrid,1:nlay,igcm_ccn_number))/ptimestep pdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_mass) = & (zq(1:ngrid,1:nlay,igcm_ccnco2_mass) - & zq0(1:ngrid,1:nlay,igcm_ccnco2_mass))/ptimestep pdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_number) = & (zq(1:ngrid,1:nlay,igcm_ccnco2_number) - & zq0(1:ngrid,1:nlay,igcm_ccnco2_number))/ptimestep pdqcloudco2(1:ngrid,1:nlay,igcm_dust_mass) = & (zq(1:ngrid,1:nlay,igcm_dust_mass) - & zq0(1:ngrid,1:nlay,igcm_dust_mass))/ptimestep pdqcloudco2(1:ngrid,1:nlay,igcm_dust_number) = & (zq(1:ngrid,1:nlay,igcm_dust_number) - & zq0(1:ngrid,1:nlay,igcm_dust_number))/ptimestep return end