subroutine improvedclouds(ngrid,nlay,ptimestep, & pplev,pplay,pt,pdt, & pq,pdq,pdqcloud,pdtcloud, & nq,tauscaling,rdust,rice,nuice, & rsedcloud,rhocloud) ! to use 'getin' USE ioipsl_getincom 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 A word about the radius growth ... c A word about nucleation and ice growth strategies ... c Authors: J.-B. Madeleine, based on the work by Franck Montmessin c (October 2011) c T. Navarro, debug,correction, new scheme (October-April 2011) c A. Spiga, optimization (February 2012) c------------------------------------------------------------------ #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) ! pression aux inter-couches (Pa) REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa) REAL pt(ngrid,nlay) ! temperature at the middle of the ! layers (K) REAL pdt(ngrid,nlay) ! tendance temperature des autres ! param. REAL pq(ngrid,nlay,nq) ! traceur (kg/kg) REAL pdq(ngrid,nlay,nq) ! tendance avant condensation ! (kg/kg.s-1) REAL tauscaling(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 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, yeah, n_derf, m_derf 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 REAl*8 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 REAL tdicho, tmax, rmin, rmax, req, rdicho REAL coeff0, coeff1, coeff2 REAL error_out(ngridmx,nlayermx) REAL error2d(ngridmx) REAL var1,var2,var3 REAL rccn, epsilon c---------------------------------- c some outputs for 1D -- TESTS 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 Mcon_out(ngridmx,nlayermx) ! mass to be condensed (not dMice !!) REAL gr_out(ngridmx,nlayermx) ! for 1d output REAL rice_out(ngridmx,nlayermx) ! ice radius change REAL rate_out(ngridmx,nlayermx) ! nucleation rate REAL satubf(ngridmx,nlayermx),satuaf(ngridmx,nlayermx) REAL ccnbf(ngridmx,nlayermx),ccnaf(ngridmx,nlayermx) INTEGER count LOGICAL test_flag ! flag for test/debuging outputs test_flag = .false. c---------------------------------- c---------------------------------- c------------------------------------------------------------------ IF (firstcall) THEN !============================================================= ! 0. Definition of the size grid !============================================================= 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*rmin_cld*rmin_cld 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*,'-----------------------------------' do i=1,nbin_cld+1 rb_cld(i) = dlog(rb_cld(i)) !! we save that so that it is not computed !! at each timestep and gridpoint enddo 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 !============================================================= ! 1. Initialisation !============================================================= c Initialize the tendencies 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),1e-30) ! FF 12/2004 zq0(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice) enddo enddo !============================================================= ! 2. Compute saturation !============================================================= dev2 = 1. / ( sqrt(2.) * sigma_ice ) error_out(:,:) = 0. 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. ) .or. ! if there is condensation & ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.)) THEN ! or sublimation !============================================================= ! 3. Nucleation !============================================================= c Expand the dust moments into a binned distribution Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig) !+ 1.e-30 No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30 Rn = rdust(ig,l) Rn = -dlog(Rn) Rm = Rn - 3. * sigma_ice*sigma_ice n_derf = derf( (rb_cld(1)+Rn) *dev2) m_derf = derf( (rb_cld(1)+Rm) *dev2) do i = 1, nbin_cld n_aer(i) = -0.5 * No * n_derf !! this ith previously computed m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed n_derf = derf( (rb_cld(i+1)+Rn) *dev2) m_derf = derf( (rb_cld(i+1)+Rm) *dev2) n_aer(i) = n_aer(i) + 0.5 * No * n_derf m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf enddo !!! MORE EFFICIENT COMPUTATIONALLY THAN ! Rn = rdust(ig,l) ! Rm = Rn * exp( 3. * sigma_ice**2. ) ! Rn = 1. / Rn ! Rm = 1. / Rm ! 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 c Update CCNs, can also be done after the radius growth c Dust particles zq(ig,l,igcm_dust_mass) = & zq(ig,l,igcm_dust_mass) - dM/ tauscaling(ig) zq(ig,l,igcm_dust_number) = & zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig) c CCNs zq(ig,l,igcm_ccn_mass) = & zq(ig,l,igcm_ccn_mass) + dM/ tauscaling(ig) zq(ig,l,igcm_ccn_number) = & zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig) !============================================================= ! 4. Ice growth: scheme for radius evolution !============================================================= 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) c nuclei radius rccn = CBRT(zq(ig,l,igcm_ccn_mass)/ & (pi*rho_dust*zq(ig,l,igcm_ccn_number)*4./3.)) c ice crystal radius rice (ig,l) = & CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) ) c enforce physical value : crystal cannot be smaller than its nuclei ! rice(ig,l) = max(rice(ig,l), rccn) c saturation at equilibrium seq = exp( 2.*sig(zt(ig,l))*mh2o / & (rho_ice*rgp*zt(ig,l)*rice(ig,l)) ) c If there is more than on nuclei, we peform ice growth var1 = zq0(ig,l,igcm_ccn_number)*tauscaling(ig) + dN IF (var1 .ge. -1) THEN if (test_flag) then print*, ' ' print*, ptimestep print*, 'satu,seq', real(satu), real(seq), ig,l print*, 'dN,dM', real(dN),real(dM) print*,'rccn', rccn print*, 'Nccn, Mccn', zq(ig,l,igcm_ccn_number)*tauscaling(ig), & zq(ig,l,igcm_ccn_mass)*tauscaling(ig) endif c crystal radius to reach saturation at equilibrium (i.e. satu=seq) req = ( zq0(ig,l,igcm_h2o_ice) + zq0(ig,l,igcm_h2o_vap) & + zq(ig,l,igcm_ccn_mass)*tauscaling(ig)* & rho_ice/rho_dust - seq * zqsat(ig,l)) & / ( zq(ig,l,igcm_ccn_number)*tauscaling(ig)* & pi*rho_ice*4./3. ) req = CBRT(req) c compute ice radius growth resistances (diffusive and latent heat resistancea) call growthrate(zt(ig,l),pplay(ig,l), & ph2o/satu,seq,req,coeff1,coeff2) coeff0 = -zqsat(ig,l) / (4.*pi*req*rho_ice & * zq(ig,l,igcm_ccn_number)*tauscaling(ig)) c compute tmax, the time needed to reach req call phi(rice(ig,l),req,coeff1,coeff2,tmax) if (test_flag) then print*, 'coeffs', coeff0,coeff1,coeff2 print*, 'req,tmax', req,tmax*coeff0 print*, 'i,rmin,rdicho,rmax,tdicho' endif c rmin is rice if r increases (satu >1) or req if it decreases (satu<1) c if req is lower than rccn (ie not enough ice to reach saturation), rmin is forced to rccn if (satu .ge. seq) then ! crystal size is increasing rmin = max(min(rice(ig,l),req),rccn) rmax = max(rice(ig,l),req) else ! crystal size is decreasing rmax = max(min(rice(ig,l),req),rccn) rmin = max(rice(ig,l),req) endif !rmax = min(rmax,1.e-3) ! au max on a des rayons de 1 mm pour la dicho ... rdicho = 0.5*(rmin+rmax) ! for output var1 = rice(ig,l) var2 = rmin var3 = rmax c Given the phi function is monotonous, we perform a simple dichotomy to find the raidus at t+1 do i = 1,10 ! dichotomy loop c compute tdicho, the time needed to reach rdicho call phi(rdicho,req,coeff1,coeff2,tdicho) !print*, tdicho,tmax tdicho = coeff0*(tdicho - tmax) if (test_flag) print*, i,rmin,rdicho,rmax,tdicho if (tdicho .ge. ptimestep) then rmax = rdicho else rmin = rdicho endif rdicho = 0.5*(rmin+rmax) enddo ! of dichotomy loop if (test_flag) then if (abs(rdicho - rccn) .ge. 1e-15) then ! to avoid infinite values epsilon = (rmax - rmin)/(2**10) error_out(ig,l) = & 100*(epsilon**3 +3*epsilon**2*rdicho +3*epsilon*rdicho**2) & / (rdicho**3-rccn**3) endif print*, 'error masse glace %', error_out(ig,l) print*, 'rice,ice,vap bf', & rice(ig,l),zq0(ig,l,igcm_h2o_ice),zq0(ig,l,igcm_h2o_vap) endif c If the initial condition is subsaturated and there is not enough ice available for sublimation c to reach equilibrium, req is neagtive. Therefore, enforce physical value. rice(ig,l) = max(rdicho,rccn) !!!!! Water ice mass is computed with radius at t+1 and their current number !!!!! Nccn is at t or t+1, depending on what has been done before ! zq(ig,l,igcm_h2o_ice) = ! & (pi*rho_ice*zq(ig,l,igcm_ccn_number)*4./3. ! & * rdicho*rdicho*rdicho - ! & zq(ig,l,igcm_ccn_mass)*rho_ice/rho_dust) ! & * tauscaling(ig) !!!!! Water ice mass is computed with radius at t+1 and their number at t+1 !!!!! that is dirty, but this enforces the use of Nccn at t+1, whatever is done before !!!!! TO BE CLEANED ONE DAY var1 = zq0(ig,l,igcm_ccn_number)*tauscaling(ig) + dN var2 = zq0(ig,l,igcm_ccn_mass)*tauscaling(ig) + dM zq(ig,l,igcm_h2o_ice) = & (pi*rho_ice*var1*4./3. & * rdicho*rdicho*rdicho - & var2*rho_ice/rho_dust) !!!! enforce realistic values (if the case of growth on Nccn(t) and condensation on Nccn(t+1)) zq(ig,l,igcm_h2o_ice) = & min(zq0(ig,l,igcm_h2o_ice) + zq0(ig,l,igcm_h2o_vap), & zq(ig,l,igcm_h2o_ice)) zq(ig,l,igcm_h2o_ice) = & max(1e-30,zq(ig,l,igcm_h2o_ice)) zq(ig,l,igcm_h2o_vap) = & zq0(ig,l,igcm_h2o_ice) + zq0(ig,l,igcm_h2o_vap) & - zq(ig,l,igcm_h2o_ice) if (test_flag) then print*, 'rice,ice,vap af', & rice(ig,l),zq(ig,l,igcm_h2o_ice),zq(ig,l,igcm_h2o_vap) print*, 'satu bf, af', & zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l), & zq(ig,l,igcm_h2o_vap)/zqsat(ig,l) endif !!!!!!!!!!!! TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST if ((zq(ig,l,igcm_h2o_ice).le. -1e-8) & .or. (zq(ig,l,igcm_h2o_vap).le. -1e-8)) then print*, 'NEGATIVE WATER' print*, 'ig,l', ig,l print*, 'satu', satu print*, 'vap, ice bf', & zq0(ig,l,igcm_h2o_vap), zq0(ig,l,igcm_h2o_ice) print*, 'vap, ice af', & zq(ig,l,igcm_h2o_vap), zq(ig,l,igcm_h2o_ice) print*, 'ccn N,M bf', & zq0(ig,l,igcm_ccn_number), zq0(ig,l,igcm_ccn_mass) print*, 'ccn N,M af', & zq(ig,l,igcm_ccn_number), zq(ig,l,igcm_ccn_mass) print*, 'tauscaling', & tauscaling(ig) print*, 'req,rccn,rice bf,rice af', & req,rccn,var1,rice(ig,l) print*, 'rmin, rmax', var2,var3 print*, 'error_out,tdicho,ptimestep', & error_out(ig,l),tdicho,ptimestep print*, ' ' endif !!!!!!!!!!!! TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST ENDIF !of if Nccn >1 !============================================================= ! 5. Dust cores released, tendancies, latent heat, etc ... !============================================================= c If all the ice particles sublimate, all the condensation c nuclei are released: if (zq(ig,l,igcm_h2o_ice).le.1e-10) then ! 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 zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap) & + zq(ig,l,igcm_h2o_ice) 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 ! 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) ! ccnbf(ig,l) = zq0(ig,l,igcm_ccn_number) * tauscaling(ig) ! ccnaf(ig,l) = zq(ig,l,igcm_ccn_number) * tauscaling(ig) ! ! satubf(ig,l) = zq0(ig,l,igcm_h2o_vap) / zqsat(ig,l) ! satuaf(ig,l) = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l) c-----update temperature (latent heat relase) 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 microphysic 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) if ((Mo.lt.1.e-20) .or. (No.le.1)) then rice(ig,l) = 1.e-8 else rice(ig,l) = & CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) ) !**(1./3.) endif 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)* & (1.+nuice_sed)*(1.+nuice_sed), & rdust(ig,l) ) rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4) ENDDO ! of ig loop ENDDO ! of nlayer loop !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS IF (test_flag) then error2d(:) = 0. DO l=1,nlay DO ig=1,ngrid error2d(ig) = max(abs(error_out(ig,l)),error2d(ig)) ENDDO ENDDO print*, 'count is ',count, ' i.e. ', & count*100/(nlay*ngrid), '% for microphys computation' IF (ngrid.ne.1) THEN ! 3D ! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3, ! & satu_out) ! call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3, ! & dM_out) ! call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3, ! & dN_out) call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2, & error2d) ! call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3, ! & zqsat) ENDIF IF (ngrid.eq.1) THEN ! 1D call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1, & error_out) call WRITEDIAGFI(ngrid,"satu_bf","satu before","kg/kg",1, & satubf) call WRITEDIAGFI(ngrid,"satu_af","satu after","kg/kg",1, & satuaf) call WRITEDIAGFI(ngrid,"vapbf","h2ovap before","kg/kg",1, & zq0(1,:,igcm_h2o_vap)) call WRITEDIAGFI(ngrid,"vapaf","h2ovap after","kg/kg",1, & zq(1,:,igcm_h2o_vap)) call WRITEDIAGFI(ngrid,"icebf","h2oice before","kg/kg",1, & zq0(1,:,igcm_h2o_ice)) call WRITEDIAGFI(ngrid,"iceaf","h2oice after","kg/kg",1, & zq(1,:,igcm_h2o_ice)) call WRITEDIAGFI(ngrid,"ccnbf","ccn before","/kg",1, & ccnbf) call WRITEDIAGFI(ngrid,"ccnaf","ccn after","/kg",1, & ccnaf) c call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1, c & gr_out) c call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1, c & rate_out) c call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1, c & dM_out) c call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1, c & dN_out) 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) ENDIF ENDIF ! endif test_flag !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0 c It is an analytical solution to the ice radius growth equation, c with the approximation of a constant 'reduced' cunningham correction factor c (lambda in growthrate.F) taken at radius req instead of rice cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine phi(rice,req,coeff1,coeff2,time) implicit none ! inputs real rice ! ice radius real req ! ice radius at equilibirum real coeff1 ! coeff for the log real coeff2 ! coeff for the arctan ! output real time !local real var ! 1.73205 is sqrt(3) var = max( & abs(rice-req) / sqrt(rice*rice + rice*req + req*req),1e-30) time = & coeff1 * & log( var ) & + coeff2 * 1.73205 * & atan( (2*rice+req) / (1.73205*req) ) return end