subroutine improvedCO2clouds(ngrid,nlay,ptimestep, & pplay,pt,pdt, & pq,pdq,pdqcloudco2,pdtcloudco2, & nq,tauscaling) ! 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 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 "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 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 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 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 l0,l1,l2,l3,l4 REAL cste REAL 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, n_derf, m_derf ! 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 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 DOUBLE PRECISION rate(nbinco2_cld) ! nucleation rate DOUBLE PRECISION rateh2o(nbinco2_cld) ! nucleation rate REAL seq DOUBLE PRECISION riceco2(ngrid,nlay) ! CO2Ice mean radius (m) 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 c Parameters of the size discretization c used by the microphysical scheme DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-9 ! Minimum radius (m) DOUBLE PRECISION, PARAMETER :: rmax_cld = 5.e-5 ! Maximum radius (m) DOUBLE PRECISION, PARAMETER :: rbmin_cld =0.099e-9 ! Minimum boundary radius (m) DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! 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 REAL sigma_iceco2 ! Variance of the ice and CCN distributions SAVE sigma_iceco2 DOUBLE PRECISION Niceco2,Qccnco2,Nccnco2 c---------------------------------- c TESTS INTEGER countcells LOGICAL test_flag ! flag for test/debuging outputs SAVE test_flag 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 c 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)) c write(*,*) 'Variance of ice & CCN distribs :', sigma_iceco2 c write(*,*) 'nuice for sedimentation:', nuiceco2_sed c write(*,*) 'Volume of a co2 molecule:', vo1co2 test_flag = .false. firstcall=.false. END IF !============================================================= ! 1. Initialisation !============================================================= !cste = 4*pi*rho_ice*ptimestep !not used for co2 res_out(:,:) = 0 rice(:,:) = 1.e-11 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 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 ) call co2sat(ngrid*nlay,zt,pplay,zqsat) !zqsat is psat(co2) countcells = 0 c Faire rice co2 update en n-1 puis a chaque microdt, mettre a jour riceco2 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 !============================================================= c call updaterccn(zq(ig,l,igcm_dust_mass), c & zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig)) IF ( satu .ge. 1d0 ) THEN ! if there is condensation write(*,*) write(*,*) "l, pco2, satu= ",l,pco2,satu c Masse_atm=mmean(ig,l)*1.e-3*pplay(ig,l)/rgp/zt(ig,l) !Kg par couche call updaterccn(zq(ig,l,igcm_dust_mass), & zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig)) write(*,*) "Improved, l,Rdust = ",l,rdust(ig,l) 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.) write(*,*) "Improved2, l,Rdust = ",l,rdust(ig,l) rdust(ig,l)=max(1.e-9,rdust(ig,l)) rdust(ig,l)=min(2e-6,rdust(ig,l)) write(*,*) "Improved3, l,Rdust = ",l,rdust(ig,l) c Expand the dust moments into a binned distribution Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig) No = zq(ig,l,igcm_dust_number)* tauscaling(ig) write(*,*) "dust number, mass = ", & zq(ig,l,igcm_dust_number)* tauscaling(ig), & zq(ig,l,igcm_dust_mass)* tauscaling(ig) c write(*,*) "No, Mo = ",No, Mo 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 m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf c write(*,*) "i, rb_cldco2(i) = ",i, rb_cldco2(i),n_aer(i) enddo sumcheck = 0 do i = 1, nbinco2_cld sumcheck = sumcheck + n_aer(i) enddo sumcheck = abs(sumcheck/No - 1) if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then print*, "WARNING, No sumcheck PROBLEM" print*, "sumcheck, No",sumcheck, No print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l print*, "Dust binned distribution", n_aer STOP endif sumcheck = 0 do i = 1, nbinco2_cld sumcheck = sumcheck + m_aer(i) enddo sumcheck = abs(sumcheck/Mo - 1) if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) & then print*, "WARNING, Mo sumcheck PROBLEM" print*, "sumcheck, Mo",sumcheck, Mo print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig,l print*, "Dust binned distribution", m_aer STOP endif c Expand the water ice moments into a binned distribution c For now the radius grid's bound are same as for dust c min=100 nm and max=10microns c might need a change if rice (water) is large (but how large?) - listo Mo = zq(ig,l,igcm_h2o_ice)* tauscaling(ig) + 1.e-30 No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30 Rn = rice(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_h2oice(i) = -0.5 * No * n_derf !! this ith previously computed m_aer_h2oice(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_h2oice(i) = n_aer(i) + 0.5 * No * n_derf ! vector not really needed - temp var - listo m_aer_h2oice(i) = m_aer(i) + 0.5 * Mo * m_derf ! vector not really needed - temp var rad_h2oice(i) = ((m_aer_h2oice(i)/rho_ice/ & n_aer_h2oice(i) + vol_cld(i))) & *0.75/pi**(1./3) c write(*,*) "before nuc, i,rad_h2o(i)= ",i,rad_h2oice(i) c & ,m_aer(i),n_aer(i) enddo c Get the rates of nucleation call nucleaCO2(dble(pco2),zt(ig,l),dble(satu) & ,n_aer,rate,n_aer_h2oice & ,rad_h2oice,rateh2o) ! regarder rateh20, et mettre = 0 si non nul pour le moment dN = 0. dM = 0. dNh2o = 0. dMh2o = 0. do i = 1, nbinco2_cld !n_aer(i) = n_aer(i)/( 1. + (rate(i)+rateh2o(i))*ptimestep) !m_aer(i) = m_aer(i)/( 1. + (rate(i)+rateh2o(i))*ptimestep) Proba=1.0-dexp(-1.*ptimestep*rate(i)) c dNh2o = dNh2o + n_aer_h2oice(i) * rateh2o(i) c & * ptimestep c dMh2o = dMh2o + m_aer_h2oice(i) * rateh2o(i) c & *ptimestep dN = dN + n_aer(i) * Proba dM = dM + m_aer(i) * Proba c write(*,*) "i, dNi, dN= ",i,n_aer(i)*Proba,dN enddo c do i=1,nbinco2_cld c write(*,*) "i n_aer m_aer = ",i,n_aer(i),m_aer(i) c enddo ! dM masse activée (kg) et dN nb particules par kg d'air c Update Dust particles c For CO2 ice : no subtraction from dust (neither for water ice particles) ! zq(ig,l,igcm_dust_mass) = ! & zq(ig,l,igcm_dust_mass) - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10) ! zq(ig,l,igcm_dust_number) = ! & zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10) c write(*,*) " nuclea dM = ",dM/tauscaling(ig), c & " nuclea dN = ", dN/tauscaling(ig) dNN= dN/tauscaling(ig) dMM= dM/tauscaling(ig) dNN=min(dNN,abs(zq0(ig,l,igcm_dust_number))) dMM=min(dMM,zq0(ig,l,igcm_dust_mass)) c Update CCNs for CO2 crystals ! WARNING dM dMh2o, interaction nuages eau-co2 -- h20 set to 0 for now 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 + enlever les CCN a la distri de dust write(*,*) "new dust_mass, number =", & zq(ig,l,igcm_dust_mass)* tauscaling(ig), & zq(ig,l,igcm_dust_number)*tauscaling(ig) write(*,*) "new ccn mass, number =", & zq(ig,l,igcm_ccnco2_mass)* tauscaling(ig) & ,zq(ig,l,igcm_ccnco2_number)*tauscaling(ig) 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. c IF ( zq(ig,l,igcm_ccnco2_number)*tauscaling(ig).ge. 1.0) THEN IF (zq(ig,l,igcm_ccnco2_number)*tauscaling(ig) .ge. threshJA) & THEN ! we trigger crystal growth call updaterccn(zq(ig,l,igcm_dust_mass), & zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig)) 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)) rdust(ig,l)=min(2e-6,rdust(ig,l)) ! rdust(ig,l)=1.e-7 IF (zq(ig,l,igcm_ccnco2_mass) .lt. 0. .or. & zq(ig,l,igcm_ccnco2_number) .lt. 0. .or. & zq(ig,l,igcm_dust_mass) .lt. 0. .or. & zq(ig,l,igcm_dust_number) .lt. 0. ) THEN write(*,*) "before growth CCN N,M = " & ,zq(ig,l,igcm_ccnco2_number)*tauscaling(ig) & ,zq(ig,l,igcm_ccnco2_mass)*tauscaling(ig) write(*,*) "before growth dust number mass = ", & zq(ig,l,igcm_dust_number)*tauscaling(ig), & zq(ig,l,igcm_dust_mass)*tauscaling(ig) STOP END IF c write(*,*) "reff dN = ",reff,dN c reff=reff/dble(dN) c if (zq(ig,l,igcm_co2_ice) .le. 1.e-20) then c riceco2(ig,l)=reff c endif c write(*,*) "Rdust in improved = ",rdust(ig,l) riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/ & (4.0*rho_ice_co2*pi*zq(ig,l,igcm_ccnco2_number) & *tauscaling(ig)) +rdust(ig,l)*rdust(ig,l) & *rdust(ig,l) )**(1.0/3.0) ! 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 write(*,*) "Rdust before growth = ",rdust(ig,l) write(*,*) "Riceco2 before growth = ",riceco2(ig,l) !! Niceco2,Qccnco2,Nccnco2 Niceco2 = zq(ig,l,igcm_co2_ice) Qccnco2 = zq(ig,l,igcm_ccnco2_mass) Nccnco2 = zq(ig,l,igcm_ccnco2_number) call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2, & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l)) write(*,*) "Riceco2 update before growth = ",riceco2(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 ! drsurdt=-1.0/(4.0*pi*riceco2(ig,l)* & riceco2(ig,l)*rho_ice_co2)*Ic_rice dMice = No * Ic_rice * ptimestep ! Kg par kg d'air, <0 si croissance ! write(*,*) "dMicev0 in improved = " , dMice if (dMice .gt. 0) dMice = min(dMice,zq0(ig,l,igcm_co2_ice)) if (dMice .lt. 0) dMice = max(dMice,-1.*zq0(ig,l,igcm_co2)) 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 c write(*,*) "zq co2 ice = ", zq(ig,l,igcm_co2_ice) countcells = countcells + 1 riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/ & (4.0*rho_ice_co2*pi*zq(ig,l,igcm_ccnco2_number) & *tauscaling(ig)) +rdust(ig,l)*rdust(ig,l) & *rdust(ig,l) )**(1.0/3.0) write(*,*) "new riceco2 = ",riceco2(ig,l) !! Niceco2,Qccnco2,Nccnco2 Niceco2 = zq(ig,l,igcm_co2_ice) Qccnco2 = zq(ig,l,igcm_ccnco2_mass) Nccnco2 = zq(ig,l,igcm_ccnco2_number) call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2, & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l)) write(*,*) "new riceco2 updaterad= ",riceco2(ig,l) ! latent heat release l0=595594. l1=903.111 l2=-11.5959 l3=0.0528288 l4=-0.000103183 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)= dMice*lw/cpp/ptimestep ! kg par couche * J par kg /J par K / s = K par seconde c write(*,*) "pdtcloudco2 = ",pdtcloudco2(ig,l) c Peut etre -1*dMice? !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 all the ice particles sublimate, all the condensation c nuclei are released: c !!! with CO2 ice nuclei: dust/H2O nuclei are not released because c they were not subtracted to dust_number c Their counter is just set to "0". c (see end of section 3.) : On peut les enlever à dust c interaction ho-co2 ici, dans la mise a jour des traceurs WARNING reflechir ENDIF !of if Nccn>1 c if (riceco2(ig,l) .lt. 1.e-9) then if (zq(ig,l,igcm_co2_ice).le.1.e-20 .or. & riceco2(ig,l) .lt. 1.e-9 .or. riceco2(ig,l) & .ge. 4.999e-4) then ! Reverser les ccn libérés dans les h2o ou dust? c c ICI: Co2 ice devient vapeur zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2) & + zq(ig,l,igcm_co2_ice) zq(ig,l,igcm_dust_mass) = & zq(ig,l,igcm_dust_mass) & + zq(ig,l,igcm_ccnco2_mass) zq(ig,l,igcm_dust_number) & = zq(ig,l,igcm_dust_number) & + zq(ig,l,igcm_ccnco2_number) c c CCNs zq(ig,l,igcm_ccnco2_mass) = 0. zq(ig,l,igcm_ccnco2_number) =0. zq(ig,l,igcm_co2_ice) = 0. riceco2(ig,l)=0. endif c write(*,*) "zq co2 end imp= ", zq(i g,l,igcm_co2_ice),satu 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 c do l=1,nlay c write(*,*) "end imp",pdqcloudco2(1,l,igcm_co2), c & pdqcloudco2(1,l,igcm_co2_ice) c enddo 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 if (scavenging) then 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 endif c call WRITEdiagfi(ngrid,"Improvedb","co2 traceur","kg/kg",1, c & zq0(1,:,igcm_co2_ice) ) c call WRITEdiagfi(ngrid,"Improveda","co2 traceur","kg/kg",1, c & zq(1,:,igcm_co2_ice) ) c call WRITEdiagfi(ngrid,"satuco2","satu co2 in improved","kg/kg",1, c & satu) !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS IF (test_flag) then ! error2d(:) = 0. DO l=1,nlay DO ig=1,ngrid ! error2d(ig) = max(abs(error_out(ig,l)),error2d(ig)) satubf(ig,l) = zq0(ig,l,igcm_co2)/zqsat(ig,l) ! att. if zqsat=mmr or psat satuaf(ig,l) = zq(ig,l,igcm_co2)/zqsat(ig,l) ENDDO ENDDO write(*,*) 'count is ',countcells, ' i.e. ', & countcells*100/(nlay*ngrid), '% for microphys computation' ! IF (ngrid.ne.1) THEN ! 3D ! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3, ! & satu_out) ! call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3, ! & dM_out) ! call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3, ! & dN_out) ! call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2, ! & error2d) ! call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3, ! & zqsat) ! ENDIF ! IF (ngrid.eq.1) THEN ! 1D ! call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1, ! & error_out) ! call WRITEdiagfi(ngrid,"resist","resistance","s/m2",1, ! & res_out) call WRITEdiagfi(ngrid,"satu_bf","satu before","kg/kg",1, & satubf) call WRITEdiagfi(ngrid,"satu_af","satu after","kg/kg",1, & satuaf) call WRITEdiagfi(ngrid,"vapbf","CO2vap before","kg/kg",1, & zq0(1,1,igcm_co2)) call WRITEdiagfi(ngrid,"vapaf","CO2vap after","kg/kg",1, & zq(1,1,igcm_co2)) call WRITEdiagfi(ngrid,"icebf","CO2ice before","kg/kg",1, & zq0(1,1,igcm_co2_ice)) call WRITEdiagfi(ngrid,"iceaf","CO2ice after","kg/kg",1, & zq(1,1,igcm_co2_ice)) call WRITEdiagfi(ngrid,"ccnbf","ccn before","/kg",1, & zq0(1,1,igcm_ccnco2_number)) call WRITEdiagfi(ngrid,"ccnaf","ccn after","/kg",1, & zq(1,1,igcm_ccnco2_number)) c call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1, c & gr_out) c call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1, c & rate_out) c call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1, c & dM_out) c call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1, c & dN_out) c call WRITEdiagfi(ngrid,"zqsat","p vap sat","kg/kg",1, c & zqsat) ! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1, ! & satu_out) ! call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1, ! & rdust) ! call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1, ! & rsedcloud) ! call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1, ! & rhocloud) ! ENDIF ENDIF ! endif test_flag !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0 c It is an analytical solution to the ice radius growth equation, c with the approximation of a constant 'reduced' cunningham correction factor c (lambda in growthrate.F) taken at radius req instead of rice cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine phi(rice,req,coeff1,coeff2,time) c c implicit none c c ! inputs c real rice ! ice radius c real req ! ice radius at equilibirum c real coeff1 ! coeff for the log c real coeff2 ! coeff for the arctan c c ! output c real time c c !local c real var c c ! 1.73205 is sqrt(3) c c var = max( c & abs(rice-req) / sqrt(rice*rice + rice*req + req*req),1e-30) c c time = c & coeff1 * c & log( var ) c & + coeff2 * 1.73205 * c & atan( (2*rice+req) / (1.73205*req) ) c c return c end