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,zlev,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, & meteo_flux_mass,meteo_flux_number, & meteo_alt IMPLICIT NONE 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 c======================================================================= c----------------------------------------------------------------------- c declarations: c ------------- !#include "dimensions.h" !#include "dimphys.h" #include "callkeys.h" !#include "tracer.h" !#include "comgeomfi.h" !#include "dimradmars.h" ! naerkind is set in scatterers.h (built when compiling with makegcm -s #) !#include"scatterers.h" #include "microphys.h" 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) :: zlev(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 Mo,No 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 LOGICAL,SAVE :: firstcall=.true. DOUBLE PRECISION Nccnco2, Niceco2,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 integer meteo_lvl real tempo_traceur_t(ngrid,nlay) real tempo_traceurs(ngrid,nlay,nq) real sav_trac(ngrid,nlay,nq) real pdqsed(ngrid,nlay,nq) DOUBLE PRECISION,allocatable,save :: memdMMccn(:,:) DOUBLE PRECISION,allocatable,save :: memdMMh2o(:,:) DOUBLE PRECISION,allocatable,save :: memdNNccn(:,:) 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 if (meteo_flux_number .ne. 0) then write(*,*) "Warning ! Meteoritic flux of dust is turned on" write(*,*) "Dust mass = ",meteo_flux_mass write(*,*) "Dust number = ",meteo_flux_number write(*,*) "Are added at the z-level (km)",meteo_alt write(*,*) "Every timestep in co2cloud.F" 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. firstcall=.false. ENDIF ! of IF (firstcall) c-----Initialization 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 add meteoritic flux of dust !convert meteo_alt (in km) to z-level !pzlay altitudes of the layers if (meteo_flux_number .ne. 0) then do ig=1, ngrid diff0=1000 meteo_lvl=1 do l=1,nlay diff=pzlay(ig,l)/1000-meteo_alt if (abs(diff) .lt. diff0) then diff0=abs(diff) meteo_lvl=l endif enddo c write(*,*) "meteo_lvl=",meteo_lvl pdq(ig,meteo_lvl,igcm_dust_mass)= & pdq(ig,meteo_lvl,igcm_dust_mass) & +dble(meteo_flux_mass)/tauscaling(ig) pdq(ig,meteo_lvl,igcm_dust_number)= & pdq(ig,meteo_lvl,igcm_dust_number) & +dble(meteo_flux_number)/tauscaling(ig) enddo endif c call WRITEDIAGFI(ngrid,"pzlay","altitude","km",3, c & pzlay) 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 tempo_traceur_t(1:ngrid,1:nlay)=0 tempo_traceurs(1:ngrid,1:nlay,1:nq)=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)= zlev(ig,l+1) - zlev(ig,l) enddo enddo c------------------------------------------------------------------- c 1. Tendencies: c------------------ 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 c tempo_traceur_t(ig,l)=tempo_traceur_t(ig,l) c & + subpdtcloudco2(ig,l) !write(*,*) 'T micro= ', tempo_traceur_t(ig,l) c tempo_traceurs(ig,l,:)=tempo_traceurs(ig,l,:) c & +subpdqcloudco2(ig,l,:) 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) tempo_traceur_t(ig,l)= pt(ig,l)+subpdt(ig,l)*microtimestep tempo_traceurs(ig,l,:)= pq(ig,l,:)+subpdq(ig,l,:) & *microtimestep !Stepped entry for sedimentation ENDDO ENDDO !RSEDCLOUD AND RICECO2 HERE DO l=1, nlay DO ig=1,ngrid 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) ! 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)) c riceco2(ig,l)= Niceco2*3.0/ c & (4.0*rho_ice_co2*pi*Nccnco2) c & +rdust(ig,l)*rdust(ig,l)*rdust(ig,l) c riceco2(ig,l)=riceco2(ig,l)**(1.0/3.0) c write(*,*) "in co2clouds, rice = ",riceco2(ig,l) c write(*,*) "in co2clouds, rho = ",rhocloudco2t(ig,l) 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)) !riceco2(ig,l)=min(5.e-5,riceco2(ig,l)) c call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2, c & tauscaling(ig),riceco2(ig,l),rhocloudco2t(ig,l)) c write(*,*) "in co2clouds, rice update = ",riceco2(ig,l) c write(*,*) "in co2clouds, rho update = " c & ,rhocloudco2t(ig,l) rsedcloudco2(ig,l)=max(riceco2(ig,l)* & (1.+nuiceco2_sed)*(1.+nuiceco2_sed)*(1.+nuiceco2_sed), & rdust(ig,l)) rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),1.e-5) c write(*,*) 'Rsedcloud = ',rsedcloudco2(ig,l) !write(*,*) 'Rhocloudco2 = ',rhocloudco2t(ig,l) ENDDO ENDDO ! Gravitational sedimentation ! sedimentation computed from radius computed from q in module radii_mod 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 do ig=1,ngrid pdqs_sedco2(ig)=pdqs_sedco2(ig)+ wq(ig,1) 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 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 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 c------------------------------------------------------------------- c 2. Main call to the different cloud schemes: c------------------------------------------------ IF (microphysco2) THEN CALL improvedCO2clouds(ngrid,nlay,microtimestep, & pplay,pt,subpdt, & pq,subpdq,subpdqcloudco2,subpdtcloudco2, & nq,tauscaling,memdMMccn,memdMMh2o,memdNNccn) ELSE write(*,*) ' no simpleCO2clouds procedure: STOP' ! listo write(*,*) ' microphysco2 and co2clouds must be .true.' ! listo STOP ENDIF c------------------------------------------------------------------- c 3. Updating tendencies after cloud scheme: c----------------------------------------------- c IF (microphysco2) THEN DO l=1,nlay DO ig=1,ngrid 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) c & +pdqsed(ig,l,igcm_ccnco2_mass) subpdq(ig,l,igcm_ccnco2_number) = & subpdq(ig,l,igcm_ccnco2_number) & + subpdqcloudco2(ig,l,igcm_ccnco2_number) c & +pdqsed(ig,l,igcm_ccnco2_number) subpdq(ig,l,igcm_co2_ice) = & subpdq(ig,l,igcm_co2_ice) & + subpdqcloudco2(ig,l,igcm_co2_ice) c & +pdqsed(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 IF (activice) THEN DO l=1,nlay DO ig=1,ngrid subpdt(ig,l) = & subpdt(ig,l) + subpdtcloudco2(ig,l) ENDDO ENDDO ENDIF 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)/ptimestep 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 ! call WRITEdiagfi(ngrid,"co2cloud00","co2 traceur","kg/kg",1, ! & pq(1,:,igcm_co2_ice) + ptimestep* ! & (pdq(1,:,igcm_co2_ice) + pdqcloudco2(1,:,igcm_co2_ice))) c IF(microphysco2) THEN 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 c ENDIF IF(scavenging) THEN 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 ENDIF c ENDIF c------- Due to stepped entry, other processes tendencies can add up to negative values c------- Therefore, enforce positive values and conserve mass IF(microphysco2) THEN DO l=1,nlay DO ig=1,ngrid IF ((pq(ig,l,igcm_ccnco2_number) + & ptimestep* (pdq(ig,l,igcm_ccnco2_number) + & pdqcloudco2(ig,l,igcm_ccnco2_number)) & .lt. 1.e-15) & .or. (pq(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) = & - pq(ig,l,igcm_ccnco2_number)/ptimestep & - pdq(ig,l,igcm_ccnco2_number)+1.e-15 pdqcloudco2(ig,l,igcm_dust_number) = & -pdqcloudco2(ig,l,igcm_ccnco2_number) pdqcloudco2(ig,l,igcm_ccnco2_mass) = & - pq(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 ENDIF IF(scavenging) THEN DO l=1,nlay DO ig=1,ngrid IF ( (pq(ig,l,igcm_dust_number) + & ptimestep* (pdq(ig,l,igcm_dust_number) + & pdqcloudco2(ig,l,igcm_dust_number)) .le. 1.e-15) & .or. (pq(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) = & - pq(ig,l,igcm_dust_number)/ptimestep & - pdq(ig,l,igcm_dust_number)+1.e-15 pdqcloudco2(ig,l,igcm_ccnco2_number) = & -pdqcloudco2(ig,l,igcm_dust_number) pdqcloudco2(ig,l,igcm_dust_mass) = & - pq(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 ENDIF !pq+ptime*(pdq+pdqc)=1 ! pdqc=1-pq/ptime-pdq 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-25) 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 ENDIF c IF (pq(ig,l,igcm_co2) + ptimestep* c & (pdq(ig,l,igcm_co2) + pdqcloudco2(ig,l,igcm_co2)) c & .lt. 1.e-10) THEN c pdqcloudco2(ig,l,igcm_co2) = c & - pq(ig,l,igcm_co2)/ptimestep - pdq(ig,l,igcm_co2) c pdqcloudco2(ig,l,igcm_co2_ice) = -pdqcloudco2(ig,l,igcm_co2) c write(*,*) "WARNING CO2 in co2cloud.F" c ENDIF c ENDDO c ENDDO c------Update the ice and dust particle size "riceco2" for output or photochemistry c------Only rsedcloudco2 is used for the co2 (cloud) cycle c IF(scavenging) THEN c DO l=1, nlay c DO ig=1,ngrid c call updaterdust( c & pq(ig,l,igcm_dust_mass) + ! dust mass c & (pdq(ig,l,igcm_dust_mass) + ! dust mass c & pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep, ! dust mass c & pq(ig,l,igcm_dust_number) + ! dust number c & (pdq(ig,l,igcm_dust_number) + ! dust number c & pdqcloudco2(ig,l,igcm_dust_number))*ptimestep, ! dust number c & rdust(ig,l)) c write(*,*) "in co2clouds, rdust(ig,l)= ",rdust(ig,l) c mdustJA= pq(ig,l,igcm_dust_mass) + c & (pdq(ig,l,igcm_dust_mass) + c & pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep c ndustJA=pq(ig,l,igcm_dust_number) + c & (pdq(ig,l,igcm_dust_number) + c & pdqcloudco2(ig,l,igcm_dust_number))*ptimestep c if ((ndustJA .lt. tauscaling(ig)) .or. (mdustJA .lt. c & 1.e-30 *tauscaling(ig))) then c rdust(ig,l)=1.e-10 c else c rdust(ig,l)=(3./4./pi/2500.*mdustJA/ndustJA)**(1./3.) c rdust(ig,l)=min(rdust(ig,l),5.e-6) c rdust(ig,l)=max(rdust(ig,l),1.e-10) c endif c ENDDO c ENDDO c ENDIF IF(microphysco2) THEN 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=pq(ig,l,igcm_co2_ice) + & (pdq(ig,l,igcm_co2_ice) + & pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep Nccnco2=max((pq(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((pq(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. 1) .or. (Qccnco2 .lt. & 1.e-30)) 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)) !rdust(ig,l)=min(5.e-5,rdust(ig,l)) endif rho_ice_co2T(ig,l)=1000.*(1.72391-2.53e-4* & (pt(ig,l)+(subpdt(ig,l)+pdtcloudco2(ig,l))*ptimestep) & -2.87e-6* & (pt(ig,l)+(subpdt(ig,l)+pdtcloudco2(ig,l))*ptimestep)* & (pt(ig,l)+(subpdt(ig,l)+pdtcloudco2(ig,l))*ptimestep)) rho_ice_co2=rho_ice_co2T(ig,l) 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) if ( (Niceco2 .le. 1.e-23) .or. & (riceco2(ig,l) .le. rdust(ig,l)) & .or. (Nccnco2 .le. 0.01))then c write(*,*) "Riceco2 co2cloud before reset=",riceco2(ig,l) c write(*,*) "Niceco2 co2cloud before reset=",Niceco2 c write(*,*) "Nccnco2 co2cloud before reset=",Nccnco2 c write(*,*) "Rdust co2cloud before reset=",rdust(ig,l) c$$$ !NO CLOUD : RESET TRACERS AND CONSERVE MASS if (pq(ig,l,igcm_co2_ice)+(pdq(ig,l,igcm_co2_ice)+ & pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep .le. 0.) then pdqcloudco2(ig,l,igcm_co2)=0. pdqcloudco2(ig,l,igcm_co2_ice)=0. else pdqcloudco2(ig,l,igcm_co2)= pq(ig,l,igcm_co2_ice) & /ptimestep+pdq(ig,l,igcm_co2_ice) pdqcloudco2(ig,l,igcm_co2_ice)=-pq(ig,l,igcm_co2_ice) & /ptimestep-pdq(ig,l,igcm_co2_ice) 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) 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) 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) else memdMMh2o(ig,l)=0. pdqcloudco2(ig,l,igcm_h2o_ice)=0. endif if (pq(ig,l,igcm_ccnco2_mass)+(pdq(ig,l,igcm_ccnco2_mass)+ & pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep & .le. 1.e-30) then pdqcloudco2(ig,l,igcm_ccnco2_mass)=0. pdqcloudco2(ig,l,igcm_ccnco2_number)=0. pdqcloudco2(ig,l,igcm_co2)=0. pdqcloudco2(ig,l,igcm_co2_ice)=0. else pdqcloudco2(ig,l,igcm_ccnco2_mass)= & -pq(ig,l,igcm_ccnco2_mass) & /ptimestep-pdq(ig,l,igcm_ccnco2_mass) endif if (pq(ig,l,igcm_ccnco2_number)+ & (pdq(ig,l,igcm_ccnco2_number)+ & pdqcloudco2(ig,l,igcm_ccnco2_number)) & *ptimestep.le. 1.e-30) & then pdqcloudco2(ig,l,igcm_ccnco2_mass)=0. pdqcloudco2(ig,l,igcm_ccnco2_number)=0. pdqcloudco2(ig,l,igcm_co2)=0. pdqcloudco2(ig,l,igcm_co2_ice)=0. else pdqcloudco2(ig,l,igcm_ccnco2_number)= & -pq(ig,l,igcm_ccnco2_number) & /ptimestep-pdq(ig,l,igcm_ccnco2_number) endif if (pq(ig,l,igcm_dust_number)+ & (pdq(ig,l,igcm_dust_number)+ & pdqcloudco2(ig,l,igcm_dust_number)) & *ptimestep.le. 1.e-30) & then pdqcloudco2(ig,l,igcm_dust_number)=0. pdqcloudco2(ig,l,igcm_dust_mass)=0. else pdqcloudco2(ig,l,igcm_dust_number)= & pq(ig,l,igcm_ccnco2_number) & /ptimestep+pdq(ig,l,igcm_ccnco2_number) & -memdNNccn(ig,l) endif if (pq(ig,l,igcm_dust_mass)+ & (pdq(ig,l,igcm_dust_mass)+ & pdqcloudco2(ig,l,igcm_dust_mass)) & *ptimestep .le. 1.e-30) & then pdqcloudco2(ig,l,igcm_dust_number)=0. pdqcloudco2(ig,l,igcm_dust_mass)=0. else pdqcloudco2(ig,l,igcm_dust_mass)= & pq(ig,l,igcm_ccnco2_mass) & /ptimestep+pdq(ig,l,igcm_ccnco2_mass) & -(memdMMh2o(ig,l)+memdMMccn(ig,l)) endif memdMMccn(ig,l)=0. memdMMh2o(ig,l)=0. memdNNccn(ig,l)=0. riceco2(ig,l)=0. pdtcloudco2(ig,l)=0. endif !update rice water call updaterice_micro( & pq(ig,l,igcm_h2o_ice) + ! ice mass & (pdq(ig,l,igcm_h2o_ice) + ! ice mass & pdqcloudco2(ig,l,igcm_h2o_ice))*ptimestep, ! ice mass & pq(ig,l,igcm_ccn_mass) + ! ccn mass & (pdq(ig,l,igcm_ccn_mass) + ! ccn mass & pdqcloudco2(ig,l,igcm_ccn_mass))*ptimestep, ! ccn mass & pq(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( & pq(ig,l,igcm_dust_mass) + ! dust mass & (pdq(ig,l,igcm_dust_mass) + ! dust mass & pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep, ! dust mass & pq(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 ELSE ! no microphys ! not of concern for co2 clouds - listo ENDIF ! of IF(microphysco2) c TO CHEK for relevancy - listo 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 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)) rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),1.e-5) ENDDO ENDDO call co2sat(ngrid*nlay,pt,pplay,zqsatco2) do ig=1,ngrid do l=1,nlay satuco2(ig,l) = (pq(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) c write(*,*) "In CO2 pt,sat ",pt(ig,l),satuco2(ig,l) enddo enddo 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) ! used for rad. transfer calculations ! nuice is constant because a lognormal distribution is prescribed c nuice(1:ngrid,1:nlay)=nuice_ref c======================================================================= END