source: trunk/LMDZ.MARS/libf/phymars/improvedCO2clouds_mod.F @ 1962

Last change on this file since 1962 was 1939, checked in by aslmd, 7 years ago

rename improvedCO2clouds in improvedCO2clouds_mod since it is a module

  • Property svn:executable set to *
File size: 29.6 KB
RevLine 
[1938]1      MODULE improvedCO2clouds_mod
2
3      IMPLICIT NONE
4
5      CONTAINS
6
[1911]7      subroutine improvedCO2clouds(ngrid,nlay,microtimestep,
8     &             pplay,pplev,pteff,sum_subpdt,
9     &             pqeff,sum_subpdq,subpdqcloudco2,subpdtcloudco2,
[1685]10     &             nq,tauscaling,
[1922]11     &             mem_Mccn_co2,mem_Mh2o_co2,mem_Nccn_co2)
[1884]12      USE comcstfi_h, only: pi, g, cpp
13      USE updaterad, only: updaterice_micro, updaterice_microco2
14      use tracer_mod, only: igcm_dust_mass, igcm_dust_number, rho_dust,
15     &                      igcm_h2o_ice, igcm_ccn_mass,
16     &                      igcm_ccn_number, nuice_sed,
17     &                      igcm_co2, igcm_co2_ice, igcm_ccnco2_mass,
18     &                      igcm_ccnco2_number, nuiceco2_sed,
19     &                      rho_ice_co2
[1617]20      use conc_mod, only: mmean
[1918]21      use datafile_mod, only: datadir
[1720]22
[1617]23      implicit none
24     
25c------------------------------------------------------------------
26c  This routine is used to form CO2 clouds when a parcel of the GCM is
27c    saturated. It includes the ability to have supersaturation, a
28c    computation of the nucleation rates, growthrates and the
29c    scavenging of dust particles by clouds.
30c  It is worth noting that the amount of dust is computed using the
31c    dust optical depth computed in aeropacity.F. That's why
32c    the variable called "tauscaling" is used to convert
33c    pq(dust_mass) and pq(dust_number), which are relative
34c    quantities, to absolute and realistic quantities stored in zq.
35c    This has to be done to convert the inputs into absolute
36c    values, but also to convert the outputs back into relative
37c    values which are then used by the sedimentation and advection
38c    schemes.
[1685]39c CO2 ice particles can nucleate on both dust and on water ice particles
40c When CO2 ice is deposited onto a water ice particles, the particle is
41c removed from the water tracers.
[1818]42c Memory of the origin of the co2 particles is kept and thus the
43c water cycle shouldn't be modified by this.
[1922]44c WARNING: no sedimentation of the water ice origin is performed
[1685]45c in the microphysical timestep in co2cloud.F.
[1617]46
47c  Authors of the water ice clouds microphysics
48c J.-B. Madeleine, based on the work by Franck Montmessin
49c           (October 2011)
50c           T. Navarro, debug,correction, new scheme (October-April 2011)
51c           A. Spiga, optimization (February 2012)
52c Adaptation for CO2 clouds by Joachim Audouard (09/16), based on the work
53c of Constantino Listowski
[1818]54c There is an energy limit to how much co2 can sublimate/condensate. It is
55c defined by the difference of the GCM temperature with the co2 condensation
56c temperature.
[1816]57c Warning:
58c If meteoritic particles are activated and turn into co2 ice particles,
59c then they will be reversed in the dust tracers if the cloud sublimates
60 
[1617]61c------------------------------------------------------------------
[1884]62      include "callkeys.h"
63      include "microphys.h"
[1617]64c------------------------------------------------------------------
[1884]65c     Arguments:
[1617]66
[1884]67      INTEGER,INTENT(in) :: ngrid,nlay
68      integer,intent(in) :: nq         ! number of tracers
[1911]69      REAL,INTENT(in) :: microtimestep     ! physics time step (s)
[1884]70      REAL,INTENT(in) :: pplay(ngrid,nlay)     ! mid-layer pressure (Pa)
71      REAL,INTENT(in) :: pplev(ngrid,nlay+1)   ! inter-layer pressure (Pa)
[1911]72      REAL,INTENT(in) :: pteff(ngrid,nlay) ! temperature at the middle of the
[1617]73                                 !   layers (K)
[1911]74      REAL,INTENT(in) :: sum_subpdt(ngrid,nlay) ! tendency on temperature from
[1884]75                                 !  previous physical parametrizations
[1911]76      REAL,INTENT(in) :: pqeff(ngrid,nlay,nq) ! tracers (kg/kg)
77      REAL,INTENT(in) :: sum_subpdq(ngrid,nlay,nq) ! tendencies on tracers
[1884]78                                 !  before condensation (kg/kg.s-1)
79      REAL,INTENT(in) :: tauscaling(ngrid) ! Convertion factor for qdust and Ndust
[1617]80c     Outputs:
[1911]81      REAL,INTENT(out) :: subpdqcloudco2(ngrid,nlay,nq) ! tendency on tracers
[1884]82                                   ! due to CO2 condensation (kg/kg.s-1)
[1617]83      ! condensation si igcm_co2_ice
[1911]84      REAL,INTENT(out) :: subpdtcloudco2(ngrid,nlay)  ! tendency on temperature due
[1884]85                                   ! to latent heat
[1617]86
87c------------------------------------------------------------------
88c     Local variables:
[1884]89      LOGICAL,SAVE :: firstcall=.true.
[1617]90      REAL*8   derf ! Error function
[1921]91      INTEGER ig,l,i
92
[1884]93      real masse (ngrid,nlay) ! Layer mass (kg.m-2)
94      REAL rice(ngrid,nlay)    ! Water Ice mass mean radius (m)
95                                ! used for nucleation of CO2 on ice-coated ccns
96      REAL rccnh2o(ngrid,nlay)    ! Water Ice mass mean radius (m)
[1617]97      REAL zq(ngrid,nlay,nq)  ! local value of tracers
98      REAL zq0(ngrid,nlay,nq) ! local initial value of tracers
99      REAL zt(ngrid,nlay)       ! local value of temperature
100      REAL zqsat(ngrid,nlay)    ! saturation vapor pressure for CO2
[1816]101      real tcond(ngrid,nlay)
102      real zqco2(ngrid,nlay)
[1617]103      REAL lw                         !Latent heat of sublimation (J.kg-1)
[1720]104      REAL,save :: l0,l1,l2,l3,l4
[1685]105      DOUBLE PRECISION dMice           ! mass of condensed ice
[1617]106      DOUBLE PRECISION sumcheck
[1816]107      DOUBLE PRECISION facteurmax!for energy limit on mass growth
108      DOUBLE PRECISION pco2,psat  ! Co2 vapor partial pressure (Pa)
[1617]109      DOUBLE PRECISION satu ! Co2 vapor saturation ratio over ice
[1921]110      DOUBLE PRECISION Mo,No,No_dust,Mo_dust
[1720]111      DOUBLE PRECISION  Rn, Rm, dev2,dev3, n_derf, m_derf
[1922]112      DOUBLE PRECISION mem_Mccn_co2(ngrid,nlay) ! Memory of CCN mass of H2O and dust used by CO2
113      DOUBLE PRECISION mem_Mh2o_co2(ngrid,nlay) ! Memory of H2O mass integred into CO2 crystal
114      DOUBLE PRECISION mem_Nccn_co2(ngrid,nlay) ! Memory of CCN number of H2O and dust used by CO2
[1930]115      DOUBLE PRECISION interm1,interm2,interm3     
116 
[1617]117!     Radius used by the microphysical scheme (m)
118      DOUBLE PRECISION n_aer(nbinco2_cld) ! number concentration volume-1 of particle/each size bin
119      DOUBLE PRECISION m_aer(nbinco2_cld) ! mass mixing ratio of particle/each size bin
[1720]120      DOUBLE PRECISION m_aer_h2oice2(nbinco2_cld) ! mass mixing ratio of particle/each size bin
[1617]121
122      DOUBLE PRECISION n_aer_h2oice(nbinco2_cld) ! Same - for CO2 nucleation
123      DOUBLE PRECISION m_aer_h2oice(nbinco2_cld) ! Same - for CO2 nucleation
124      DOUBLE PRECISION rad_h2oice(nbinco2_cld)
125
126c      REAL*8 sigco2      ! Co2-ice/air surface tension  (N.m)
127c      EXTERNAL sigco2
128
[1685]129      DOUBLE PRECISION dN,dM, dNh2o, dMh2o, dNN,dMM,dNNh2o,dMMh2o
[1720]130      DOUBLE PRECISION dMh2o_ice,dMh2o_ccn
131
[1617]132      DOUBLE PRECISION rate(nbinco2_cld)  ! nucleation rate
133      DOUBLE PRECISION rateh2o(nbinco2_cld)  ! nucleation rate
134      REAL seq
[1685]135      DOUBLE PRECISION rho_ice_co2T(ngrid,nlay)
[1617]136      DOUBLE PRECISION riceco2(ngrid,nlay)      ! CO2Ice mean radius (m)
[1685]137      REAL rhocloud(ngrid,nlay) ! Cloud density (kg.m-3)
138                 
[1617]139      REAL rhocloudco2(ngrid,nlay)  ! Cloud density (kg.m-3)
140      REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m)
141
142c      REAL res      ! Resistance growth
143      DOUBLE PRECISION Ic_rice      ! Mass transfer rate CO2 ice crystal
[1685]144      DOUBLE PRECISION ratioh2o_ccn
145      DOUBLE PRECISION vo2co2
[1921]146
147c     Parameters of the size discretization used by the microphysical scheme
148      DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-9   ! Minimum radius (m)
149      DOUBLE PRECISION, PARAMETER :: rmax_cld = 5.e-6   ! Maximum radius (m)
150      DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-10  ! Minimum bounary radius (m)
151      DOUBLE PRECISION, PARAMETER :: rbmax_cld = 2.e-4  ! Maximum boundary radius (m)
152      DOUBLE PRECISION vrat_cld                         ! Volume ratio
153      DOUBLE PRECISION rb_cldco2(nbinco2_cld+1)         ! boundary values of each rad_cldco2 bin (m)
[1617]154      SAVE rb_cldco2
[1921]155      DOUBLE PRECISION dr_cld(nbinco2_cld)              ! width of each rad_cldco2 bin (m)
156      DOUBLE PRECISION vol_cld(nbinco2_cld)             ! particle volume for each bin (m3)
[1617]157
[1685]158      DOUBLE PRECISION Proba,Masse_atm,drsurdt,reff,Probah2o
[1921]159      REAL sigma_iceco2   ! Variance of the co2 ice and CCN distributions
[1816]160      SAVE sigma_iceco2
[1921]161      REAL sigma_ice      ! Variance of the h2o ice and CCN distributions
[1720]162      SAVE sigma_ice
[1632]163      DOUBLE PRECISION Niceco2,Qccnco2,Nccnco2
[1921]164
[1816]165!     Variables for the meteoritic flux:
[1884]166      integer,parameter :: nbin_meteor=100
167      integer,parameter :: nlev_meteor=130
168      double precision meteor_ccn(ngrid,nlay,100) !100=nbinco2_cld !!!
169      double precision,save :: meteor(130,100)
170      double precision mtemp(100),pression_meteor(130)
[1816]171      logical file_ok
[1884]172      integer read_ok
[1816]173      integer nelem,lebon1,lebon2
174      double precision :: ltemp1(130),ltemp2(130)
[1913]175      integer ibin,j
176      integer,parameter :: uMeteor=666
[1720]177
[1617]178      IF (firstcall) THEN
179!=============================================================
180! 0. Definition of the size grid
181!=============================================================
182c       rad_cldco2 is the primary radius grid used for microphysics computation.
183c       The grid spacing is computed assuming a constant volume ratio
184c       between two consecutive bins; i.e. vrat_cld.
185c       vrat_cld is determined from the boundary values of the size grid:
186c       rmin_cld and rmax_cld.
187c       The rb_cldco2 array contains the boundary values of each rad_cldco2 bin.
188c       dr_cld is the width of each rad_cldco2 bin.
189
[1816]190        vrat_cld = dlog(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3.
191        vrat_cld = dexp(vrat_cld)
[1617]192        rb_cldco2(1)  = rbmin_cld
193        rad_cldco2(1) = rmin_cld
194        vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld
195        do i=1,nbinco2_cld-1
196          rad_cldco2(i+1)  = rad_cldco2(i) * vrat_cld**(1./3.)
197          vol_cld(i+1)  = vol_cld(i) * vrat_cld
[1816]198        enddo       
[1617]199        do i=1,nbinco2_cld
200          rb_cldco2(i+1)= ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) *
201     &      rad_cldco2(i)
202          dr_cld(i)  = rb_cldco2(i+1) - rb_cldco2(i)
203        enddo
204        rb_cldco2(nbinco2_cld+1) = rbmax_cld
205        dr_cld(nbinco2_cld)   = rb_cldco2(nbinco2_cld+1) -
206     &       rb_cldco2(nbinco2_cld)
207        print*, ' '
208        print*,'Microphysics co2: size bin information:'
209        print*,'i,rb_cldco2(i), rad_cldco2(i),dr_cld(i)'
210        print*,'-----------------------------------'
211        do i=1,nbinco2_cld
[1884]212          write(*,'(i3,3x,3(e13.6,4x))') i,rb_cldco2(i), rad_cldco2(i),
[1617]213     &      dr_cld(i)
214        enddo
[1884]215        write(*,'(i3,3x,e13.6)') nbinco2_cld+1,rb_cldco2(nbinco2_cld+1)
[1617]216        print*,'-----------------------------------'
217        do i=1,nbinco2_cld+1
[1816]218            rb_cldco2(i) = dlog(rb_cldco2(i))  !! we save that so that it is not computed
[1617]219                                         !! at each timestep and gridpoint
220        enddo
221c       Contact parameter of co2 ice on dst ( m=cos(theta) )
222c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
223c        mteta  = 0.952
224c        mtetaco2 = 0.952
[1921]225c        write(*,*) 'co2_param contact parameter:', mtetaco2
226
[1617]227c       Volume of a co2 molecule (m3)
[1921]228        vo1co2 = m0co2 / dble(rho_ice_co2) ! m0co2 et non mco2
229
[1617]230c       Variance of the ice and CCN distributions
231        sigma_iceco2 = sqrt(log(1.+nuiceco2_sed))
[1720]232        sigma_ice = sqrt(log(1.+nuice_sed))
233
[1685]234 
235        write(*,*) 'Variance of ice & CCN CO2 distribs :', sigma_iceco2
236        write(*,*) 'nuice for co2 ice sedimentation:', nuiceco2_sed
237        write(*,*) 'Volume of a co2 molecule:', vo1co2
[1720]238       
239        if (co2useh2o) then
240           write(*,*)
241           write(*,*) "co2useh2o=.true. in callphys.def"
242           write(*,*) "This means water ice particles can "
243           write(*,*) "serve as CCN for CO2 microphysics"
244        endif
[1921]245
[1720]246        if (meteo_flux) then
247           write(*,*)
248           write(*,*) "meteo_flux=.true. in callphys.def"
249           write(*,*) "meteoritic dust particles are available"
250           write(*,*) "for co2 ice nucleation! "
[1816]251           write(*,*) "Flux given by J. Plane (pressions,size bins)"
[1720]252           ! Initialisation of the flux: it is constant and is it saved
[1816]253           !We must interpolate the table to the GCM pressures
[1918]254           INQUIRE(FILE=TRIM(datadir)//
[1720]255     &       '/Meteo_flux_Plane.dat', EXIST=file_ok)
256           IF (.not. file_ok) THEN
257              write(*,*) 'file Meteo_flux_Plane.dat should be in '
[1918]258     &             ,trim(datadir)
[1720]259              STOP
260           endif
261!used Variables
[1918]262           open(unit=uMeteor,file=trim(datadir)//
[1720]263     &          '/Meteo_flux_Plane.dat'
264     &          ,FORM='formatted')
[1816]265!13000 records (130 pressions x 100 bin sizes)
[1884]266           read(uMeteor,*) !skip 1 line
[1720]267           do i=1,130
[1884]268              read(uMeteor,*,iostat=read_ok) pression_meteor(i)
269              if (read_ok==0) then
270                write(*,*) pression_meteor(i)
271              else
272                write(*,*) 'Error reading Meteo_flux_Plane.dat'
273                call abort_physic("CO2clouds",
274     &                    "Error reading Meteo_flux_Plane.dat",1)
275              endif
[1816]276           enddo
[1884]277           read(uMeteor,*)       !skip 1 line
[1816]278           do i=1,130
[1720]279              do j=1,100        ! les mêmes 100 bins size que la distri nuclea: on touche pas
[1884]280                 read(uMeteor,'(F12.6)',iostat=read_ok) meteor(i,j)
281                 if (read_ok/=0) then
282                   write(*,*) 'Error reading Meteo_flux_Plane.dat'
283                   call abort_physic("CO2clouds",
284     &                    "Error reading Meteo_flux_Plane.dat",1)
285                 endif
[1720]286              enddo
[1816]287!On doit maintenant réinterpoler le tableau(130,100) sur les pressions du GCM (nlay,100)
[1720]288           enddo
[1884]289           close(uMeteor)
[1816]290        write(*,*) "File meteo_flux read, end of firstcall in co2 micro"
291        endif                     !of if meteo_flux
292     
293        ! Parameter values for Latent heat computation
[1720]294        l0=595594d0     
295        l1=903.111d0   
296        l2=-11.5959d0   
297        l3=0.0528288d0
298        l4=-0.000103183d0
[1617]299        firstcall=.false.
300      END IF
301!=============================================================
302! 1. Initialisation
303!=============================================================
[1921]304
[1884]305      meteor_ccn(:,:,:)=0.
[1685]306      rice(:,:) = 1.e-8
[1617]307      riceco2(:,:) = 1.e-11
308
309c     Initialize the tendencies
[1911]310      subpdqcloudco2(1:ngrid,1:nlay,1:nq)=0.
311      subpdtcloudco2(1:ngrid,1:nlay)=0.
[1617]312     
[1911]313c pteff temperature layer; sum_subpdt dT.s-1
314c pqeff traceur kg/kg; sum_subpdq tendance idem .s-1
[1617]315      zt(1:ngrid,1:nlay) =
[1911]316     &      pteff(1:ngrid,1:nlay) +
317     &      sum_subpdt(1:ngrid,1:nlay) * microtimestep
[1617]318      zq(1:ngrid,1:nlay,1:nq) =
[1911]319     &      pqeff(1:ngrid,1:nlay,1:nq) +
320     &      sum_subpdq(1:ngrid,1:nlay,1:nq) * microtimestep
[1617]321      WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
[1816]322     &     zq(1:ngrid,1:nlay,1:nq) = 1.e-30
323         
324         zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
[1617]325!=============================================================
326! 2. Compute saturation
327!=============================================================
328      dev2 = 1. / ( sqrt(2.) * sigma_iceco2 )
[1720]329      dev3 = 1. / ( sqrt(2.) * sigma_ice )
[1617]330      call co2sat(ngrid*nlay,zt,pplay,zqsat) !zqsat is psat(co2)
[1816]331      zqco2=zq(:,:,igcm_co2)+zq(:,:,igcm_co2_ice)
332      CALL tcondco2(ngrid,nlay,pplay,zqco2,tcond)
[1720]333!=============================================================
[1816]334! 3. Bonus: additional meteoritic particles for nucleation
[1720]335!=============================================================
336      if (meteo_flux) then
[1816]337         !pression_meteo(130)
338         !pplev(ngrid,nlay+1)
[1720]339         !meteo(130,100)
340         !resultat: meteo_ccn(ngrid,nlay,100)
341         do l=1,nlay
342            do ig=1,ngrid
[1816]343               masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
[1884]344               ltemp1=abs(pression_meteor(:)-pplev(ig,l))
345               ltemp2=abs(pression_meteor(:)-pplev(ig,l+1))
[1720]346               lebon1=minloc(ltemp1,DIM=1)
347               lebon2=minloc(ltemp2,DIM=1)
348               nelem=lebon2-lebon1+1.
349               mtemp(:)=0d0     !mtemp(100) : valeurs pour les 100bins
350               do ibin=1,100
[1884]351                  mtemp(ibin)=sum(meteor(lebon1:lebon2,ibin))
[1720]352               enddo
[1884]353               meteor_ccn(ig,l,:)=mtemp(:)/nelem/masse(ig,l) !Par kg air
[1816]354csi par m carre, x epaisseur/masse pour par kg/air
355               !write(*,*) "masse air ig l=",masse(ig,l)
356               !check original unit with J. Plane
[1720]357            enddo
358         enddo
359      endif
[1816]360c ------------------------------------------------------------------------
361c ---------  Actual microphysics : Main loop over the GCM's grid ---------
362c ------------------------------------------------------------------------
[1685]363       DO l=1,nlay
364         DO ig=1,ngrid
[1617]365c       Get the partial pressure of co2 vapor and its saturation ratio
366           pco2 = zq(ig,l,igcm_co2) * (mmean(ig,l)/44.01) * pplay(ig,l)
367           satu = pco2 / zqsat(ig,l)
[1816]368
[1685]369           rho_ice_co2T(ig,l)=1000.*(1.72391-2.53e-4*zt(ig,l)
[1816]370     &          -2.87e-6*zt(ig,l)*zt(ig,l)) !T-dependant CO2 ice density
[1720]371           vo2co2 = m0co2 / dble(rho_ice_co2T(ig,l))
[1685]372           rho_ice_co2=rho_ice_co2T(ig,l)
[1720]373
[1816]374!=============================================================
375!4. Nucleation
376!=============================================================
377           IF ( satu .ge. 1 ) THEN ! if there is condensation
[1617]378
[1921]379c We do rdust computation "by hand" because we don't want to
380c change the mininumum rdust value in updaterad... It would
381c mess up various part of the GCM !
[1816]382             
[1617]383              rdust(ig,l)= zq(ig,l,igcm_dust_mass)
[1921]384     &             *0.75/(pi*rho_dust
385     &             * zq(ig,l,igcm_dust_number))
[1617]386              rdust(ig,l)= rdust(ig,l)**(1./3.)
[1816]387              if (zq(ig,l,igcm_dust_mass)*tauscaling(ig) .le. 1.e-20
388     &             .or.  zq(ig,l,igcm_dust_number)*tauscaling(ig) .le.1
389     &             .or. rdust(ig,l) .le. 1.e-9) then
390                 rdust(ig,l)=1.e-9
391              endif
392              rdust(ig,l)=min(5.e-4,rdust(ig,l))
[1685]393     
[1921]394c Expand the dust moments into a binned distribution
395             
396              n_aer(:)=0.
397              m_aer(:)=0.
398
[1685]399              Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)+1.e-30
400              No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+1.e-30
[1921]401
402              No_dust=No
403              Mo_dust=Mo
404
[1617]405              Rn = rdust(ig,l)
[1816]406              Rn = -dlog(Rn)
[1617]407              Rm = Rn - 3. * sigma_iceco2*sigma_iceco2 
[1816]408              n_derf = derf( (rb_cldco2(1)+Rn) *dev2)
409              m_derf = derf( (rb_cldco2(1)+Rm) *dev2)
[1921]410
[1617]411              do i = 1, nbinco2_cld
412                 n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
413                 m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
414                 n_derf = derf((rb_cldco2(i+1)+Rn) *dev2)
415                 m_derf = derf((rb_cldco2(i+1)+Rm) *dev2)
[1921]416                 n_aer(i) = n_aer(i) + 0.5 * No * n_derf
[1617]417                 m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
[1921]418              enddo
419
420c Ajout meteor_ccn particles aux particules de poussière background
421              if (meteo_flux) then
422                 do i = 1, nbinco2_cld
423                    n_aer(i) = n_aer(i) + meteor_ccn(ig,l,i)
424                    m_aer(i) = m_aer(i) + 4./3.*pi*rho_dust
[1884]425     &                *meteor_ccn(ig,l,i)*rad_cldco2(i)*rad_cldco2(i)
[1720]426     &                *rad_cldco2(i)
[1921]427                 enddo
428              endif
[1617]429             
[1921]430c Same but with h2o particles as CCN only if co2useh2o=.true.
431             
432              n_aer_h2oice(:)=0.
433              m_aer_h2oice(:)=0.
434
435              if (co2useh2o) then
436                call updaterice_micro(zq(ig,l,igcm_h2o_ice),
437     &               zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
438     &                 tauscaling(ig),rice(ig,l),rhocloud(ig,l))
439                Mo = zq(ig,l,igcm_h2o_ice) +
440     &               zq(ig,l,igcm_ccn_mass)*tauscaling(ig)+1.e-30
441                     ! Total mass of H20 crystals,CCN included
442                No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
443                Rn = rice(ig,l)
444                Rn = -dlog(Rn)
445                Rm = Rn - 3. * sigma_ice*sigma_ice 
446                n_derf = derf( (rb_cldco2(1)+Rn) *dev3)
447                m_derf = derf( (rb_cldco2(1)+Rm) *dev3)
448                do i = 1, nbinco2_cld
449                  n_aer_h2oice(i) = -0.5 * No * n_derf
450                  m_aer_h2oice(i) = -0.5 * Mo * m_derf
451                  n_derf = derf( (rb_cldco2(i+1)+Rn) *dev3)
452                  m_derf = derf( (rb_cldco2(i+1)+Rm) *dev3)
453                  n_aer_h2oice(i) = n_aer_h2oice(i) + 0.5 * No * n_derf
454                  m_aer_h2oice(i) = m_aer_h2oice(i) + 0.5 * Mo * m_derf
455                  rad_h2oice(i) = rad_cldco2(i)
456                enddo
457              endif
458             
459
460! Call to nucleation routine
[1617]461              call nucleaCO2(dble(pco2),zt(ig,l),dble(satu)
462     &             ,n_aer,rate,n_aer_h2oice
[1685]463     &             ,rad_h2oice,rateh2o,vo2co2)
[1617]464              dN = 0.
465              dM = 0.
466              dNh2o = 0.
467              dMh2o = 0.
468              do i = 1, nbinco2_cld
[1921]469                 Proba    =1.0-dexp(-1.*microtimestep*rate(i))
[1617]470                 dN       = dN + n_aer(i) * Proba
[1816]471                 dM       = dM + m_aer(i) * Proba             
[1617]472              enddo
[1921]473              if (co2useh2o) then
474                 do i = 1, nbinco2_cld
475                    Probah2o = 1.0-dexp(-1.*microtimestep*rateh2o(i))
476                    dNh2o    = dNh2o + n_aer_h2oice(i) * Probah2o
477                    dMh2o    = dMh2o + m_aer_h2oice(i) * Probah2o
478                 enddo
479              endif
480
481! dM  masse activée (kg) et dN nb particules par  kg d'air
482! Now increment CCN tracers and update dust tracers
[1617]483              dNN= dN/tauscaling(ig)
484              dMM= dM/tauscaling(ig)
[1720]485              dNN=min(dNN,zq(ig,l,igcm_dust_number))
486              dMM=min(dMM,zq(ig,l,igcm_dust_mass))
[1617]487              zq(ig,l,igcm_ccnco2_mass)   =
488     &             zq(ig,l,igcm_ccnco2_mass)   + dMM
[1720]489              zq(ig,l,igcm_ccnco2_number) =
[1617]490     &             zq(ig,l,igcm_ccnco2_number) + dNN
[1720]491              zq(ig,l,igcm_dust_mass)= zq(ig,l,igcm_dust_mass)-dMM
492              zq(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number)-dNN
493
[1921]494
495c Update CCN for CO2 nucleating on H2O CCN :
[1685]496              ! Warning: must keep memory of it
[1816]497              if (co2useh2o) then
498                 dNNh2o=dNh2o/tauscaling(ig)
499                 dNNh2o=min(dNNh2o,zq(ig,l,igcm_ccn_number))
500                 ratioh2o_ccn=1./(zq(ig,l,igcm_h2o_ice)
501     &                +zq(ig,l,igcm_ccn_mass)*tauscaling(ig)) 
502                 dMh2o_ice=dMh2o*zq(ig,l,igcm_h2o_ice)*ratioh2o_ccn
503                 dMh2o_ccn=dMh2o*zq(ig,l,igcm_ccn_mass)*
504     &                tauscaling(ig)*ratioh2o_ccn
505                 dMh2o_ccn=dMh2o_ccn/tauscaling(ig)
506                 dMh2o_ccn=min(dMh2o_ccn,zq(ig,l,igcm_ccn_mass))
507                 dMh2o_ice=min(dMh2o_ice,zq(ig,l,igcm_h2o_ice))
508                 zq(ig,l,igcm_ccnco2_mass)   =
509     &                zq(ig,l,igcm_ccnco2_mass)  + dMh2o_ice+dMh2o_ccn
510                 zq(ig,l,igcm_ccnco2_number) =
511     &                zq(ig,l,igcm_ccnco2_number) + dNNh2o
512                zq(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)-dNNh2o
513                zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)-dMh2o_ice
514                zq(ig,l,igcm_ccn_mass)= zq(ig,l,igcm_ccn_mass)-dMh2o_ccn
[1922]515                mem_Mh2o_co2(ig,l)=mem_Mh2o_co2(ig,l)+dMh2o_ice
516                mem_Mccn_co2(ig,l)=mem_Mccn_co2(ig,l)+dMh2o_ccn
517                mem_Nccn_co2(ig,l)=mem_Nccn_co2(ig,l)+dNNh2o
[1921]518             endif ! of if co2useh2o
519           ENDIF   ! of is satu >1
520
[1617]521!=============================================================
[1921]522! 5. Ice growth: scheme for radius evolution
[1617]523!=============================================================
524
525c We trigger crystal growth if and only if there is at least one nuclei (N>1).
526c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
527c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
[1816]528           No   = zq(ig,l,igcm_ccnco2_number)* tauscaling(ig)+1.e-30
529           IF (No .ge. 1)THEN   ! we trigger crystal growth
[1930]530              interm1 = DBLE(zq(ig,l,igcm_co2_ice))
531              interm2 = DBLE(zq(ig,l,igcm_ccnco2_mass))
532              interm3 = DBLE(zq(ig,l,igcm_ccnco2_number))
533              call updaterice_microco2(interm1,
534     &            interm2,interm3,
[1816]535     &            tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
[1617]536
[1816]537              Ic_rice=0.
538              lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)**2 +
539     &             l3 * zt(ig,l)**3 + l4 * zt(ig,l)**4 !J.kg-1
540              facteurmax=abs(Tcond(ig,l)-zt(ig,l))*cpp/lw
541              !specific heat of co2 ice = 1000 J.kg-1.K-1
542              !specific heat of atm cpp = 744.5 J.kg-1.K-1
[1720]543
[1921]544c call scheme of microphys. mass growth for CO2
[1617]545              call massflowrateCO2(pplay(ig,l),zt(ig,l),
[1816]546     &             satu,riceco2(ig,l),mmean(ig,l),Ic_rice)
[1921]547c Ic_rice Mass transfer rate (kg/s) for a rice particle >0 si croissance !
[1816]548             
549              if (isnan(Ic_rice) .or. Ic_rice .eq. 0.) then
[1720]550                 Ic_rice=0.
[1911]551                 subpdtcloudco2(ig,l)=-sum_subpdt(ig,l)
[1816]552                 dMice=0
553                 
554              else
[1911]555                 dMice=zq(ig,l,igcm_ccnco2_number)*Ic_rice*microtimestep
[1816]556     &                *tauscaling(ig) ! Kg par kg d'air, >0 si croissance !
557                 !kg.s-1 par particule * nb particule par kg air*s
558                 ! = kg par kg air
559             
560              dMice = max(dMice,max(-facteurmax,-zq(ig,l,igcm_co2_ice)))
561              dMice = min(dMice,min(facteurmax,zq(ig,l,igcm_co2)))
[1921]562! facteurmax maximum quantity of CO2 that can sublime/condense according to available thermal energy
[1816]563! latent heat release       >0 if growth i.e. if dMice >0
[1911]564              subpdtcloudco2(ig,l)=dMice*lw/cpp/microtimestep
[1816]565! kgco2/kgair* J/kgco2 * 1/(J.kgair-1.K-1)/s= K par seconde
566              !Now update tracers
567              zq(ig,l,igcm_co2_ice) = zq(ig,l,igcm_co2_ice)+dMice
568              zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)-dMice
[1720]569              endif
[1617]570
571!=============================================================
[1921]572! 6. Dust cores releasing if no more co2 ice :
[1617]573!=============================================================
574
[1816]575              if (zq(ig,l,igcm_co2_ice).le. 1.e-25)THEN
[1921]576! On sublime tout
[1816]577                 if (co2useh2o) then
[1922]578                   if (mem_Mccn_co2(ig,l) .gt. 0) then
[1685]579                    zq(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass)
[1922]580     &                   +mem_Mccn_co2(ig,l)
[1921]581                   endif
[1922]582                   if (mem_Mh2o_co2(ig,l) .gt. 0) then
[1685]583                    zq(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
[1922]584     &                   +mem_Mh2o_co2(ig,l)
[1921]585                   endif
[1685]586                 
[1922]587                   if (mem_Nccn_co2(ig,l) .gt. 0) then
[1685]588                    zq(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)
[1922]589     &                   +mem_Nccn_co2(ig,l)
[1921]590                   endif
[1685]591                 endif
592                    zq(ig,l,igcm_dust_mass) =
593     &                   zq(ig,l,igcm_dust_mass)
594     &                   + zq(ig,l,igcm_ccnco2_mass)-
[1922]595     &                   (mem_Mh2o_co2(ig,l)+mem_Mccn_co2(ig,l))
[1685]596                    zq(ig,l,igcm_dust_number) =
597     &                   zq(ig,l,igcm_dust_number)
[1922]598     &                   + zq(ig,l,igcm_ccnco2_number)
599     &                   -mem_Nccn_co2(ig,l)
[1685]600                 
601                    zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)
602     &                   + zq(ig,l,igcm_co2_ice)
603                 
604                 zq(ig,l,igcm_ccnco2_mass)=0.
605                 zq(ig,l,igcm_co2_ice)=0.
606                 zq(ig,l,igcm_ccnco2_number)=0.
[1922]607                 mem_Nccn_co2(ig,l)=0.
608                 mem_Mh2o_co2(ig,l)=0.
609                 mem_Mccn_co2(ig,l)=0.
[1685]610                 riceco2(ig,l)=0.
[1720]611
[1816]612              endif !of if co2_ice <1e-25
613
[1921]614              ENDIF            ! of if NCCN > 1
[1720]615          ENDDO                ! of ig loop
[1921]616        ENDDO                  ! of nlayer loop 
[1816]617
618!=============================================================
[1921]619! 7. END: get cloud tendencies
[1816]620!=============================================================
621
622          ! Get cloud tendencies
[1911]623        subpdqcloudco2(1:ngrid,1:nlay,igcm_co2) =
[1816]624     &       (zq(1:ngrid,1:nlay,igcm_co2) -
[1911]625     &       zq0(1:ngrid,1:nlay,igcm_co2))/microtimestep
626        subpdqcloudco2(1:ngrid,1:nlay,igcm_co2_ice) =
[1816]627     &       (zq(1:ngrid,1:nlay,igcm_co2_ice) -
[1911]628     &       zq0(1:ngrid,1:nlay,igcm_co2_ice))/microtimestep
[1921]629       
630        if (co2useh2o) then
631           subpdqcloudco2(1:ngrid,1:nlay,igcm_h2o_ice) =
632     &         (zq(1:ngrid,1:nlay,igcm_h2o_ice) -
633     &         zq0(1:ngrid,1:nlay,igcm_h2o_ice))/microtimestep
634           subpdqcloudco2(1:ngrid,1:nlay,igcm_ccn_mass) =
635     &         (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
636     &         zq0(1:ngrid,1:nlay,igcm_ccn_mass))/microtimestep
637           subpdqcloudco2(1:ngrid,1:nlay,igcm_ccn_number) =
638     &         (zq(1:ngrid,1:nlay,igcm_ccn_number) -
639     &         zq0(1:ngrid,1:nlay,igcm_ccn_number))/microtimestep
640        endif
641
[1911]642        subpdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_mass) =
[1816]643     &       (zq(1:ngrid,1:nlay,igcm_ccnco2_mass) -
[1911]644     &       zq0(1:ngrid,1:nlay,igcm_ccnco2_mass))/microtimestep
645        subpdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_number) =
[1816]646     &       (zq(1:ngrid,1:nlay,igcm_ccnco2_number) -
[1911]647     &       zq0(1:ngrid,1:nlay,igcm_ccnco2_number))/microtimestep
648        subpdqcloudco2(1:ngrid,1:nlay,igcm_dust_mass) =
[1816]649     &       (zq(1:ngrid,1:nlay,igcm_dust_mass) -
[1911]650     &       zq0(1:ngrid,1:nlay,igcm_dust_mass))/microtimestep
651        subpdqcloudco2(1:ngrid,1:nlay,igcm_dust_number) =
[1816]652     &       (zq(1:ngrid,1:nlay,igcm_dust_number) -
[1911]653     &       zq0(1:ngrid,1:nlay,igcm_dust_number))/microtimestep
[1884]654
[1921]655
[1922]656c     TEST D.BARDET
657      call WRITEDIAGFI(ngrid,"No_dust","Nombre particules de poussière"
658     &        ,"part/kg",3,No_dust)
659      call WRITEDIAGFI(ngrid,"Mo_dust","Masse particules de poussière"
660     &        ,"kg/kg ",3,Mo_dust)     
661
662        END SUBROUTINE improvedCO2clouds
[1938]663
664        END MODULE improvedCO2clouds_mod
665
Note: See TracBrowser for help on using the repository browser.