source: trunk/LMDZ.MARS/libf/phymars/improvedco2clouds_mod.F90 @ 2541

Last change on this file since 2541 was 2494, checked in by cmathe, 4 years ago

Mars GCM:
co2_ice as scatterer in radiative transfert. Need co2clouds and

activeco2ice .eqv. True. Files involved:

  • aeropacity_mod.F
  • callradite_mod.F
  • physiq_mod.F
  • updatereffrad_mod.F
  • suaer.F90
  • determine co2_ice density from temperature. Used in riceco2 computation.

Files involved:

  • co2cloud.F90
  • improvedco2clouds_mod.F90
  • updaterad.F90
  • updatereffrad_mod.F
  • co2condens_mod4micro.F: variable initialization
  • initracer.F: add nuiceco2_ref = 0.2
  • phyredem.F: remove co2_ice from qsurf since co2_ice => co2ice
  • watercloud_mod.F: tiny typo

CM
  • Property svn:executable set to *
File size: 27.5 KB
Line 
1!======================================================================================================================!
2! Module: Scheme of co2 cloud formation ===============================================================================!
3!----------------------------------------------------------------------------------------------------------------------!
4! Authors: Joaquim Audouard, Constantino Listowski, Anni Määttänen
5! Date: 09/2016
6!----------------------------------------------------------------------------------------------------------------------!
7! Contains subroutines:
8!     - improvedco2clouds_mod: nucleation
9!======================================================================================================================!
10module improvedco2clouds_mod
11
12implicit none
13
14contains
15!======================================================================================================================!
16! SUBROUTINE: improvedco2clouds =======================================================================================!
17!======================================================================================================================!
18! Subject:
19!---------
20!  This routine is used to form CO2 clouds when a parcel of the GCM is saturated.
21!----------------------------------------------------------------------------------------------------------------------!
22! Comments:
23!----------
24!  Adaptation for CO2 clouds based on improvedclouds_mod.F
25!
26!  It includes the ability to have supersaturation, a computation of the nucleation rates, growthrates and the
27!  scavenging of dust particles by clouds. It is worth noting that the amount of dust is computed using the dust
28!  optical depth computed in aeropacity.F.
29!  That's why the variable called "tauscaling" is used to convert pq(dust_mass) and pq(dust_number), which are relative
30!  quantities, to absolute and realistic quantities stored in zq. This has to be done to convert the inputs into
31!  absolute values, but also to convert the outputs back into relative values which are then used by the sedimentation
32!  and advection schemes.
33!  CO2 ice particles can nucleate on both dust and water ice particles. When CO2 ice is deposited onto a water ice
34!  particles, the particle is removed from the water tracers. Memory of the origin of the co2 particles is kept and
35!  thus the water cycle shouldn't be modified by this.
36!  There is an energy limit to how much co2 can sublimate/condensate. It is defined by the difference of the GCM
37!  temperature with the co2 condensation temperature.
38!
39!                                                /!\ WARNING /!\
40!  No sedimentation of the water ice origin is performed in the microphysical timestep in co2cloud.F.
41!
42!  If meteoritic particles are activated and turn into co2 ice particles, then they will be reversed in the dust
43!  tracers if the cloud sublimates.
44!----------------------------------------------------------------------------------------------------------------------!
45! Paper:
46!-------
47!    see co2cloud.F90
48!----------------------------------------------------------------------------------------------------------------------!
49! Algorithm:
50!-----------
51!   0. Firstcall
52!     0.1. Bonus: meteoritic component, extract data
53!   1. Initialization
54!   2. Compute saturation
55!   3. Bonus: additional meteoritic particles for nucleation
56!   4. Actual microphysics: Main loop over the GCM's grid
57!     4.1 Nucleation
58!     4.2. Ice growth: scheme for radius evolution
59!     4.3 Dust cores releasing if no more co2 ice
60!   5. Get cloud tendencies
61!======================================================================================================================!
62  subroutine improvedCO2clouds(ngrid, nlay, microtimestep, pplay, pplev, pteff, sum_subpdt, pqeff, sum_subpdq, &
63                               subpdqcloudco2, subpdtcloudco2, nq, tauscaling, mem_Mccn_co2, mem_Mh2o_co2, &
64                               mem_Nccn_co2, rb_cldco2, sigma_iceco2, dev2)
65
66  use comcstfi_h,   only: pi, g, cpp
67
68  use updaterad,    only: updaterice_micro, updaterice_microco2, updaterccnCO2
69
70  use tracer_mod,   only: igcm_dust_mass, igcm_dust_number, rho_dust, igcm_h2o_ice, igcm_ccn_mass, igcm_ccn_number, &
71                          nuice_sed, igcm_co2, igcm_co2_ice, igcm_ccnco2_mass, igcm_ccnco2_number, nuiceco2_sed, &
72                          nuiceco2_ref
73
74  use conc_mod,     only: mmean
75
76  use datafile_mod, only: datadir
77
78  use density_co2_ice_mod, only: density_co2_ice
79
80  implicit none
81
82  include "callkeys.h"
83  include "microphys.h"
84!----------------------------------------------------------------------------------------------------------------------!
85! VARIABLES DECLARATION
86!----------------------------------------------------------------------------------------------------------------------!
87! Input arguments:
88!-----------------
89  integer, intent(in) :: &
90     nq,    &! number of tracers
91     ngrid, &! number of point grid
92     nlay    ! number of layer
93
94  real, intent(in) :: &
95     microtimestep,           &! physics time step (s)
96     pplay(ngrid,nlay),       &! mid-layer pressure (Pa)
97     pplev(ngrid,nlay+1),     &! inter-layer pressure (Pa)
98     pteff(ngrid,nlay),       &! temperature at the middle of the layers (K)
99     sum_subpdt(ngrid,nlay),  &! tendency on temperature from previous physical parametrizations
100     pqeff(ngrid,nlay,nq),    &! tracers (kg/kg)
101     tauscaling(ngrid),       &! convertion factor for qdust and Ndust
102     sum_subpdq(ngrid,nlay,nq) ! tendencies on tracers before condensation (kg/kg.s-1)
103
104  real,  intent(in) :: &
105     sigma_iceco2 ! Variance of the co2 ice and CCN distributions
106
107  double precision, intent(in) :: &
108     rb_cldco2(nbinco2_cld+1), & ! boundary values of each rad_cldco2
109     dev2
110!----------------------------------------------------------------------------------------------------------------------!
111! Output arguments:
112!------------------
113  real, intent(out) :: &
114     subpdtcloudco2(ngrid,nlay),  &! tendency on tracers due to CO2 condensation (K/s)
115     subpdqcloudco2(ngrid,nlay,nq) ! tendency on tracers due to CO2 condensation (kg/kg.s-1)
116!----------------------------------------------------------------------------------------------------------------------!
117! Local:
118!-------
119!----1) Parameters:
120!------------------
121  integer, parameter :: &
122! ---Meteoritic flux input file
123     nbin_meteor = 100, &! Dimension 1
124     nlev_meteor = 130, &! Dimension 2
125     uMeteor = 666,     &! File unit
126! ---Latent heat computation
127     l0 = 595594d0,     &! coeff from: Azreg-Aïnou (2005)
128     l1 = 903.111d0,    &!   Title: "Low-temperature data for carbon dioxide"
129     l2 = -11.5959d0,   &!   Pulication: eprint arXiv:1403.4403
130     l3 = 0.0528288d0,  &!
131     l4 = -0.000103183d0 !
132
133  real, parameter :: &
134     threshold = 1e-30 ! limit value
135
136  character(len=20), parameter:: &
137     file_meteoritic_flux = 'Meteo_flux_Plane.dat'
138!----------------------------------------------------------------------------------------------------------------------!
139!----2) Saved:
140!-------------
141  real, save :: &
142     sigma_ice  ! Variance of the h2o ice and CCN distributions
143
144  double precision, save :: &
145     meteor(nlev_meteor,nbin_meteor), &! Meteoritic flux read from file uMeteor
146     dev3                              ! 1. / ( sqrt(2.) * sigma_ice )
147
148  logical, save :: &
149     firstcall = .true. ! Used to compute saved variables
150!----------------------------------------------------------------------------------------------------------------------!
151!----3) Variables:
152!-----------------
153  integer :: &
154     ig,     &! loop on ngrid
155     l,      &! loop on nlay
156     i,      &! loop on nbinco2
157! ---Variables for meteoritic flux
158     ibin,   &! loop on nbin_meteor
159     nelem,  &! nb of points during interpolation of meteoritic flux
160     lebon1, &! index where P_meteor is the nearest of pplev(ig,l)
161     lebon2, &! index where P_meteor is the nearest of pplev(ig,l+1)
162     read_ok  ! file uMeteor iostat
163
164  real :: &
165    masse(ngrid,nlay),     &! mass layer (kg.m-2)
166    rice(ngrid,nlay),      &! water ice mass mean radius (m): used for nucleation of CO2 on ice-coated ccns
167    zq(ngrid,nlay,nq),     &! local value of tracers (kg/kg)
168    zq0(ngrid,nlay,nq),    &! local init value of tracers (kg/kg)
169    zt(ngrid,nlay),        &! local value of temperature (K)
170    zqsat(ngrid,nlay),     &! saturation vapor pressure for CO2 (K)
171    tcond(ngrid,nlay),     &! condensation temperature of CO2 (K)
172    lw,                    &! Latent heat of sublimation (J.kg-1)
173    rdust(ngrid,nlay),     &! Dust geometric mean radius (m)
174    rhocloud(ngrid,nlay),  &! Cloud density (kg.m-3)
175    rhocloudco2(ngrid,nlay) ! Cloud density (kg.m-3)
176
177  real(kind=8) :: &
178    derf ! Error function
179 
180  double precision :: &
181    dMice,                             &! mass of condensed ice
182    facteurmax,                        &! for energy limit on mass growth
183    pco2,                              &! Co2 vapor partial pressure (Pa)
184    satu,                              &! Co2 vapor saturation ratio over ice
185    Mo,                                &! mass of aerosol particles
186    No,                                &! number of aerosol particles
187    Rn,                                &! logarithm of rdust/rice
188    Rm,                                &! Rn * variance of ice and CCN distribution
189    n_derf,                            &! derf( (rb_cldco2(1)+Rn) *dev3)
190    m_derf,                            &! derf( (rb_cldco2(1)+Rm) *dev2)
191    mem_Mccn_co2(ngrid,nlay),          &! Memory of CCN mass of H2O and dust used by CO2
192    mem_Mh2o_co2(ngrid,nlay),          &! Memory of H2O mass integred into CO2 crystal
193    mem_Nccn_co2(ngrid,nlay),          &! Memory of CCN number of H2O and dust used by CO2
194    n_aer(nbinco2_cld),                &! Radius used by the microphysical scheme (m)
195    m_aer(nbinco2_cld),                &! number concentration V-1 of particle/each size bin
196    n_aer_h2oice(nbinco2_cld),         &! mass mixing ratio of particle/each size bin
197    m_aer_h2oice(nbinco2_cld),         &! Same - for CO2 nucleation
198    rad_h2oice(nbinco2_cld),           &! Same - for CO2 nucleation
199    Ic_rice,                           &! Mass transfer rate CO2 ice crystal
200    ratioh2o_ccn,                      &! 1./(zq(ig,l,igcm_h2o_ice)  + zq(ig,l,igcm_ccn_mass)*tauscaling(ig))
201    vo2co2,                            &! volume of co2 ice particle
202    dN,                                &! number of particle of dust used as ccn
203    dM,                                &! mass of dN
204    dNh2o,                             &! number of particle of h2o ice used as ccn
205    dMh2o,                             &! mass of dNh2o
206    dNN,                               &! min(dN,zq(ig,l,igcm_dust_number))
207    dMM,                               &! min(dM,zq(ig,l,igcm_dust_mass))
208    dNNh2o,                            &! min(dNNh2o,zq(ig,l,igcm_ccn_number))
209    dMh2o_ice,                         &! min(dMh2o*zq(ig,l,igcm_h2o_ice)*ratioh2o_ccn, zq(ig,l,igcm_h2o_ice))
210    dMh2o_ccn,                         &! min(dMh2o_ccn,zq(ig,l,igcm_ccn_mass))
211    rate(nbinco2_cld),                 &! nucleation rate
212    rateh2o(nbinco2_cld),              &! nucleation rate for h2o
213    rho_ice_co2T,                      &! density of co2 ice Temperature-dependent
214    riceco2(ngrid,nlay),               &! CO2 ice mean radius (m)
215    vrat_cld,                          &! Volume ratio
216    Proba,                             &! 1.0 - exp(-1.*microtimestep*rate(i))
217    Probah2o,                          &! 1.0 - exp(-1.*microtimestep*rateh2o(i))
218    mtemp(nbinco2_cld),                &! sum(meteor(lebon1:lebon2,ibin))
219    pression_meteor(nlev_meteor),      &! pressure from meteoritic flux input file
220    ltemp1(nlev_meteor),               &! abs(pression_meteor(:)-pplev(ig,l))
221    ltemp2(nlev_meteor),               &! abs(pression_meteor(:)-pplev(ig,l+1))
222    meteor_ccn(ngrid,nlay,nbinco2_cld)  ! nbinco2_cld = 100
223
224  logical :: &
225    file_ok ! test if meteoritic input file exists
226!======================================================================================================================!
227! BEGIN ===============================================================================================================!
228!======================================================================================================================!
229! 0. Firstcall
230!----------------------------------------------------------------------------------------------------------------------!
231  if (firstcall) then
232    firstcall = .false.
233
234!   Variance of the ice and CCN distributions
235    sigma_ice = sqrt(log(1.+nuice_sed))
236    dev3 = 1. / ( sqrt(2.) * sigma_ice )
237!----------------------------------------------------------------------------------------------------------------------!
238! 0.1. Bonus: meteoritic component, extract data
239!----------------------------------------------------------------------------------------------------------------------!
240    if (meteo_flux) then
241      ! Check if file exists
242      inquire(file=trim(datadir)//'/'//file_meteoritic_flux, exist=file_ok)
243      if (.not. file_ok) then
244        call abort_physic("CO2clouds", 'file '//file_meteoritic_flux//' should be in'//trim(datadir), 1)
245      end if
246
247      ! open file
248      open(unit=uMeteor,file=trim(datadir)//'/'//file_meteoritic_flux, FORM='formatted')
249
250      !skip 1 line
251      read(uMeteor,*)
252
253      ! extract pressure_meteor
254      do i = 1, nlev_meteor
255        read(uMeteor,*)pression_meteor(i)
256      end do
257
258      !skip 1 line
259      read(uMeteor,*)
260
261      ! extract meteor flux
262      do i = 1, nlev_meteor
263        ! les mêmes 100 bins size que la distri nuclea : on touche pas
264        do ibin = 1, nbin_meteor
265          read(uMeteor,'(F12.6)') meteor(i,ibin)
266        end do
267      end do
268
269      ! close file
270      close(uMeteor)
271    end if ! of if meteo_flux
272  end if ! firstcall
273!----------------------------------------------------------------------------------------------------------------------!
274! 1. Initialization
275!----------------------------------------------------------------------------------------------------------------------!
276  rdust(:,:) = 0.
277  meteor_ccn(:,:,:) = 0.
278  rice(:,:) = 1.e-8
279  riceco2(:,:) = 1.e-11
280
281  ! Initialize the tendencies
282  subpdqcloudco2(:,:,:) = 0.
283  subpdtcloudco2(:,:) = 0.
284
285  ! pteff temperature layer; sum_subpdt dT.s-1
286  zt(1:ngrid,1:nlay) = 0.
287  zt(:,:) = pteff(:,:) + sum_subpdt(:,:) * microtimestep
288
289  ! pqeff traceur kg/kg; sum_subpdq tendance idem .s-1
290  zq(:,:,:) = pqeff(:,:,:) + sum_subpdq(:,:,:) * microtimestep
291  where( zq(:,:,:) < threshold ) zq(:,:,:) = threshold
292
293  zq0(:,:,:) = zq(:,:,:)
294  zqsat(:,:) = 0.
295!----------------------------------------------------------------------------------------------------------------------!
296! 2. Compute saturation
297!----------------------------------------------------------------------------------------------------------------------!
298  call co2sat(ngrid*nlay,zt,zqsat)
299  call tcondco2(ngrid,nlay,pplay, zq(:,:,igcm_co2) + zq(:,:,igcm_co2_ice),tcond)
300!----------------------------------------------------------------------------------------------------------------------!
301! 3. Bonus: additional meteoritic particles for nucleation
302!----------------------------------------------------------------------------------------------------------------------!
303  ! TODO: instead of intepolation, used only the nearest pplev(ig,l)
304  if (meteo_flux) then
305    do l = 1, nlay
306      do ig = 1, ngrid
307        masse(ig,l) = (pplev(ig,l) - pplev(ig,l+1)) / g
308
309        ltemp1 = abs(pression_meteor(:)-pplev(ig,l))
310        ltemp2 = abs(pression_meteor(:)-pplev(ig,l+1))
311
312        lebon1 = minloc(ltemp1,DIM=1)
313        lebon2 = minloc(ltemp2,DIM=1)
314
315        nelem = lebon2-lebon1+1.
316
317        mtemp(:) = 0d0
318
319        do ibin = 1, nbin_meteor
320          mtemp(ibin) = sum(meteor(lebon1:lebon2,ibin))
321        end do
322
323        ! Par kg air csi par m carre, x epaisseur/masse pour par kg/air. Check original unit with J. Plane
324        meteor_ccn(ig,l,:)=mtemp(:)/nelem/masse(ig,l)
325      end do
326    end do
327  end if
328!----------------------------------------------------------------------------------------------------------------------!
329! 4. Actual microphysics: Main loop over the GCM's grid
330!----------------------------------------------------------------------------------------------------------------------!
331  do l = 1, nlay
332    do ig = 1, ngrid
333
334      ! Get the partial pressure of co2 vapor and its saturation ratio
335      pco2 = zq(ig,l,igcm_co2) * (mmean(ig,l)/(mco2*1e3)) * pplay(ig,l) ! mco2 (kg/mol) => mmean (g/mol)
336      satu = pco2 / zqsat(ig,l)
337
338      !T-dependant CO2 ice density
339      call density_co2_ice(dble(zt(ig,l)), rho_ice_co2T)
340
341      vo2co2 = m0co2 / rho_ice_co2T
342!----------------------------------------------------------------------------------------------------------------------!
343! 4.1 Nucleation
344!----------------------------------------------------------------------------------------------------------------------!
345      ! if there is condensation
346      if ( satu >= 1 ) then
347        call updaterccnCO2(zq(ig,l,igcm_dust_mass), zq(ig,l,igcm_dust_number), rdust(ig,l), tauscaling(ig))
348
349        ! Expand the dust moments into a binned distribution
350        n_aer(:) = 0d0 ! number of aerosol particles
351        m_aer(:) = 0d0 ! mass of aerosol particles
352
353        No = zq(ig,l,igcm_dust_number) * tauscaling(ig)
354
355        Mo = (4./3.) * pi * rho_dust * No * rdust(ig,l)**3 *exp(9.*nuiceco2_ref/2.)
356
357        if (No > threshold) then
358          Rn = -log(rdust(ig,l))
359
360          Rm = Rn - 3. * sigma_iceco2 * sigma_iceco2
361
362          n_derf = derf( (rb_cldco2(1)+Rn) *dev2)
363          m_derf = derf( (rb_cldco2(1)+Rm) *dev2)
364
365          do i = 1, nbinco2_cld
366            n_aer(i) = -0.5 * No * n_derf
367            m_aer(i) = -0.5 * Mo * m_derf
368
369            n_derf = derf((rb_cldco2(i+1)+Rn) *dev2)
370            m_derf = derf((rb_cldco2(i+1)+Rm) *dev2)
371
372            n_aer(i) = n_aer(i) + 0.5 * No * n_derf
373            m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
374          end do
375
376          ! Ajout meteor_ccn particles aux particules de poussière background
377          if (meteo_flux) then
378            do i = 1, nbinco2_cld
379              n_aer(i) = n_aer(i) + meteor_ccn(ig,l,i)
380
381              m_aer(i) = m_aer(i) + (4./3.) * pi * rho_dust *meteor_ccn(ig,l,i) * rad_cldco2(i)**3
382             end do
383          end if
384
385          ! Same but with h2o particles as CCN only if co2useh2o = .true.
386          n_aer_h2oice(:)=0.
387          m_aer_h2oice(:)=0.
388
389          if (co2useh2o) then
390            call updaterice_micro(zq(ig,l,igcm_h2o_ice), zq(ig,l,igcm_ccn_mass), zq(ig,l,igcm_ccn_number), &
391                                  tauscaling(ig), rice(ig,l), rhocloud(ig,l))
392
393            Mo = zq(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_ccn_mass) * tauscaling(ig) + threshold
394
395            ! Total mass of H20 crystals,CCN included
396            No = zq(ig,l,igcm_ccn_number) * tauscaling(ig) + threshold
397
398            Rn = -log(rice(ig,l))
399
400            Rm = Rn - 3. * sigma_ice * sigma_ice
401
402            n_derf = derf( (rb_cldco2(1)+Rn) *dev3)
403            m_derf = derf( (rb_cldco2(1)+Rm) *dev3)
404
405            do i = 1, nbinco2_cld
406              n_aer_h2oice(i) = -0.5 * No * n_derf
407              m_aer_h2oice(i) = -0.5 * Mo * m_derf
408
409              n_derf = derf( (rb_cldco2(i+1)+Rn) *dev3)
410              m_derf = derf( (rb_cldco2(i+1)+Rm) *dev3)
411
412              n_aer_h2oice(i) = n_aer_h2oice(i) + 0.5 * No * n_derf
413              m_aer_h2oice(i) = m_aer_h2oice(i) + 0.5 * Mo * m_derf
414
415              rad_h2oice(i) = rad_cldco2(i)
416            end do
417          end if
418
419          ! Call to nucleation routine
420          call nucleaCO2(dble(pco2), zt(ig,l), dble(satu), n_aer, rate, n_aer_h2oice, rad_h2oice, rateh2o, vo2co2)
421
422          dN = 0.
423          dM = 0.
424          dNh2o = 0.
425          dMh2o = 0.
426
427          do i = 1, nbinco2_cld
428            Proba = 1.0 - exp(-1.*microtimestep*rate(i))
429            dN = dN + n_aer(i) * Proba
430            dM = dM + m_aer(i) * Proba
431          end do
432
433          if (co2useh2o) then
434            do i = 1, nbinco2_cld
435              Probah2o = 1.0 - exp(-1.*microtimestep*rateh2o(i))
436              dNh2o = dNh2o + n_aer_h2oice(i) * Probah2o
437              dMh2o = dMh2o + m_aer_h2oice(i) * Probah2o
438            end do
439          end if
440
441          ! Now increment CCN tracers and update dust tracers
442          dNN = min(dN,zq(ig,l,igcm_dust_number)) ! dNN est devenu DN
443          dMM = min(dM,zq(ig,l,igcm_dust_mass))  ! idem dans le min
444
445          zq(ig,l,igcm_ccnco2_mass) = zq(ig,l,igcm_ccnco2_mass) + dMM /tauscaling(ig)
446
447          zq(ig,l,igcm_ccnco2_number) = zq(ig,l,igcm_ccnco2_number) + dNN /tauscaling(ig)
448
449          zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - dMM /tauscaling(ig)
450
451          zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - dNN /tauscaling(ig)
452
453          ! Update CCN for CO2 nucleating on H2O CCN : Warning: must keep memory of it
454          if (co2useh2o) then
455            dNNh2o = dNh2o/tauscaling(ig)
456            dNNh2o = min(dNNh2o,zq(ig,l,igcm_ccn_number))
457
458            ratioh2o_ccn = 1./(zq(ig,l,igcm_h2o_ice)  + zq(ig,l,igcm_ccn_mass)*tauscaling(ig))
459
460            dMh2o_ccn = dMh2o * zq(ig,l,igcm_ccn_mass) *  tauscaling(ig) * ratioh2o_ccn
461            dMh2o_ccn = dMh2o_ccn/tauscaling(ig)
462            dMh2o_ccn = min(dMh2o_ccn,zq(ig,l,igcm_ccn_mass))
463
464            dMh2o_ice = dMh2o * zq(ig,l,igcm_h2o_ice) * ratioh2o_ccn
465            dMh2o_ice = min(dMh2o_ice,zq(ig,l,igcm_h2o_ice))
466
467            zq(ig,l,igcm_ccnco2_mass) = zq(ig,l,igcm_ccnco2_mass)  + dMh2o_ice + dMh2o_ccn
468
469            zq(ig,l,igcm_ccnco2_number) =  zq(ig,l,igcm_ccnco2_number) + dNNh2o
470
471            zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number)  - dNNh2o
472
473            zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)  - dMh2o_ice
474
475            zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) - dMh2o_ccn
476
477            mem_Mh2o_co2(ig,l) = mem_Mh2o_co2(ig,l) + dMh2o_ice
478            mem_Mccn_co2(ig,l) = mem_Mccn_co2(ig,l) + dMh2o_ccn
479            mem_Nccn_co2(ig,l) = mem_Nccn_co2(ig,l) + dNNh2o
480          end if ! of if co2useh2o
481        end if ! of if No > 1e-30
482      end if ! of is satu > 1
483!----------------------------------------------------------------------------------------------------------------------!
484! 4.2. Ice growth: scheme for radius evolution
485!----------------------------------------------------------------------------------------------------------------------!
486!    We trigger crystal growth if and only if there is at least one nuclei (N>1). Indeed, if we are supersaturated
487!      and still don't have at least one nuclei, we should better wait to avoid unrealistic value for nuclei radius
488!      and so on for cases that remain negligible.
489!----------------------------------------------------------------------------------------------------------------------!
490      ! we trigger crystal growth
491      if (zq(ig,l,igcm_ccnco2_number) * tauscaling(ig) + threshold >= 1) then
492
493        call updaterice_microco2(dble(zq(ig,l,igcm_co2_ice)), dble(zq(ig,l,igcm_ccnco2_mass)), &
494                                 dble(zq(ig,l,igcm_ccnco2_number)), zt(ig,l), tauscaling(ig), riceco2(ig,l), &
495                                 rhocloudco2(ig,l))
496        Ic_rice = 0.
497
498        ! J.kg-1
499        lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)**2 + l3 * zt(ig,l)**3 + l4 * zt(ig,l)**4
500
501        facteurmax = abs(tcond(ig,l)-zt(ig,l)) * (cpp/lw)
502
503        ! call scheme of microphys. mass growth for CO2 (evaporation/condensation)
504        call massflowrateCO2(pplay(ig,l), zt(ig,l), satu, riceco2(ig,l), mmean(ig,l), Ic_rice)
505
506        ! Ic_rice Mass transfer rate (kg/s) for a rice particle > 0 si croissance !
507        if (isnan(Ic_rice) .or. Ic_rice == 0.) then
508          Ic_rice = 0.
509          subpdtcloudco2(ig,l) = -sum_subpdt(ig,l)
510          dMice = 0
511        else
512          ! Kg par kg d'air, >0 si croissance !
513          ! kg.s-1 par particule * nb particule par kg air*s = kg par kg air
514          dMice = zq(ig,l,igcm_ccnco2_number) * Ic_rice * microtimestep * tauscaling(ig)
515
516          ! facteurmax maximum quantity of CO2 that can sublime/condense according to available thermal energy
517          ! latent heat release > 0 if growth i.e. if dMice > 0
518          dMice = max(dMice,max(-facteurmax,-zq(ig,l,igcm_co2_ice)))
519          dMice = min(dMice,min(facteurmax,zq(ig,l,igcm_co2)))
520
521          ! kgco2/kgair* J/kgco2 * 1/(J.kgair-1.K-1)/s = K /s
522          subpdtcloudco2(ig,l) = dMice * lw / cpp / microtimestep
523
524          !Now update tracers
525          zq(ig,l,igcm_co2_ice) = zq(ig,l,igcm_co2_ice) + dMice
526          zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2) - dMice
527        end if
528      end if ! if zq(ccnco2_number) >= 1
529!----------------------------------------------------------------------------------------------------------------------!
530! 4.3 Dust cores releasing if no more co2 ice
531!----------------------------------------------------------------------------------------------------------------------!
532      ! On sublime tout
533      if ((zq(ig,l,igcm_co2_ice) <= threshold).or.(zq(ig,l,igcm_ccnco2_number) * tauscaling(ig) < 1.)) then
534
535        if (co2useh2o) then
536
537          if (mem_Mccn_co2(ig,l) > 0) then
538            zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + mem_Mccn_co2(ig,l)
539          end if
540
541          if (mem_Mh2o_co2(ig,l) > 0) then
542            zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice) + mem_Mh2o_co2(ig,l)
543          end if
544
545          if (mem_Nccn_co2(ig,l) > 0) then
546            zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + mem_Nccn_co2(ig,l)
547          end if
548
549        end if
550
551        zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + zq(ig,l,igcm_ccnco2_mass) - ( mem_Mh2o_co2(ig,l) + &
552                                     mem_Mccn_co2(ig,l) )
553
554        zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + zq(ig,l,igcm_ccnco2_number)  - mem_Nccn_co2(ig,l)
555
556        zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)  + zq(ig,l,igcm_co2_ice)
557
558        zq(ig,l,igcm_co2_ice) = 0.
559        zq(ig,l,igcm_ccnco2_mass) = 0.
560        zq(ig,l,igcm_ccnco2_number) = 0.
561        mem_Nccn_co2(ig,l) = 0.
562        mem_Mh2o_co2(ig,l) = 0.
563        mem_Mccn_co2(ig,l) = 0.
564        riceco2(ig,l) = 0.
565      end if !of if co2_ice < threshold or zq(ccnco2_number) < 1
566    end do ! of ig loop
567  end do ! of nlayer loop
568!----------------------------------------------------------------------------------------------------------------------!
569! 5. Get cloud tendencies
570!----------------------------------------------------------------------------------------------------------------------!
571  subpdqcloudco2(:,:,igcm_co2) =   ( zq(:,:,igcm_co2) - zq0(:,:,igcm_co2) ) / microtimestep
572
573  subpdqcloudco2(:,:,igcm_co2_ice) =  ( zq(:,:,igcm_co2_ice) - zq0(:,:,igcm_co2_ice) ) / microtimestep
574
575  subpdqcloudco2(:,:,igcm_ccnco2_mass) =  ( zq(:,:,igcm_ccnco2_mass) - zq0(:,:,igcm_ccnco2_mass) ) / microtimestep
576
577  subpdqcloudco2(:,:,igcm_ccnco2_number) = ( zq(:,:,igcm_ccnco2_number) - zq0(:,:,igcm_ccnco2_number))/microtimestep
578
579  subpdqcloudco2(:,:,igcm_dust_mass) = ( zq(:,:,igcm_dust_mass) - zq0(:,:,igcm_dust_mass) )  / microtimestep
580
581  subpdqcloudco2(:,:,igcm_dust_number) = ( zq(:,:,igcm_dust_number) - zq0(:,:,igcm_dust_number) )   / microtimestep
582
583  if (co2useh2o) then
584    subpdqcloudco2(:,:,igcm_h2o_ice) = ( zq(:,:,igcm_h2o_ice) - zq0(:,:,igcm_h2o_ice) )  / microtimestep
585
586    subpdqcloudco2(:,:,igcm_ccn_mass) = ( zq(:,:,igcm_ccn_mass) - zq0(:,:,igcm_ccn_mass) )    / microtimestep
587
588    subpdqcloudco2(:,:,igcm_ccn_number) = ( zq(:,:,igcm_ccn_number) - zq0(:,:,igcm_ccn_number) ) / microtimestep
589  end if
590!======================================================================================================================!
591! END =================================================================================================================!
592!======================================================================================================================!
593  end subroutine improvedCO2clouds
594
595end module improvedCO2clouds_mod
596
Note: See TracBrowser for help on using the repository browser.