SUBROUTINE co2cloud(ngrid,nlay,ptimestep, & pplev,pplay,pdpsrf,pzlay,pt,pdt, & pq,pdq,pdqcloudco2,pdtcloudco2, & nq,tau,tauscaling,rdust,rice,riceco2,nuice, & rsedcloudco2,rhocloudco2,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_ccnco2_mass, igcm_ccnco2_number, & rho_dust, nuiceco2_sed, nuiceco2_ref, & rho_ice_co2,r3n_q, & 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 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 ------ ! 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) 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: 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 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" firstcall=.false. ENDIF ! of IF (firstcall) c-----Initialization c add meteoritic flux of dust !convert meteo_alt (in km) to z-level !pzlay altitudes of the layers 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 pq(ig,meteo_lvl,igcm_dust_mass)=pq(ig,meteo_lvl,igcm_dust_mass) & +dble(meteo_flux_mass) pq(ig,meteo_lvl,igcm_dust_number)= & pq(ig,meteo_lvl,igcm_dust_number) & +dble(meteo_flux_number) enddo call WRITEDIAGFI(ngrid,"pzlay","altitude","km",3, & pzlay) 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 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) 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 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 ((ndustJA .lt. tauscaling(ig)) .or. (mdustJA .lt. & 1.e-30 *tauscaling(ig))) then rdust(ig,l)=1.e-10 else rdust(ig,l)=(3./4./pi/2500.*mdustJA/ndustJA)**(1./3.) rdust(ig,l)=min(rdust(ig,l),5.e-6) rdust(ig,l)=max(rdust(ig,l),1.e-9) end if c rhocloudco2t(ig,l) = (Niceco2 *rho_ice_co2 c & + Qccnco2*rho_dust) c & / (Niceco2 + Qccnco2) riceco2(ig,l)= Niceco2*3.0/ & (4.0*rho_ice_co2*pi*Nccnco2) & +rdust(ig,l)*rdust(ig,l)*rdust(ig,l) 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) call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2, & 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),5.e-4) 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) ELSE write(*,*) ' no simpleCO2clouds procedure: STOP' ! listo STOP c CALL simpleclouds(ngrid,nlay,microtimestep, ! for water-ice clouds c & pplay,pzlay,pt,subpdt, c & pq,subpdq,subpdqcloud,subpdtcloud, c & nq,tau,riceco2) 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) ENDDO ENDDO !ici ! call WRITEdiagfi(ngrid,"co2cloud000","co2 traceur","kg/kg",1, ! & pq(1,:,igcm_co2_ice) + ptimestep* ! & ( subpdq(1,:,igcm_co2_ice))) 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)/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)/imicro & - pdq(ig,l,igcm_co2_ice) pdqcloudco2(ig,l,igcm_co2) = & subpdq(ig,l,igcm_co2)/imicro & - pdq(ig,l,igcm_co2) 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))) IF(microphysco2) THEN DO l=1,nlay DO ig=1,ngrid pdqcloudco2(ig,l,igcm_ccnco2_mass) = & subpdq(ig,l,igcm_ccnco2_mass)/imicro & - pdq(ig,l,igcm_ccnco2_mass) pdqcloudco2(ig,l,igcm_ccnco2_number) = & subpdq(ig,l,igcm_ccnco2_number)/imicro & - pdq(ig,l,igcm_ccnco2_number) ENDDO ENDDO 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. 0) & .or. (pq(ig,l,igcm_ccnco2_mass) + & ptimestep* (pdq(ig,l,igcm_ccnco2_mass) + & pdqcloudco2(ig,l,igcm_ccnco2_mass)) & .lt. 0)) THEN pdqcloudco2(ig,l,igcm_ccnco2_number) = & - pq(ig,l,igcm_ccnco2_number)/ptimestep & - pdq(ig,l,igcm_ccnco2_number) +0 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)+0 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)) .lt. 0.) & .or. (pq(ig,l,igcm_dust_mass)+ & ptimestep* (pdq(ig,l,igcm_dust_mass) + & pdqcloudco2(ig,l,igcm_dust_mass)) & .lt. 0.)) then pdqcloudco2(ig,l,igcm_dust_number) = & - pq(ig,l,igcm_dust_number)/ptimestep & - pdq(ig,l,igcm_dust_number)+0 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) +0 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 DO l=1,nlay DO ig=1,ngrid IF (pq(ig,l,igcm_co2_ice) + ptimestep* & (pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice)) & .lt. 1.e-25) THEN pdqcloudco2(ig,l,igcm_co2_ice) = & - pq(ig,l,igcm_co2_ice)/ptimestep - pdq(ig,l,igcm_co2_ice) & +1.e-25 pdqcloudco2(ig,l,igcm_co2) = -pdqcloudco2(ig,l,igcm_co2_ice) ENDIF ENDDO ENDDO c------Update the ice and dust particle size "riceco2" for output or photochemistry c------Only rsedcloudco2 is used for the co2 (cloud) cycle IF(scavenging) THEN DO l=1, nlay 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) mdustJA= pq(ig,l,igcm_dust_mass) + & (pdq(ig,l,igcm_dust_mass) + & pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep ndustJA=pq(ig,l,igcm_dust_number) + & (pdq(ig,l,igcm_dust_number) + & pdqcloudco2(ig,l,igcm_dust_number))*ptimestep if ((ndustJA .lt. tauscaling(ig)) .or. (mdustJA .lt. & 1.e-30 *tauscaling(ig))) then rdust(ig,l)=1.e-10 else rdust(ig,l)=(3./4./pi/2500.*mdustJA/ndustJA)**(1./3.) rdust(ig,l)=min(rdust(ig,l),5.e-4) rdust(ig,l)=max(rdust(ig,l),1.e-10) endif ENDDO ENDDO 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) c rhocloudco2t(ig,l) = (Niceco2 *rho_ice_co2 + Qccnco2*rho_dust) c & / (Niceco2 + Qccnco2) c rhocloudco2(ig,l) = min(max(rhocloudco2t,rho_ice_co2),rho_dust) c write(*,*) "test, nccnco2 =",nccnco22 riceco2(ig,l)= Niceco2*3.0/ & (4.0*rho_ice_co2*pi*Nccnco2) & +rdust(ig,l)*rdust(ig,l)*rdust(ig,l) riceco2(ig,l)=riceco2(ig,l)**(1.0/3.0) c write(*,*) "In co2cloud, after loop, riceco2 =",riceco2(ig,l) c write(*,*) "In co2cloud, after loop, rhoco2 =" c & ,rhocloudco2t(ig,l) call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2, & tauscaling(ig),riceco2(ig,l),rhocloudco2t(ig,l)) c write(*,*) "In co2cloud, after loop and update, riceco2 =" c & ,riceco2(ig,l) c write(*,*) "In co2cloud, after loop and update, rhoco2 =" c & ,rhocloudco2t(ig,l) if ( Niceco2 & .le. 1.e-23 .or. riceco2(ig,l) .le. 1.e-10 .or. & riceco2(ig,l) .ge. 4.999e-4) then ! .or. riceco2(ig,l) .gt. 1.e-4 ) then riceco2(ig,l)=0. !NO CLOUD : RESET TRACER AND CONSERVE MASS 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) pdqcloudco2(ig,l,igcm_ccnco2_mass)= & -pq(ig,l,igcm_ccnco2_mass) & /ptimestep-pdq(ig,l,igcm_ccnco2_mass) pdqcloudco2(ig,l,igcm_ccnco2_number)= & -pq(ig,l,igcm_ccnco2_number) & /ptimestep-pdq(ig,l,igcm_ccnco2_number) pdqcloudco2(ig,l,igcm_dust_number)= & pq(ig,l,igcm_ccnco2_number) & /ptimestep+pdq(ig,l,igcm_ccnco2_number) pdqcloudco2(ig,l,igcm_dust_mass)= & pq(ig,l,igcm_ccnco2_mass) & /ptimestep+pdq(ig,l,igcm_ccnco2_mass) c$$$ c$$$ endif c write(*,*) "in co2clouds, riceco2(ig,l)v2= ",riceco2(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 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-4) 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)* & (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