SUBROUTINE co2cloud(ngrid,nlay,ptimestep, & pplev,pplay,pdpsrf,pzlay,pt,pdt, & pq,pdq,pdqcloudco2,pdtcloudco2, & nq,tau,tauscaling,rdust,rice,riceco2,nuice, & rsedcloudco2,rhocloudco2, & rsedcloud,rhocloud,pzlev,pdqs_sedco2) ! to use 'getin' use dimradmars_mod, only: naerkind USE comcstfi_h USE ioipsl_getincom USE updaterad use conc_mod, only: mmean use tracer_mod, only: nqmx, igcm_co2, igcm_co2_ice, & igcm_dust_mass, igcm_dust_number,igcm_h2o_ice, & igcm_ccn_mass,igcm_ccn_number, & igcm_ccnco2_mass, igcm_ccnco2_number, & rho_dust, nuiceco2_sed, nuiceco2_ref, & rho_ice_co2,r3n_q,rho_ice,nuice_sed IMPLICIT NONE #include "datafile.h" #include "callkeys.h" #include "microphys.h" c======================================================================= c CO2 clouds formation c c There is a time loop specific to cloud formation c due to timescales smaller than the GCM integration timestep. c microphysics subroutine is improvedCO2clouds.F c c The co2 clouds tracers (co2_ice, ccn mass and concentration) are c sedimented at each microtimestep. pdqs_sedco2 keeps track of the c CO2 flux at the surface c c Authors: 09/2016 Joachim Audouard & Constantino Listowski c Adaptation of the water ice clouds scheme (with specific microphysics) c of Montmessin, Navarro & al. c c 07/2017 J.Audouard c Several logicals can be set to .true. in callphys.def c co2clouds=.true. call this routine c co2useh2o=.true. allow the use of water ice particles as CCN for CO2 c meteo_flux=.true. supply meteoritic particles c CLFvaryingCO2=.true. allows a subgrid temperature distribution c of amplitude spantCO2(=integer in callphys.def) c c======================================================================= c----------------------------------------------------------------------- c declarations: c ------------- c Inputs: c ------ INTEGER ngrid,nlay INTEGER nq ! nombre de traceurs REAL ptimestep ! pas de temps physique (s) REAL pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa) REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa) REAL pdpsrf(ngrid) ! tendence surf pressure REAL pzlay(ngrid,nlay) ! altitude at the middle of the layers REAL pt(ngrid,nlay) ! temperature at the middle of the layers (K) REAL pdt(ngrid,nlay) ! tendence temperature des autres param. real,intent(in) :: pzlev(ngrid,nlay+1) ! altitude at the boundaries of the layers real pq(ngrid,nlay,nq) ! traceur (kg/kg) real pdq(ngrid,nlay,nq) ! tendance avant condensation (kg/kg.s-1) real rice(ngrid,nlay) ! Water Ice mass mean radius (m) ! used for nucleation of CO2 on ice-coated ccns DOUBLE PRECISION rho_ice_co2T(ngrid,nlay) REAL tau(ngrid,naerkind) ! Column dust optical depth at each point REAL tauscaling(ngrid) ! Convertion factor for dust amount real rdust(ngrid,nlay) ! Dust geometric mean radius (m) c Outputs: c ------- real pdqcloudco2(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1) REAL pdtcloudco2(ngrid,nlay) ! tendence temperature due ! a la chaleur latente DOUBLE PRECISION riceco2(ngrid,nlay) ! Ice mass mean radius (m) ! (r_c in montmessin_2004) REAL nuice(ngrid,nlay) ! Estimated effective variance ! of the size distribution real rsedcloudco2(ngrid,nlay) ! Cloud sedimentation radius real rhocloudco2(ngrid,nlay) ! Cloud density (kg.m-3) real rhocloudco2t(ngrid,nlay) ! Cloud density (kg.m-3) real pdqs_sedco2(ngrid) ! CO2 flux at the surface c local: c ------ !water real rsedcloud(ngrid,nlay) ! Cloud sedimentation radius real rhocloud(ngrid,nlay) ! Cloud density (kg.m-3) ! for ice radius computation REAl ccntyp ! for time loop INTEGER microstep ! time subsampling step variable INTEGER imicro ! time subsampling for coupled water microphysics & sedimentation SAVE imicro REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation SAVE microtimestep ! tendency given by clouds (inside the micro loop) REAL subpdqcloudco2(ngrid,nlay,nq) ! cf. pdqcloud REAL subpdtcloudco2(ngrid,nlay) ! cf. pdtcloud ! global tendency (clouds+physics) REAL subpdq(ngrid,nlay,nq) ! cf. pdqcloud REAL subpdt(ngrid,nlay) ! cf. pdtcloud real wq(ngrid,nlay+1) ! ! displaced tracer mass (kg.m-2) during microtimestep because sedim (?/m-2) REAL satuco2(ngrid,nlay) ! co2 satu ratio for output REAL zqsatco2(ngrid,nlay) ! saturation co2 INTEGER iq,ig,l,i LOGICAL,SAVE :: firstcall=.true. DOUBLE PRECISION Nccnco2, Niceco2,Nco2,mdustJA,ndustJA DOUBLE PRECISION Qccnco2 real :: beta real epaisseur (ngrid,nlay) ! Layer thickness (m) real masse (ngrid,nlay) ! Layer mass (kg.m-2) double precision diff,diff0 real tempo_traceur_t(ngrid,nlay) real tempo_traceurs(ngrid,nlay,nq) real sav_trac(ngrid,nlay,nq) real pdqsed(ngrid,nlay,nq) REAL lw !Latent heat of sublimation (J.kg-1) REAL,save :: l0,l1,l2,l3,l4 !Coeffs for lw DOUBLE PRECISION,allocatable,save :: memdMMccn(:,:) ! Nb particules intégré DOUBLE PRECISION,allocatable,save :: memdMMh2o(:,:) DOUBLE PRECISION,allocatable,save :: memdNNccn(:,:) DOUBLE PRECISION :: myT ! What we need for Qext reading and tau computation 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, 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 boundary radius (m) DOUBLE PRECISION, PARAMETER :: rbmax_cld = 5.e-6 ! Maximum boundary radius (m) 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) REAL sigma_iceco2 ! Variance of the ice and CCN distributions logical :: file_ok double precision :: radv(10000),Qextv1mic(10000) double precision :: Qext1bins(100),Qtemp double precision :: ltemp1(10000),ltemp2(10000) integer :: nelem,lebon1,lebon2,uQext save Qext1bins character(len=100) scanline DOUBLE PRECISION n_aer(nbinco2_cld),Rn,No,n_derf,dev2 DOUBLE PRECISION Qext1bins2(ngrid,nlay) DOUBLE PRECISION tau1mic(ngrid) !co2 ice column opacity at 1µm ! For sub grid T distribution REAL zt(ngrid,nlay) ! local value of temperature REAL :: zq(ngrid, nlay,nq) REAL :: tcond(ngrid,nlay) REAL :: zqvap(ngrid,nlay) REAL :: zqice(ngrid,nlay) REAL :: spant,zdelt ! delta T for the temperature distribution REAL :: zteff(ngrid, nlay)! effective temperature in the cloud,neb REAL :: pqeff(ngrid, nlay, nq)! effective tracers quantities in the cloud REAL :: cloudfrac(ngrid,nlay) ! cloud fraction REAL :: mincloud ! min cloud frac c logical :: CLFvaryingCO2 c ** un petit test de coherence c -------------------------- IF (firstcall) THEN if (nq.gt.nqmx) then write(*,*) 'stop in co2cloud (nq.gt.nqmx)!' write(*,*) 'nq=',nq,' nqmx=',nqmx stop endif write(*,*) "co2cloud.F: rho_ice_co2 = ",rho_ice_co2 write(*,*) "co2cloud: igcm_co2=",igcm_co2 write(*,*) " igcm_co2_ice=",igcm_co2_ice write(*,*) "time subsampling for microphysic ?" #ifdef MESOSCALE imicro = 2 #else imicro = 30 #endif call getin("imicro",imicro) write(*,*)"imicro = ",imicro microtimestep = ptimestep/real(imicro) write(*,*)"Physical timestep is",ptimestep write(*,*)"CO2 Microphysics timestep is",microtimestep ! Values for latent heat computation l0=595594d0 l1=903.111d0 l2=-11.5959d0 l3=0.0528288d0 l4=-0.000103183d0 c$$$ c$$$ if (meteo_flux_number .ne. 0) then c$$$ c$$$ write(*,*) "Warning ! Meteoritic flux of dust is turned on" c$$$ write(*,*) "Dust mass = ",meteo_flux_mass c$$$ write(*,*) "Dust number = ",meteo_flux_number c$$$ write(*,*) "Are added at the z-level (km)",meteo_alt c$$$ write(*,*) "Every timestep in co2cloud.F" c$$$ endif if (.not. allocated(memdMMccn)) allocate(memdMMccn(ngrid,nlay)) if (.not. allocated(memdNNccn)) allocate(memdNNccn(ngrid,nlay)) if (.not. allocated(memdMMh2o)) allocate(memdMMh2o(ngrid,nlay)) memdMMccn(:,:)=0. memdMMh2o(:,:)=0. memdNNccn(:,:)=0. c Compute the size bins of the distribution of CO2 ice particles c --> used for opacity calculations 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. sigma_iceco2 = sqrt(log(1.+nuiceco2_sed)) c Volume ratio between two adjacent bins ! vrat_cld vrat_cld = log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3. vrat_cld = exp(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 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) c read the Qext values INQUIRE(FILE=datafile(1:LEN_TRIM(datafile))// & '/optprop_co2ice_1mic.dat', EXIST=file_ok) IF (.not. file_ok) THEN write(*,*) 'file optprop_co2ice_1mic.dat should be in ' & ,datafile STOP endif open(newunit=uQext,file=trim(datafile)// & '/optprop_co2ice_1mic.dat' & ,FORM='formatted') read(uQext,*) !skip 1 line do i=1,10000 read(uQext,'(E11.5)') radv(i) enddo read(uQext,*) !skip 1 line do i=1,10000 read(uQext,'(E11.5)') Qextv1mic(i) enddo close(uQext) c innterpol the Qext values !rice_out=rad_cldco2 do i=1,nbinco2_cld ltemp1=abs(radv(:)-rb_cldco2(i)) ltemp2=abs(radv(:)-rb_cldco2(i+1)) lebon1=minloc(ltemp1,DIM=1) lebon2=minloc(ltemp2,DIM=1) nelem=lebon2-lebon1+1. Qtemp=0d0 do l=0,nelem Qtemp=Qtemp+Qextv1mic(lebon1+l) !mean value in the interval enddo Qtemp=Qtemp/nelem Qext1bins(i)=Qtemp enddo Qext1bins(:)=Qext1bins(:)*rad_cldco2(:)*rad_cldco2(:)*pi ! The actuall tau computation and output is performed in co2cloud.F print*,'--------------------------------------------' print*,'Microphysics co2: size bin-Qext information:' print*,' i, rad_cldco2(i), Qext1bins(i)' do i=1,nbinco2_cld write(*,'(i3,3x,3(e12.6,4x))') i, rad_cldco2(i), & Qext1bins(i) enddo 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 if (CLFvaryingCO2) then write(*,*) write(*,*) "CLFvaryingCO2 is set to true is callphys.def" write(*,*) "The temperature field is enlarged to +/-",spantCO2 write(*,*) "for the CO2 microphysics " endif firstcall=.false. ENDIF ! of IF (firstcall) c-----Initialization dev2 = 1. / ( sqrt(2.) * sigma_iceco2 ) beta=0.85 subpdq(1:ngrid,1:nlay,1:nq) = 0 subpdt(1:ngrid,1:nlay) = 0 subpdqcloudco2(1:ngrid,1:nlay,1:nq) = 0 subpdtcloudco2(1:ngrid,1:nlay) = 0 c$$$ c$$$ call WRITEDIAGFI(ngrid,"pzlay","pzlay","km",3, c$$$ & pzlay) c$$$ call WRITEDIAGFI(ngrid,"pplay","pplay","Pa",3, c$$$ & pplay) wq(:,:)=0 ! default value if no ice rhocloudco2(1:ngrid,1:nlay) = rho_dust rhocloudco2t(1:ngrid,1:nlay) = rho_dust epaisseur(1:ngrid,1:nlay)=0 masse(1:ngrid,1:nlay)=0 sav_trac(1:ngrid,1:nlay,1:nq)=0 pdqsed(1:ngrid,1:nlay,1:nq)=0 do l=1,nlay do ig=1, ngrid masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l) enddo enddo !CLFvaryingCO2=.true. c------------------------------------------------------------------- c 0. Representation of sub-grid water ice clouds c------------------ IF (CLFvaryingCO2) THEN spant=spantCO2 ! delta T for the temprature distribution mincloud=0.1 ! min cloudfrac when there is ice ! IF (flagcloudco2) THEN ! write(*,*) "CLFCO2 Delta T", spant ! write(*,*) "CLFCO2 mincloud", mincloud ! flagcloudco2=.false. ! END IF c-----Tendencies DO l=1,nlay DO ig=1,ngrid zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep ENDDO ENDDO DO l=1,nlay DO ig=1,ngrid DO iq=1,nq zq(ig,l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep ENDDO ENDDO ENDDO zqvap=zq(:,:,igcm_co2) zqice=zq(:,:,igcm_co2_ice) CALL tcondco2(ngrid,nlay,pplay,zqvap,tcond) ! A tester: CALL tcondco2(ngrid,nlay,pplay,zqvap,tcond) zdelt=spant DO l=1,nlay DO ig=1,ngrid IF (tcond(ig,l) .ge. (zt(ig,l)+zdelt)) THEN !Toute la fraction est saturée zteff(ig,l)=zt(ig,l) cloudfrac(ig,l)=1. ELSE IF (tcond(ig,l) .le. (zt(ig,l)-zdelt)) THEN !Rien n'est saturé zteff(ig,l)=zt(ig,l)-zdelt cloudfrac(ig,l)=mincloud ELSE cloudfrac(ig,l)=(tcond(ig,l)-zt(ig,l)+zdelt)/ & (2.0*zdelt) zteff(ig,l)=(tcond(ig,l)+zt(ig,l)-zdelt)/2. !Ja Temperature moyenne de la fraction nuageuse END IF zteff(ig,l)=zteff(ig,l)-pdt(ig,l)*ptimestep IF (cloudfrac(ig,l).le. mincloud) THEN cloudfrac(ig,l)=mincloud ELSE IF (cloudfrac(ig,l).ge. 1) THEN cloudfrac(ig,l)=1. END IF ENDDO ENDDO ! Totalcloud frac of the column missing here c----------------------- c-----No sub-grid cloud representation (CLFvarying=false) ELSE DO l=1,nlay DO ig=1,ngrid zteff(ig,l)=pt(ig,l) END DO END DO END IF ! end if (CLFvaryingco2)-------------------------------------------- c 1. Tendencies: c------------------ c------ Effective tracers quantities in the cloud fraction IF (CLFvaryingCO2) THEN pqeff(:,:,:)=pq(:,:,:) ! prevent from buggs (A. Pottier) c pqeff(:,:,igcm_ccn_mass) =pq(:,:,igcm_ccn_mass)/ c & cloudfrac(:,:) c pqeff(:,:,igcm_ccn_number)= c & pq(:,:,igcm_ccn_number)/cloudfrac(:,:) c pqeff(:,:,igcm_h2o_ice)= pq(:,:,igcm_h2o_ice)/ c & cloudfrac(:,:) pqeff(:,:,igcm_ccnco2_mass) =pq(:,:,igcm_ccnco2_mass)/ & cloudfrac(:,:) pqeff(:,:,igcm_ccnco2_number)= & pq(:,:,igcm_ccnco2_number)/cloudfrac(:,:) pqeff(:,:,igcm_co2_ice)= pq(:,:,igcm_co2_ice)/ & cloudfrac(:,:) ELSE pqeff(:,:,:)=pq(:,:,:) c pqeff(:,:,igcm_ccn_mass)= pq(:,:,igcm_ccn_mass) c pqeff(:,:,igcm_ccn_number)= pq(:,:,igcm_ccn_number) c pqeff(:,:,igcm_h2o_ice)= pq(:,:,igcm_h2o_ice) c pqeff(:,:,igcm_ccnco2_mass)= pq(:,:,igcm_ccnco2_mass) c pqeff(:,:,igcm_ccnco2_number)= pq(:,:,igcm_ccnco2_number) c pqeff(:,:,igcm_co2_ice)= pq(:,:,igcm_co2_ice) END IF tempo_traceur_t(1:ngrid,1:nlay)=zteff(1:ngrid,1:nlay) tempo_traceurs(1:ngrid,1:nlay,1:nq)=pqeff(1:ngrid,1:nlay,1:nq) c------------------------------------------------------------------ c Time subsampling for microphysics c------------------------------------------------------------------ DO microstep=1,imicro c------ Temperature tendency subpdt ! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt ! If imicro=1 subpdt is the same as pdt DO l=1,nlay DO ig=1,ngrid subpdt(ig,l) = subpdt(ig,l) & + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry subpdq(ig,l,igcm_dust_mass) = & subpdq(ig,l,igcm_dust_mass) & + pdq(ig,l,igcm_dust_mass) subpdq(ig,l,igcm_dust_number) = & subpdq(ig,l,igcm_dust_number) & + pdq(ig,l,igcm_dust_number) subpdq(ig,l,igcm_ccnco2_mass) = & subpdq(ig,l,igcm_ccnco2_mass) & + pdq(ig,l,igcm_ccnco2_mass) subpdq(ig,l,igcm_ccnco2_number) = & subpdq(ig,l,igcm_ccnco2_number) & + pdq(ig,l,igcm_ccnco2_number) subpdq(ig,l,igcm_co2_ice) = & subpdq(ig,l,igcm_co2_ice) & + pdq(ig,l,igcm_co2_ice) subpdq(ig,l,igcm_co2) = & subpdq(ig,l,igcm_co2) & + pdq(ig,l,igcm_co2) subpdq(ig,l,igcm_h2o_ice) = & subpdq(ig,l,igcm_h2o_ice) & + pdq(ig,l,igcm_h2o_ice) subpdq(ig,l,igcm_ccn_mass) = & subpdq(ig,l,igcm_ccn_mass) & + pdq(ig,l,igcm_ccn_mass) subpdq(ig,l,igcm_ccn_number) = & subpdq(ig,l,igcm_ccn_number) & + pdq(ig,l,igcm_ccn_number) ENDDO ENDDO c add meteoritic flux of dust (old version) cNew version (John Plane values) are added in improvedCO2clouds !convert meteo_alt (in km) to z-level !pzlay altitudes of the layers c$$$!zlev altitudes (nlayl+1) of the boundaries c$$$ if (meteo_flux_number .ge. 0) then c$$$ do ig=1, ngrid c$$$ l=1 c$$$ do while ( (((meteo_alt .ge. pplev(ig,l)) .and. c$$$ & (meteo_alt .le. pplev(ig,l+1))) .eq. .false.) c$$$ & .and. (l .lt. nlay) ) c$$$ l=l+1 c$$$ enddo c$$$ meteo_lvl=l c$$$ subpdq(ig,meteo_lvl,igcm_dust_mass)= c$$$ & subpdq(ig,meteo_lvl,igcm_dust_mass) c$$$ & +meteo_flux_mass c$$$ subpdq(ig,meteo_lvl,igcm_dust_number)= c$$$ & subpdq(ig,meteo_lvl,igcm_dust_number) c$$$ & +meteo_flux_number c$$$ enddo c$$$ endif c------------------------------------------------------------------- c 2. Main call to the cloud schemes: c------------------------------------------------ CALL improvedCO2clouds(ngrid,nlay,microtimestep, & pplay,pzlev,zteff,subpdt, & pqeff,subpdq,subpdqcloudco2,subpdtcloudco2, & nq,tauscaling,memdMMccn,memdMMh2o,memdNNccn) c------------------------------------------------------------------- c 3. Updating tendencies after cloud scheme: c----------------------------------------------- DO l=1,nlay DO ig=1,ngrid subpdt(ig,l) = & subpdt(ig,l) + subpdtcloudco2(ig,l) subpdq(ig,l,igcm_dust_mass) = & subpdq(ig,l,igcm_dust_mass) & + subpdqcloudco2(ig,l,igcm_dust_mass) subpdq(ig,l,igcm_dust_number) = & subpdq(ig,l,igcm_dust_number) & + subpdqcloudco2(ig,l,igcm_dust_number) subpdq(ig,l,igcm_ccnco2_mass) = & subpdq(ig,l,igcm_ccnco2_mass) & + subpdqcloudco2(ig,l,igcm_ccnco2_mass) subpdq(ig,l,igcm_ccnco2_number) = & subpdq(ig,l,igcm_ccnco2_number) & + subpdqcloudco2(ig,l,igcm_ccnco2_number) subpdq(ig,l,igcm_co2_ice) = & subpdq(ig,l,igcm_co2_ice) & + subpdqcloudco2(ig,l,igcm_co2_ice) subpdq(ig,l,igcm_co2) = & subpdq(ig,l,igcm_co2) & + subpdqcloudco2(ig,l,igcm_co2) subpdq(ig,l,igcm_h2o_ice) = & subpdq(ig,l,igcm_h2o_ice) & + subpdqcloudco2(ig,l,igcm_h2o_ice) subpdq(ig,l,igcm_ccn_mass) = & subpdq(ig,l,igcm_ccn_mass) & + subpdqcloudco2(ig,l,igcm_ccn_mass) subpdq(ig,l,igcm_ccn_number) = & subpdq(ig,l,igcm_ccn_number) & + subpdqcloudco2(ig,l,igcm_ccn_number) ENDDO ENDDO !Sedimentation !Update traceurs, compute rsed DO l=1, nlay DO ig=1,ngrid tempo_traceur_t(ig,l)=zteff(ig,l)+subpdt(ig,l) & *microtimestep tempo_traceurs(ig,l,:)=pqeff(ig,l,:) & +subpdq(ig,l,:)*microtimestep rho_ice_co2T(ig,l)=1000.*(1.72391-2.53e-4* & tempo_traceur_t(ig,l)-2.87e-6* & tempo_traceur_t(ig,l)*tempo_traceur_t(ig,l)) rho_ice_co2=rho_ice_co2T(ig,l) Niceco2=max(tempo_traceurs(ig,l,igcm_co2_ice),1.e-30) Nccnco2=max(tempo_traceurs(ig,l,igcm_ccnco2_number), & 1.e-30) Qccnco2=max(tempo_traceurs(ig,l,igcm_ccnco2_mass), & 1.e-30) mdustJA= tempo_traceurs(ig,l,igcm_dust_mass) ndustJA=tempo_traceurs(ig,l,igcm_dust_number) if ((Nccnco2 .lt. tauscaling(ig)) .or. (Qccnco2 .lt. & 1.e-30 *tauscaling(ig))) then rdust(ig,l)=1.e-10 else rdust(ig,l)=(3./4./pi/2500.*Qccnco2/Nccnco2)**(1./3.) rdust(ig,l)=max(rdust(ig,l),1.e-10) c rdust(ig,l)=min(rdust(ig,l),5.e-6) end if rhocloudco2t(ig,l) = (Niceco2 *rho_ice_co2 & + Qccnco2*tauscaling(ig)*rho_dust) & / (Niceco2 + Qccnco2*tauscaling(ig)) rhocloudco2t(ig,l)=min(max(rhocloudco2t(ig,l) & ,rho_ice_co2),rho_dust) riceco2(ig,l)=(Niceco2*3.0/ & (4.0*rho_ice_co2*pi*Nccnco2 & *tauscaling(ig)) +rdust(ig,l)*rdust(ig,l) & *rdust(ig,l))**(1.0/3.0) riceco2(ig,l)=max(1.e-10,riceco2(ig,l)) rsedcloudco2(ig,l)=max(riceco2(ig,l)* & (1.+nuiceco2_sed)*(1.+nuiceco2_sed)*(1.+nuiceco2_sed), & rdust(ig,l)) ENDDO ENDDO ! Gravitational sedimentation sav_trac(:,:,igcm_co2_ice)=tempo_traceurs(:,:,igcm_co2_ice) sav_trac(:,:,igcm_ccnco2_mass)= & tempo_traceurs(:,:,igcm_ccnco2_mass) sav_trac(:,:,igcm_ccnco2_number)= & tempo_traceurs(:,:,igcm_ccnco2_number) call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay, & microtimestep,pplev,masse,epaisseur,tempo_traceur_t, & rsedcloudco2,rhocloudco2t, & tempo_traceurs(:,:,igcm_co2_ice),wq,beta) ! 3 traceurs ! sedim at the surface of co2 ice : keep track of it for physiq_mod do ig=1,ngrid pdqs_sedco2(ig)=pdqs_sedco2(ig)+ wq(ig,1)/microtimestep end do call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay, & microtimestep,pplev,masse,epaisseur,tempo_traceur_t, & rsedcloudco2,rhocloudco2t, & tempo_traceurs(:,:,igcm_ccnco2_mass),wq,beta) call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay, & microtimestep,pplev,masse,epaisseur,tempo_traceur_t, & rsedcloudco2,rhocloudco2t, & tempo_traceurs(:,:,igcm_ccnco2_number),wq,beta) DO l = 1, nlay !Compute tendencies DO ig=1,ngrid pdqsed(ig,l,igcm_ccnco2_mass)= & (tempo_traceurs(ig,l,igcm_ccnco2_mass)- & sav_trac(ig,l,igcm_ccnco2_mass))/microtimestep pdqsed(ig,l,igcm_ccnco2_number)= & (tempo_traceurs(ig,l,igcm_ccnco2_number)- & sav_trac(ig,l,igcm_ccnco2_number))/microtimestep pdqsed(ig,l,igcm_co2_ice)= & (tempo_traceurs(ig,l,igcm_co2_ice)- & sav_trac(ig,l,igcm_co2_ice))/microtimestep ENDDO ENDDO !pdqsed est la tendance due a la sedimentation DO l=1,nlay DO ig=1,ngrid subpdq(ig,l,igcm_ccnco2_mass) = & subpdq(ig,l,igcm_ccnco2_mass) & +pdqsed(ig,l,igcm_ccnco2_mass) subpdq(ig,l,igcm_ccnco2_number) = & subpdq(ig,l,igcm_ccnco2_number) & +pdqsed(ig,l,igcm_ccnco2_number) subpdq(ig,l,igcm_co2_ice) = & subpdq(ig,l,igcm_co2_ice) & +pdqsed(ig,l,igcm_co2_ice) ENDDO ENDDO ENDDO ! of DO microstep=1,imicro c------------------------------------------------------------------- c 6. Compute final tendencies after time loop: c------------------------------------------------ c CO2 flux at surface (kg.m-2.s-1) do ig=1,ngrid pdqs_sedco2(ig)=pdqs_sedco2(ig)/real(imicro) enddo c------ Temperature tendency pdtcloud DO l=1,nlay DO ig=1,ngrid pdtcloudco2(ig,l) = & subpdt(ig,l)/real(imicro)-pdt(ig,l) ENDDO ENDDO c------ Tracers tendencies pdqcloud DO l=1,nlay DO ig=1,ngrid pdqcloudco2(ig,l,igcm_co2_ice) = & subpdq(ig,l,igcm_co2_ice)/real(imicro) & - pdq(ig,l,igcm_co2_ice) pdqcloudco2(ig,l,igcm_co2) = & subpdq(ig,l,igcm_co2)/real(imicro) & - pdq(ig,l,igcm_co2) pdqcloudco2(ig,l,igcm_h2o_ice) = & subpdq(ig,l,igcm_h2o_ice)/real(imicro) & - pdq(ig,l,igcm_h2o_ice) ENDDO ENDDO DO l=1,nlay DO ig=1,ngrid pdqcloudco2(ig,l,igcm_ccnco2_mass) = & subpdq(ig,l,igcm_ccnco2_mass)/real(imicro) & - pdq(ig,l,igcm_ccnco2_mass) pdqcloudco2(ig,l,igcm_ccnco2_number) = & subpdq(ig,l,igcm_ccnco2_number)/real(imicro) & - pdq(ig,l,igcm_ccnco2_number) pdqcloudco2(ig,l,igcm_ccn_mass) = & subpdq(ig,l,igcm_ccn_mass)/real(imicro) & - pdq(ig,l,igcm_ccn_mass) pdqcloudco2(ig,l,igcm_ccn_number) = & subpdq(ig,l,igcm_ccn_number)/real(imicro) & - pdq(ig,l,igcm_ccn_number) ENDDO ENDDO DO l=1,nlay DO ig=1,ngrid pdqcloudco2(ig,l,igcm_dust_mass) = & subpdq(ig,l,igcm_dust_mass)/real(imicro) & - pdq(ig,l,igcm_dust_mass) pdqcloudco2(ig,l,igcm_dust_number) = & subpdq(ig,l,igcm_dust_number)/real(imicro) & - pdq(ig,l,igcm_dust_number) ENDDO ENDDO c------- Due to stepped entry, other processes tendencies can add up to negative values c------- Therefore, enforce positive values and conserve mass DO l=1,nlay DO ig=1,ngrid IF ((pqeff(ig,l,igcm_ccnco2_number) + & ptimestep* (pdq(ig,l,igcm_ccnco2_number) + & pdqcloudco2(ig,l,igcm_ccnco2_number)) & .lt. 1.e-20) & .or. (pqeff(ig,l,igcm_ccnco2_mass) + & ptimestep* (pdq(ig,l,igcm_ccnco2_mass) + & pdqcloudco2(ig,l,igcm_ccnco2_mass)) & .lt. 1.e-30)) THEN pdqcloudco2(ig,l,igcm_ccnco2_number) = & - pqeff(ig,l,igcm_ccnco2_number)/ptimestep & - pdq(ig,l,igcm_ccnco2_number)+1.e-20 pdqcloudco2(ig,l,igcm_dust_number) = & -pdqcloudco2(ig,l,igcm_ccnco2_number) pdqcloudco2(ig,l,igcm_ccnco2_mass) = & - pqeff(ig,l,igcm_ccnco2_mass)/ptimestep & - pdq(ig,l,igcm_ccnco2_mass)+1.e-30 pdqcloudco2(ig,l,igcm_dust_mass) = & -pdqcloudco2(ig,l,igcm_ccnco2_mass) ENDIF ENDDO ENDDO DO l=1,nlay DO ig=1,ngrid IF ( (pqeff(ig,l,igcm_dust_number) + & ptimestep* (pdq(ig,l,igcm_dust_number) + & pdqcloudco2(ig,l,igcm_dust_number)) .le. 1.e-30) & .or. (pqeff(ig,l,igcm_dust_mass)+ & ptimestep* (pdq(ig,l,igcm_dust_mass) + & pdqcloudco2(ig,l,igcm_dust_mass)) & .le. 1.e-30)) then pdqcloudco2(ig,l,igcm_dust_number) = & - pqeff(ig,l,igcm_dust_number)/ptimestep & - pdq(ig,l,igcm_dust_number)+1.e-30 pdqcloudco2(ig,l,igcm_ccnco2_number) = & -pdqcloudco2(ig,l,igcm_dust_number) pdqcloudco2(ig,l,igcm_dust_mass) = & - pqeff(ig,l,igcm_dust_mass)/ptimestep & - pdq(ig,l,igcm_dust_mass) +1.e-30 pdqcloudco2(ig,l,igcm_ccnco2_mass) = & -pdqcloudco2(ig,l,igcm_dust_mass) ENDIF ENDDO ENDDO !pq+ptime*(pdq+pdqc)=1 ! pdqc=1-pq/ptime-pdq c$$$ c$$$ c$$$ DO l=1,nlay c$$$ DO ig=1,ngrid c$$$ IF (pq(ig,l,igcm_co2_ice) + ptimestep* c$$$ & (pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice)) c$$$ & .lt. 1.e-30) THEN c$$$ pdqcloudco2(ig,l,igcm_co2_ice) = c$$$ & - pq(ig,l,igcm_co2_ice)/ptimestep - pdq(ig,l,igcm_co2_ice) c$$$ pdqcloudco2(ig,l,igcm_co2) = -pdqcloudco2(ig,l,igcm_co2_ice) c$$$ !write(*,*) "WARNING CO2 ICE in co2cloud.F" c$$$ c$$$ ENDIF c$$$ IF (pq(ig,l,igcm_co2) + ptimestep* c$$$ & (pdq(ig,l,igcm_co2) + pdqcloudco2(ig,l,igcm_co2)) c$$$ & .lt. 0.5) THEN c$$$ pdqcloudco2(ig,l,igcm_co2) = c$$$ & - pdq(ig,l,igcm_co2_ice) !- pdq(ig,l,igcm_co2) c$$$c pdqcloudco2(ig,l,igcm_co2_ice) = -pdqcloudco2(ig,l,igcm_co2) c$$$ pdqcloudco2(ig,l,igcm_co2_ice) = c$$$ & - pq(ig,l,igcm_co2_ice)/ptimestep - pdq(ig,l,igcm_co2_ice) c$$$ ! write(*,*) "WARNING CO2 in co2cloud.F" c$$$ ENDIF c$$$ ENDDO c$$$ ENDDO c$$$ DO l=1, nlay DO ig=1,ngrid c call updaterice_microco2( c & pq(ig,l,igcm_co2_ice) + ! ice mass c & (pdq(ig,l,igcm_co2_ice) + ! ice mass c & pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep, ! ice mass c & pq(ig,l,igcm_ccnco2_mass) + ! ccn mass c & (pdq(ig,l,igcm_ccnco2_mass) + ! ccn mass c & pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep, ! ccn mass c & pq(ig,l,igcm_ccnco2_number) + ! ccn number c & (pdq(ig,l,igcm_ccnco2_number) + ! ccn number c & pdqcloudco2(ig,l,igcm_ccnco2_number))*ptimestep, ! ccn number c & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l)) c write(*,*) "in co2clouds, riceco2(ig,l)= ",riceco2(ig,l) Niceco2=pqeff(ig,l,igcm_co2_ice) + & (pdq(ig,l,igcm_co2_ice) + & pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep Nco2=pqeff(ig,l,igcm_co2) + & (pdq(ig,l,igcm_co2) + & pdqcloudco2(ig,l,igcm_co2))*ptimestep Nccnco2=max((pqeff(ig,l,igcm_ccnco2_number) + & (pdq(ig,l,igcm_ccnco2_number) + & pdqcloudco2(ig,l,igcm_ccnco2_number))*ptimestep)* & tauscaling(ig),1.e-30) Qccnco2=max((pqeff(ig,l,igcm_ccnco2_mass) + & (pdq(ig,l,igcm_ccnco2_mass) + & pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep)* & tauscaling(ig),1.e-30) if (Nccnco2 .lt. 0.1) then rdust(ig,l)=1.e-10 else rdust(ig,l)= Qccnco2 & *0.75/pi/rho_dust & / Nccnco2 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)) endif myT=zteff(ig,l)+(pdt(ig,l)+pdtcloudco2(ig,l))*ptimestep rho_ice_co2T(ig,l)=1000.*(1.72391-2.53e-4* & myT-2.87e-6* myT* myT) rho_ice_co2=rho_ice_co2T(ig,l) lw = l0 + l1 * myT + l2 *myT**2 + & l3 * myT**3 + l4 * myT**4 !J.kg-1 riceco2(ig,l)=(Niceco2*3.0/ & (4.0*rho_ice_co2*pi*Nccnco2) & +rdust(ig,l)*rdust(ig,l) & *rdust(ig,l) )**(1.0/3.0) c & .or. (riceco2(ig,l) .le. rdust(ig,l)) if ( (Niceco2 .le. 1.e-25).or. (Nccnco2 .le. 0.1) )THEN !anciennement 200 c riceco2(ig,l)=0. c & .or. (Nccnco2 .le. 1.) c endif !Flux chaleur latente <0 quand sublimation pdtcloudco2(ig,l)= pdtcloudco2(ig,l)-Niceco2*lw/cpp/ptimestep c$$$ !NO CLOUD : RESET TRACERS AND CONSERVE MASS c if (pq(ig,l,igcm_co2_ice)+(pdq(ig,l,igcm_co2_ice)+ c & pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep .le. 0.) then c pdqcloudco2(ig,l,igcm_co2)=0. c pdqcloudco2(ig,l,igcm_co2_ice)=0. c else pdqcloudco2(ig,l,igcm_co2_ice)=-1.*pqeff(ig,l,igcm_co2_ice) & /ptimestep-pdq(ig,l,igcm_co2_ice) pdqcloudco2(ig,l,igcm_co2)=-1.* & pdqcloudco2(ig,l,igcm_co2_ice) c endif ! Reverse h2o ccn and ices into h2o tracers if (memdMMccn(ig,l) .gt. 0) then pdqcloudco2(ig,l,igcm_ccn_mass)=memdMMccn(ig,l)/ptimestep else memdMMccn(ig,l)=0. pdqcloudco2(ig,l,igcm_ccn_mass)=0. endif if (memdNNccn(ig,l) .gt. 0) then pdqcloudco2(ig,l,igcm_ccn_number)=memdNNccn(ig,l)/ptimestep else memdNNccn(ig,l)=0. pdqcloudco2(ig,l,igcm_ccn_number)=0. endif if (memdMMh2o(ig,l) .gt. 0) then pdqcloudco2(ig,l,igcm_h2o_ice)=memdMMh2o(ig,l)/ptimestep else memdMMh2o(ig,l)=0. pdqcloudco2(ig,l,igcm_h2o_ice)=0. endif c if (pq(ig,l,igcm_ccnco2_mass)+(pdq(ig,l,igcm_ccnco2_mass)+ c & pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep c & .le. 1.e-30) then c pdqcloudco2(ig,l,igcm_ccnco2_mass)=0. c pdqcloudco2(ig,l,igcm_ccnco2_number)=0. c pdqcloudco2(ig,l,igcm_co2)=0. c pdqcloudco2(ig,l,igcm_co2_ice)=0. c else pdqcloudco2(ig,l,igcm_ccnco2_mass)= & -1.*pqeff(ig,l,igcm_ccnco2_mass) & /ptimestep-pdq(ig,l,igcm_ccnco2_mass) c endif c if (pq(ig,l,igcm_ccnco2_number)+ c & (pdq(ig,l,igcm_ccnco2_number)+ c & pdqcloudco2(ig,l,igcm_ccnco2_number)) c & *ptimestep.le. 1.e-30) c & then c pdqcloudco2(ig,l,igcm_ccnco2_mass)=0. c pdqcloudco2(ig,l,igcm_ccnco2_number)=0. c pdqcloudco2(ig,l,igcm_co2)=0. c pdqcloudco2(ig,l,igcm_co2_ice)=0. c else pdqcloudco2(ig,l,igcm_ccnco2_number)= & -1.*pqeff(ig,l,igcm_ccnco2_number) & /ptimestep-pdq(ig,l,igcm_ccnco2_number) c endif c if (pq(ig,l,igcm_dust_number)+ c & (pdq(ig,l,igcm_dust_number)+ c & pdqcloudco2(ig,l,igcm_dust_number)) c & *ptimestep.le. 1.e-30) c & then c pdqcloudco2(ig,l,igcm_dust_number)=0. c pdqcloudco2(ig,l,igcm_dust_mass)=0. c else pdqcloudco2(ig,l,igcm_dust_number)= & pqeff(ig,l,igcm_ccnco2_number) & /ptimestep+pdq(ig,l,igcm_ccnco2_number) & -memdNNccn(ig,l)/ptimestep c endif c if (pq(ig,l,igcm_dust_mass)+ c & (pdq(ig,l,igcm_dust_mass)+ c & pdqcloudco2(ig,l,igcm_dust_mass)) c & *ptimestep .le. 1.e-30) c & then c pdqcloudco2(ig,l,igcm_dust_number)=0. c pdqcloudco2(ig,l,igcm_dust_mass)=0. c else pdqcloudco2(ig,l,igcm_dust_mass)= & pqeff(ig,l,igcm_ccnco2_mass) & /ptimestep+pdq(ig,l,igcm_ccnco2_mass) & -(memdMMh2o(ig,l)+memdMMccn(ig,l))/ptimestep c endif memdMMccn(ig,l)=0. memdMMh2o(ig,l)=0. memdNNccn(ig,l)=0. riceco2(ig,l)=0. endif c Compute opacities No=Nccnco2 Rn=-log(riceco2(ig,l)) n_derf = erf( (rb_cldco2(1)+Rn) *dev2) Qext1bins2(ig,l)=0. do i = 1, nbinco2_cld !Qext below 50 is negligible n_aer(i) = -0.5 * No * n_derf !! this ith previously computed n_derf = derf((rb_cldco2(i+1)+Rn) *dev2) n_aer(i) = n_aer(i) + 0.5 * No * n_derf Qext1bins2(ig,l)=Qext1bins2(ig,l)+Qext1bins(i)*n_aer(i) enddo !l'opacité de la case ig est la somme sur l de Qext1bins2 !update rice water call updaterice_micro( & pqeff(ig,l,igcm_h2o_ice) + ! ice mass & (pdq(ig,l,igcm_h2o_ice) + ! ice mass & pdqcloudco2(ig,l,igcm_h2o_ice))*ptimestep, ! ice mass & pqeff(ig,l,igcm_ccn_mass) + ! ccn mass & (pdq(ig,l,igcm_ccn_mass) + ! ccn mass & pdqcloudco2(ig,l,igcm_ccn_mass))*ptimestep, ! ccn mass & pqeff(ig,l,igcm_ccn_number) + ! ccn number & (pdq(ig,l,igcm_ccn_number) + ! ccn number & pdqcloudco2(ig,l,igcm_ccn_number))*ptimestep, ! ccn number & tauscaling(ig),rice(ig,l),rhocloud(ig,l)) call updaterdust( & pqeff(ig,l,igcm_dust_mass) + ! dust mass & (pdq(ig,l,igcm_dust_mass) + ! dust mass & pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep, ! dust mass & pqeff(ig,l,igcm_dust_number) + ! dust number & (pdq(ig,l,igcm_dust_number) + ! dust number & pdqcloudco2(ig,l,igcm_dust_number))*ptimestep, ! dust number & rdust(ig,l)) ENDDO ENDDO tau1mic(:)=0. do l=1,nlay do ig=1,ngrid tau1mic(ig)=tau1mic(ig)+Qext1bins2(ig,l) enddo enddo c$$$ c$$$ if (riceco2(725,22) .ge. 1.e-10) then c$$$ c$$$ write(*,*) 'DIAGJA co2 ',pqeff(725,22,igcm_co2) + c$$$ & (pdq(725,22,igcm_co2) + c$$$ & pdqcloudco2(725,22,igcm_co2))*ptimestep c$$$ write(*,*) 'DIAGJA co2_ice',pqeff(725,22,igcm_co2_ice) + c$$$ & (pdq(725,22,igcm_co2_ice) + c$$$ & pdqcloudco2(725,22,igcm_co2_ice))*ptimestep c$$$ c$$$ write(*,*) 'DIAGJA riceco2',riceco2(725,22) c$$$ write(*,*) 'DIAGJA T',zteff(725,22) + c$$$ & (pdt(725,22) + pdtcloudco2(725,22))*ptimestep c$$$ write(*,*) 'DIAG pdtcloud',pdtcloudco2(725,22)*ptimestep c$$$ write(*,*) 'DIAGJA ccnNco2',pqeff(725,22,igcm_ccnco2_number)+ c$$$ & (pdq(725,22,igcm_ccnco2_number) + c$$$ & pdqcloudco2(725,22,igcm_ccnco2_number))*ptimestep c$$$ c$$$ write(*,*) 'DIAGJA dustN',pqeff(725,22,igcm_dust_number) + c$$$ & (pdq(725,22,igcm_dust_number) + c$$$ & pdqcloudco2(725,22,igcm_dust_number))*ptimestep c$$$ ENDIF c$$$ c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c Then that should not affect the ice particle radius do ig=1,ngrid if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3))) & riceco2(ig,2)=riceco2(ig,3) riceco2(ig,1)=riceco2(ig,2) end if end do DO l=1,nlay DO ig=1,ngrid rsedcloud(ig,l)=max(rice(ig,l)* & (1.+nuice_sed)*(1.+nuice_sed)*(1.+nuice_sed), & rdust(ig,l)) ! rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4) ENDDO ENDDO DO l=1,nlay DO ig=1,ngrid rsedcloudco2(ig,l)=max(riceco2(ig,l)* & (1.+nuiceco2_sed)*(1.+nuiceco2_sed)*(1.+nuiceco2_sed), & rdust(ig,l)) c rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),1.e-5) ENDDO ENDDO call co2sat(ngrid*nlay,zteff+(pdt+pdtcloudco2)*ptimestep & ,pplay,zqsatco2) do l=1,nlay do ig=1,ngrid satuco2(ig,l) = (pqeff(ig,l,igcm_co2) + & (pdq(ig,l,igcm_co2) + & pdqcloudco2(ig,l,igcm_co2))*ptimestep)* & (mmean(ig,l)/44.01)*pplay(ig,l)/zqsatco2(ig,l) !write(*,*) "In CO2 pt,sat ",pt(ig,l),satuco2(ig,l) enddo enddo !Tout ce qui est modifié par la microphysique de CO2 doit être rapporté a cloudfrac IF (CLFvaryingCO2) THEN DO l=1,nlay DO ig=1,ngrid pdqcloudco2(ig,l,igcm_ccn_mass)= & pdqcloudco2(ig,l,igcm_ccn_mass)*cloudfrac(ig,l) pdqcloudco2(ig,l,igcm_ccnco2_mass)= & pdqcloudco2(ig,l,igcm_ccnco2_mass)*cloudfrac(ig,l) pdqcloudco2(ig,l,igcm_ccn_number)= & pdqcloudco2(ig,l,igcm_ccn_number)*cloudfrac(ig,l) pdqcloudco2(ig,l,igcm_ccnco2_number)= & pdqcloudco2(ig,l,igcm_ccnco2_number)*cloudfrac(ig,l) pdqcloudco2(ig,l,igcm_dust_mass)= & pdqcloudco2(ig,l,igcm_dust_mass)*cloudfrac(ig,l) pdqcloudco2(ig,l,igcm_dust_number)= & pdqcloudco2(ig,l,igcm_dust_number)*cloudfrac(ig,l) pdqcloudco2(ig,l,igcm_h2o_ice)= & pdqcloudco2(ig,l,igcm_h2o_ice)*cloudfrac(ig,l) pdqcloudco2(ig,l,igcm_co2_ice)= & pdqcloudco2(ig,l,igcm_co2_ice)*cloudfrac(ig,l) pdqcloudco2(ig,l,igcm_co2)= & pdqcloudco2(ig,l,igcm_co2)*cloudfrac(ig,l) pdtcloudco2(ig,l)=pdtcloudco2(ig,l)*cloudfrac(ig,l) ENDDO ENDDO ENDIF call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",3, & satuco2) call WRITEdiagfi(ngrid,"riceco2","ice radius","m" & ,3,riceco2) ! or output in diagfi.nc (for testphys1d) c call WRITEDIAGFI(ngrid,'ps','Surface pressure','Pa',0,ps) c call WRITEDIAGFI(ngrid,'temp','Temperature ', c & 'K JA',1,pt) call WRITEdiagfi(ngrid,"rsedcloudco2","rsed co2","m",3, & rsedcloudco2) call WRITEdiagfi(ngrid,"tau1mic","co2 ice opacity 1 micron"," ",2, & tau1mic) call WRITEdiagfi(ngrid,"cloudfrac","co2 cloud fraction"," ",3, & cloudfrac) ! used for rad. transfer calculations ! nuice is constant because a lognormal distribution is prescribed c nuice(1:ngrid,1:nlay)=nuice_ref c======================================================================= END