subroutine improvedclouds(ngrid,nlay,ptimestep, & pplev,pplay,pt,pdt, & pq,pdq,pdqcloud,pdqscloud,pdtcloud, & nq,tauscaling,rdust,rice,nuice, & rsedcloud,rhocloud) implicit none c------------------------------------------------------------------ c This routine is used to form 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: J.-B. Madeleine, based on the work by Franck Montmessin c (October 2011) c------------------------------------------------------------------ #include "dimensions.h" #include "dimphys.h" #include "comcstfi.h" #include "callkeys.h" #include "tracer.h" #include "comgeomfi.h" #include "dimradmars.h" #include "microphys.h" c------------------------------------------------------------------ c Inputs: INTEGER ngrid,nlay integer nq ! nombre de traceurs REAL ptimestep ! pas de temps physique (s) REAL pplev(ngrid,nlay+1),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(ngridmx) ! Convertion factor for qdust and Ndust REAL rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m) c Outputs: REAL rice(ngrid,nlay) ! Ice mass mean radius (m) ! (r_c in montmessin_2004) REAL nuice(ngrid,nlay) ! Estimated effective variance ! of the size distribution REAL rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius REAL rhocloud(ngridmx,nlayermx) ! Cloud density (kg.m-3) REAL pdqcloud(ngrid,nlay,nq) ! tendance de la condensation ! H2O(kg/kg.s-1) REAL pdqscloud(ngrid,nq) ! flux en surface (kg.m-2.s-1) REAL pdtcloud(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 CBRT EXTERNAL CBRT INTEGER ig,l,i REAL zq(ngridmx,nlayermx,nqmx) ! local value of tracers REAL zq0(ngridmx,nlayermx,nqmx) ! local initial value of tracers REAL zt(ngridmx,nlayermx) ! local value of temperature REAL zqsat(ngridmx,nlayermx) ! saturation REAL lw !Latent heat of sublimation (J.kg-1) REAL Cste REAL dMice ! mass of condensated ice REAL sumcheck REAL*8 ph2o ! Water vapor partial pressure (Pa) REAL*8 satu ! Water vapor saturation ratio over ice REAL*8 Mo,No REAL*8 dN,dM,newvap REAL*8 Rn, Rm, dev2 REAL*8 n_aer(nbin_cld) ! number conc. of particle/each size bin REAL*8 m_aer(nbin_cld) ! mass mixing ratio of particle/each size bin REAL*8 rate(nbin_cld) ! nucleation rate REAL*8 up,dwn,Ctot,gr,seq REAL*8 sig ! Water-ice/air surface tension (N.m) EXTERNAL sig c Parameters of the size discretization c used by the microphysical scheme DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m) DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m) DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.0001e-6 ! 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_cld(nbin_cld+1)! boundary values of each rad_cld bin (m) SAVE rb_cld DOUBLE PRECISION dr_cld(nbin_cld)! width of each rad_cld bin (m) DOUBLE PRECISION vol_cld(nbin_cld) ! particle volume for each bin (m3) REAL sigma_ice ! Variance of the ice and CCN distributions SAVE sigma_ice c---------------------------------- c some outputs for 1D -- TEMPORARY REAL satu_out(ngridmx,nlayermx) ! satu ratio for output REAL dN_out(ngridmx,nlayermx) ! mass variation for output REAL dM_out(ngridmx,nlayermx) ! number variation for output REAL dM_col(ngridmx) ! total mass condensed in column REAL dN_col(ngridmx) ! total mass condensed in column REAL Mcon_out(ngridmx,nlayermx) ! mass to be condensed (not dMice !!) REAL gr_out(ngridmx,nlayermx) ! for 1d output REAL newvap_out(ngridmx,nlayermx) ! for 1d output REAL Mdust_col(ngridmx) ! total column dust mass REAL Ndust_col(ngridmx) ! total column dust mass REAL Mccn_col(ngridmx) ! total column ccn mass REAL Nccn_col(ngridmx) ! total column ccn mass REAL rate_out(ngridmx,nlayermx) ! nucleation rate INTEGER count LOGICAL output_sca ! scavenging outputs flag for tests output_sca = .false. c---------------------------------- c---------------------------------- c------------------------------------------------------------------ IF (firstcall) THEN c Definition of the size grid c ~~~~~~~~~~~~~~~~~~~~~~~~~~~ c rad_cld 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_cld array contains the boundary values of each rad_cld bin. c dr_cld is the width of each rad_cld bin. c Volume ratio between two adjacent bins vrat_cld = dlog(rmax_cld/rmin_cld) / float(nbin_cld-1) *3. vrat_cld = dexp(vrat_cld) write(*,*) "vrat_cld", vrat_cld rb_cld(1) = rbmin_cld rad_cld(1) = rmin_cld vol_cld(1) = 4./3. * dble(pi) * rmin_cld**3. do i=1,nbin_cld-1 rad_cld(i+1) = rad_cld(i) * vrat_cld**(1./3.) vol_cld(i+1) = vol_cld(i) * vrat_cld enddo do i=1,nbin_cld rb_cld(i+1)= ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) * & rad_cld(i) dr_cld(i) = rb_cld(i+1) - rb_cld(i) enddo rb_cld(nbin_cld+1) = rbmax_cld dr_cld(nbin_cld) = rb_cld(nbin_cld+1) - rb_cld(nbin_cld) print*, ' ' print*,'Microphysics: size bin information:' print*,'i,rb_cld(i), rad_cld(i),dr_cld(i)' print*,'-----------------------------------' do i=1,nbin_cld write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i), & dr_cld(i) enddo write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1) print*,'-----------------------------------' c Contact parameter of water ice on dust ( m=cos(theta) ) c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! mteta = 0.95 write(*,*) 'water_param contact parameter:', mteta c Volume of a water molecule (m3) vo1 = mh2o / dble(rho_ice) c Variance of the ice and CCN distributions sigma_ice = sqrt(log(1.+nuice_sed)) write(*,*) 'Variance of ice & CCN distribs :', sigma_ice write(*,*) 'nuice for sedimentation:', nuice_sed write(*,*) 'Volume of a water molecule:', vo1 firstcall=.false. END IF c----------------------------------------------------------------------- c 1. Initialization c----------------------------------------------------------------------- c Initialize the tendencies pdqscloud(1:ngrid,1:nq)=0 pdqcloud(1:ngrid,1:nlay,1:nq)=0 pdtcloud(1:ngrid,1:nlay)=0 c Update the needed variables do l=1,nlay do ig=1,ngrid c Temperature zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep c Dust mass mixing ratio zq(ig,l,igcm_dust_mass) = & pq(ig,l,igcm_dust_mass) + & pdq(ig,l,igcm_dust_mass) * ptimestep zq0(ig,l,igcm_dust_mass)=zq(ig,l,igcm_dust_mass) c Dust particle number zq(ig,l,igcm_dust_number) = & pq(ig,l,igcm_dust_number) + & pdq(ig,l,igcm_dust_number) * ptimestep zq0(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number) c Update rdust from last tendencies rdust(ig,l)= & CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/ & max(zq(ig,l,igcm_dust_number),0.01)) rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6) c CCN mass mixing ratio zq(ig,l,igcm_ccn_mass)= & pq(ig,l,igcm_ccn_mass)+pdq(ig,l,igcm_ccn_mass)*ptimestep zq(ig,l,igcm_ccn_mass)=max(zq(ig,l,igcm_ccn_mass),1.E-30) zq0(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass) c CCN particle number zq(ig,l,igcm_ccn_number)= & pq(ig,l,igcm_ccn_number)+pdq(ig,l,igcm_ccn_number)*ptimestep zq(ig,l,igcm_ccn_number)=max(zq(ig,l,igcm_ccn_number),1.E-30) zq0(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number) c Water vapor zq(ig,l,igcm_h2o_vap)= & pq(ig,l,igcm_h2o_vap)+pdq(ig,l,igcm_h2o_vap)*ptimestep zq(ig,l,igcm_h2o_vap)=max(zq(ig,l,igcm_h2o_vap),1.E-30) ! FF 12/2004 zq0(ig,l,igcm_h2o_vap)=zq(ig,l,igcm_h2o_vap) c Water ice zq(ig,l,igcm_h2o_ice)= & pq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep zq(ig,l,igcm_h2o_ice)=max(zq(ig,l,igcm_h2o_ice),0.) ! FF 12/2004 zq0(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice) enddo enddo c------------------------------------------------------------------ c Cloud microphysical scheme c------------------------------------------------------------------ Cste = ptimestep * 4. * pi * rho_ice call watersat(ngridmx*nlayermx,zt,pplay,zqsat) count = 0 c Main loop over the GCM's grid DO l=1,nlay DO ig=1,ngrid c Get the partial pressure of water vapor and its saturation ratio ph2o = zq(ig,l,igcm_h2o_vap) * (44./18.) * pplay(ig,l) satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l) IF ((satu .ge. 1)! ) THEN ! if we have condensation & .or. ( zq(ig,l,igcm_ccn_number) & .ge. 10) ) THEN ! or sublimation 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)+ 1.e-30 Rn = rdust(ig,l) Rm = Rn * exp( 3. * sigma_ice**2. ) Rn = 1. / Rn Rm = 1. / Rm dev2 = 1. / ( sqrt(2.) * sigma_ice ) do i = 1, nbin_cld n_aer(i) = 0.5 * No * ( derf( dlog(rb_cld(i+1)*Rn) * dev2 ) & -derf( dlog(rb_cld(i) * Rn) * dev2 ) ) m_aer(i) = 0.5 * Mo * ( derf( dlog(rb_cld(i+1)*Rm) * dev2 ) & -derf( dlog(rb_cld(i) * Rm) * dev2 ) ) enddo ! sumcheck = 0 ! do i = 1, nbin_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 ! endif ! ! sumcheck = 0 ! do i = 1, nbin_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 ! endif c Get the rates of nucleation call nuclea(ph2o,zt(ig,l),satu,n_aer,rate) dN = 0. dM = 0. rate_out(ig,l) = 0 do i = 1, nbin_cld n_aer(i) = n_aer(i) / ( 1. + rate(i)*ptimestep ) m_aer(i) = m_aer(i) / ( 1. + rate(i)*ptimestep ) dN = dN + n_aer(i) * rate(i) * ptimestep dM = dM + m_aer(i) * rate(i) * ptimestep rate_out(ig,l)=rate_out(ig,l)+rate(i) enddo ! dN = min( max(dN,-zq(ig,l,igcm_ccn_number) ), ! & zq(ig,l,igcm_dust_number) ) ! ! dM = min( max(dM,-zq(ig,l,igcm_ccn_mass) ), ! & zq(ig,l,igcm_dust_mass) ) Mo = zq0(ig,l,igcm_h2o_ice) + & zq0(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30 No = zq0(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30 !write(*,*) "l,cloud particles,cloud mass",l, No, Mo rhocloud(ig,l) = zq0(ig,l,igcm_h2o_ice) / Mo * rho_ice & +zq0(ig,l,igcm_ccn_mass)* tauscaling(ig) & / Mo * rho_dust rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust) rice(ig,l) = & ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.) c nuice(ig,l)=nuice_ref ! used for rad. transfer calculations if ((Mo.lt.1.e-20) .or. (No.le.1)) rice(ig,l) = 1.e-8 seq = exp( 2.*sig(zt(ig,l))*mh2o / & (rho_ice*rgp*zt(ig,l)*rice(ig,l)) ) call growthrate(ptimestep,zt(ig,l),pplay(ig,l), & ph2o,ph2o/satu,seq,rice(ig,l),gr) up = Cste * gr * rice(ig,l) * No * seq + & zq(ig,l,igcm_h2o_vap) dwn = Cste * gr * rice(ig,l) * No / zqsat(ig,l)+ 1. Ctot = zq0(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_h2o_vap) newvap = min(up/dwn,Ctot) gr = gr * ( newvap/zqsat(ig,l) - seq ) dMice = min( max(Cste * No * rice(ig,l) * gr, & -zq(ig,l,igcm_h2o_ice) ), & zq(ig,l,igcm_h2o_vap) ) c----------- TESTS 1D output --------- satu_out(ig,l) = satu Mcon_out(ig,l) = dMice newvap_out(ig,l) = newvap gr_out(ig,l) = gr dN_out(ig,l) = dN dM_out(ig,l) = dM c------------------------------------- c Water ice zq(ig,l,igcm_h2o_ice) = zq0(ig,l,igcm_h2o_ice) + & dMice c Water vapor zq(ig,l,igcm_h2o_vap) = zq0(ig,l,igcm_h2o_vap) - & dMice c If all the ice particles sublimate, all the condensation c nuclei are released: if (zq(ig,l,igcm_h2o_ice).le.1e-30) then c for coherence dM = 0 dN = 0 dN_out(ig,l) = - zq(ig,l,igcm_ccn_number)*tauscaling(ig) dM_out(ig,l) = - zq(ig,l,igcm_ccn_mass)*tauscaling(ig) c Water ice particles zq(ig,l,igcm_h2o_ice) = 0. c Dust particles zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + & zq(ig,l,igcm_ccn_mass) zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + & zq(ig,l,igcm_ccn_number) c CCNs zq(ig,l,igcm_ccn_mass) = 0. zq(ig,l,igcm_ccn_number) = 0. endif dN = dN/ tauscaling(ig) dM = dM/ tauscaling(ig) c Dust particles zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - dM zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - dN c CCNs zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + dM zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + dN pdqcloud(ig,l,igcm_dust_mass)=(zq(ig,l,igcm_dust_mass) & -zq0(ig,l,igcm_dust_mass))/ptimestep pdqcloud(ig,l,igcm_dust_number)=(zq(ig,l,igcm_dust_number) & -zq0(ig,l,igcm_dust_number))/ptimestep pdqcloud(ig,l,igcm_ccn_mass)=(zq(ig,l,igcm_ccn_mass) & -zq0(ig,l,igcm_ccn_mass))/ptimestep pdqcloud(ig,l,igcm_ccn_number)=(zq(ig,l,igcm_ccn_number) & -zq0(ig,l,igcm_ccn_number))/ptimestep pdqcloud(ig,l,igcm_h2o_vap)=(zq(ig,l,igcm_h2o_vap) & -zq0(ig,l,igcm_h2o_vap))/ptimestep pdqcloud(ig,l,igcm_h2o_ice)=(zq(ig,l,igcm_h2o_ice) & -zq0(ig,l,igcm_h2o_ice))/ptimestep count = count +1 ELSE ! for coherence (rdust, rice computations etc ..) zq(ig,l,igcm_dust_mass) = zq0(ig,l,igcm_dust_mass) zq(ig,l,igcm_dust_number) = zq0(ig,l,igcm_dust_number) zq(ig,l,igcm_ccn_mass) = zq0(ig,l,igcm_ccn_mass) zq(ig,l,igcm_ccn_number) = zq0(ig,l,igcm_ccn_number) zq(ig,l,igcm_h2o_ice) = zq0(ig,l,igcm_h2o_ice) zq(ig,l,igcm_h2o_vap) = zq0(ig,l,igcm_h2o_vap) ! pour les sorties de test satu_out(ig,l) = satu Mcon_out(ig,l) = 0 newvap_out(ig,l) = zq(ig,l,igcm_h2o_vap) gr_out(ig,l) = 0 dN_out(ig,l) = 0 dM_out(ig,l) = 0 ENDIF ! end if (saturation ratio > 1) or (there is h2o_ice) c-----update temperature lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)**2)*1.e+3 pdtcloud(ig,l)=-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp c----- update rice & rhocloud AFTER scavenging Mo = zq(ig,l,igcm_h2o_ice) + & zq(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30 No = zq(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30 rhocloud(ig,l) = zq(ig,l,igcm_h2o_ice) / Mo * rho_ice & +zq(ig,l,igcm_ccn_mass)* tauscaling(ig) & / Mo * rho_dust rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust) rice(ig,l) = & ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.) if ((Mo.lt.1.e-20) .or. (No.le.1)) rice(ig,l) = 1.e-8 nuice(ig,l)=nuice_ref ! used for rad. transfer calculations c----- update rdust and sedimentation radius rdust(ig,l)= & CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/ & max(zq(ig,l,igcm_dust_number),0.01)) rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6) rsedcloud(ig,l)=max( rice(ig,l)*(1.+nuice_sed)**3., & rdust(ig,l) ) rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4) ENDDO ENDDO c------------------------------------------------------------------ c------------------------------------------------------------------ c------------------------------------------------------------------ c------------------------------------------------------------------ c------------------------------------------------------------------ c TESTS IF (output_sca) then print*, 'count is ',count, ' i.e. ', & count*100/(nlay*ngrid), '% for microphys computation' dM_col(:) = 0 dN_col(:) = 0 Mdust_col(:) = 0 Ndust_col(:) = 0 Mccn_col(:) = 0 Nccn_col(:) = 0 do l=1, nlay do ig=1,ngrid dM_col(ig) = dM_col(ig) + & dM_out(ig,l)*(pplev(ig,l) - pplev(ig,l+1)) / g dN_col(ig) = dN_col(ig) + & dN_out(ig,l)*(pplev(ig,l) - pplev(ig,l+1)) / g Mdust_col(ig) = Mdust_col(ig) + & zq(ig,l,igcm_dust_mass)*tauscaling(ig) & *(pplev(ig,l) - pplev(ig,l+1)) / g Ndust_col(ig) = Ndust_col(ig) + & zq(ig,l,igcm_dust_number)*tauscaling(ig) & *(pplev(ig,l) - pplev(ig,l+1)) / g Mccn_col(ig) = Mccn_col(ig) + & zq(ig,l,igcm_ccn_mass)*tauscaling(ig) & *(pplev(ig,l) - pplev(ig,l+1)) / g Nccn_col(ig) = Nccn_col(ig) + & zq(ig,l,igcm_ccn_number)*tauscaling(ig) & *(pplev(ig,l) - pplev(ig,l+1)) / g enddo ! of do ig=1,ngrid enddo ! of do l=1,nlay IF (ngrid.ne.1) THEN ! 3D call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3, & satu_out) call WRITEDIAGFI(ngrid,"dM_col","dM column","kg",2, & dM_col) call WRITEDIAGFI(ngrid,"dN_col","dN column","N",2, & dN_col) call WRITEDIAGFI(ngrid,"Ndust_col","M column","N",2, & Ndust_col) call WRITEDIAGFI(ngrid,"Mdust_col","N column","kg",2, & Mdust_col) call WRITEDIAGFI(ngrid,"Nccn_col","M column","N",2, & Nccn_col) call WRITEDIAGFI(ngrid,"Mccn_col","N column","kg",2, & Mccn_col) call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3, & dM_out) call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3, & dN_out) ENDIF IF (ngrid.eq.1) THEN ! 1D call WRITEDIAGFI(ngrid,"newvap","h2o newvap","kg",1, & newvap_out) call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1, & gr_out) call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1, & rate_out) call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1, & dM_out) call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1, & dN_out) call WRITEDIAGFI(ngrid,"mcond","h2o condensed mass","kg",1, & Mcon_out) call WRITEDIAGFI(ngrid,"zqsat","p vap sat","kg/kg",1, & zqsat) call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1, & satu_out) call WRITEDIAGFI(ngrid,"rice_sca","ice radius","m",1, & rice) 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) call WRITEDIAGFI(ngrid,"dM_col","dM column","kg",0, & dM_col) call WRITEDIAGFI(ngrid,"dN_col","dN column","N",0, & dN_col) call WRITEDIAGFI(ngrid,"Ndust_col","M column","N",0, & Ndust_col) call WRITEDIAGFI(ngrid,"Mdust_col","N column","kg",0, & Mdust_col) call WRITEDIAGFI(ngrid,"Nccn_col","M column","N",0, & Nccn_col) call WRITEDIAGFI(ngrid,"Mccn_col","N column","kg",0, & Mccn_col) ENDIF ENDIF ! endif output_sca c------------------------------------------------------------------ return end