source: trunk/LMDZ.MARS/libf/phymars/co2cloud.F90 @ 2350

Last change on this file since 2350 was 2350, checked in by aslmd, 5 years ago

CO2 cloud microphysics
  • Property svn:executable set to *
File size: 48.5 KB
Line 
1!======================================================================================================================!
2! Module: CO2 clouds formation ========================================================================================!
3!----------------------------------------------------------------------------------------------------------------------!
4! Authors: Joachim Audouard, Constantino Listowski, Anni Määttänen
5! Date: 09/2016
6!----------------------------------------------------------------------------------------------------------------------!
7! Contains subroutines:
8!      - co2cloud: of co2 cloud microphysics
9!
10!      - ini_co2cloud: (only used in phys_state_var_init_mod.F90)
11!
12!      - end_co2cloud: (only used in phys_state_var_init_mod.F90)
13!======================================================================================================================!
14module co2cloud_mod
15
16implicit none
17
18double precision, allocatable, save :: &
19   mem_Mccn_co2(:,:), &! Memory of CCN mass of H2O and dust used by CO2
20   mem_Mh2o_co2(:,:), &! Memory of H2O mass integred into CO2 crystal
21   mem_Nccn_co2(:,:)   ! Memory of CCN number of H2O and dust used by CO2
22
23contains
24!======================================================================================================================!
25! SUBROUTINE: co2cloud ================================================================================================!
26!======================================================================================================================!
27! Subject:
28!---------
29!   Main subroutine of co2 cloud microphysics
30!----------------------------------------------------------------------------------------------------------------------!
31! Comments:
32!----------
33!     - Adaptation of the water ice clouds scheme (with specific microphysics) of Montmessin, Navarro et al.
34!
35!     - Microphysics subroutine is improvedCO2clouds.F
36!
37!     - There is a time loop specific to cloud formation due to timescales smaller than the GCM integration timestep
38!
39!     - The microphysics time step is a fraction of the physical one
40!
41!     - The CO2 clouds tracers (co2_ice, ccn mass and concentration) are sedimented at each microtimestep. pdqs_sedco2
42!       keeps track of the CO2 flux at the surface
43!
44!     - The subgrid Temperature distribution is modulated (0 or 1) by Spiga et al. (GRL 2012)
45!
46!     - Saturation Index to account for GW propagation or dissipation upwards
47!
48!     - 4D and column opacities are computed using Qext values at 1µm
49!----------------------------------------------------------------------------------------------------------------------!
50! Papers:
51!--------
52!   "Near-pure vapor condensation in the Martian atmosphere: CO2 ice crystal growth", Listowski et al. (2013), JGRE
53!   "Modeling the microphysics of CO2 ice clouds within wave-induced cold pockets in the martian mesosphere", Listowski
54!     et al. (2014), Icarus
55!   "Global climate modeling of the Martian water cycle with improved microphysics and radiatively active water ice
56!     clouds", Navarro et al. (2014), JGRE
57!   "Martian GCM with complete CO2 clouds microphysics", Audouard et al. (2017), EPSC abstract
58!----------------------------------------------------------------------------------------------------------------------!
59! Algorithm:
60!-----------
61!   0. Firstcall
62!     0.1. Initialization of microtimestep from imicroco2
63!     0.2. Compute the radius grid of CO2 ice particles (rb_cldco2)
64!     0.3. Read file 'optprop_co2ice_1mic.dat' to extract optical properties of CO2 ice at 1 micron (Qext)
65!     0.4. Interpole the radius grid (rb_cldco2) to get the corresponding exctinction coefficients (Qext1bins)
66!     0.5. Save the radius grid of CO2 particles (rb_cldco2)
67!   1. Initialization
68!   2. Compute mass and thickness layers
69!   3. Define the sub-grid cloud (CLFvaryingCO2)
70!     3.1. Representation of sub-grid CO2 ice clouds (CLFvaryingCO2 = True)
71!       3.1.a. Saturation index CO2
72!       3.1.b. Compute tcond
73!       3.1.c. Compute cloud fraction in cells
74!     3.2. No sub-grid cloud representation (CLFvarying = False)
75!   4. Microphysics of CO2 cloud formation
76!     4.1. Stepped entry for tendancies: At each micro timestep we add pdt in order to have a stepped entry
77!     4.2. Effective tracers quantities in the cloud
78!     4.3. Gravitational sedimentation
79!       4.3.a. Compute cloud density
80!       4.3.b. Save actualized tracer values to compute sedimentation tendancies
81!       4.3.c. Sedimentation of co2 ice
82!       4.3.d. Sedimentation for other tracers
83!       4.3.e. Compute tendencies due to the sedimation process
84!     4.4. Main call to the cloud scheme
85!     4.5. Updating tendencies after cloud scheme
86!   5. Compute final tendencies after time loop
87!   6. Update clouds physical values in the cloud (for output)
88!     6.1. Update density of co2 ice, riceco2 and opacity
89!     6.2. Update rice and rdust
90!   7. Correction if a lot of subliming CO2 fills the 1st layer
91!   8. Compute water cloud sedimentation radius
92!   9. CO2 saturated quantities
93!     9.1 Compute CO2 saturation in layers
94!     9.2 Compute CO2 saturated quantities in layers
95!   10. Everything modified by CO2 microphysics must be wrt co2cloudfrac
96!   11. Compute opacity at 1 micron
97!   12. Write outputs in diagfi.nc
98!======================================================================================================================!
99  subroutine co2cloud(ngrid, nlay, ptimestep, pplev, pplay, pdpsrf, pzlay, pt, pdt, pq, pdq, pdqcloudco2, pdtcloudco2, &
100                      nq, tau, tauscaling, rdust, rice, riceco2, nuice, rsedcloud, rhocloud, pzlev, pdqs_sedco2, pdu, &
101                      pu, pcondicea, co2ice)
102
103  use ioipsl_getincom, only: getin
104  use dimradmars_mod,  only: naerkind
105  use comcstfi_h,      only: pi, g, cpp
106  use updaterad,       only: updaterice_microco2, updaterice_micro, updaterdust
107  use conc_mod,        only: mmean, rnew
108  use tracer_mod,      only: igcm_co2, igcm_co2_ice, igcm_dust_mass, igcm_dust_number, igcm_h2o_ice, &
109                             igcm_ccn_mass, igcm_ccn_number, igcm_ccnco2_mass, igcm_ccnco2_number, rho_dust, &
110                             nuiceco2_sed, nuiceco2_ref, rho_ice_co2, r3n_q, rho_ice, nuice_sed
111
112  use newsedim_mod,    only: newsedim
113
114  use datafile_mod,    only: datadir
115
116  use improvedCO2clouds_mod, only: improvedCO2clouds, density_co2_ice
117
118  implicit none
119
120  include "callkeys.h"
121  include "microphys.h"
122!----------------------------------------------------------------------------------------------------------------------!
123! VARIABLES DECLARATION
124!----------------------------------------------------------------------------------------------------------------------!
125! Input arguments:
126!-----------------
127  integer, intent(in) ::&
128     ngrid, &! Number of grid points
129     nlay,  &! Number of layers
130     nq      ! Number of tracers
131
132  real, intent(in) :: &
133     ptimestep,           &! Physical time step (s)
134     pplev(ngrid,nlay+1), &! Inter-layer pressures (Pa)
135     pplay(ngrid,nlay),   &! Mid-layer pressures (Pa)
136     pdpsrf(ngrid),       &! Tendency on surface pressure
137     pzlay(ngrid,nlay),   &! Altitude at the middle of the layers
138     pt(ngrid,nlay),      &! Temperature at the middle of the layers (K)
139     pdt(ngrid,nlay),     &! Tendency on temperature from other parametrizations
140     pq(ngrid,nlay,nq),   &! Tracers (kg/kg)
141     pdq(ngrid,nlay,nq),  &! Tendencies before condensation (kg/kg.s-1)
142     tau(ngrid,naerkind), &! Column dust optical depth at each point
143     tauscaling(ngrid),   &! Convertion factor for dust amount
144     pu(ngrid,nlay),      &! Zonal Wind: zu = pu + (pdu * ptimestep)
145     pdu(ngrid,nlay),     &! Tendency of zonal wind before condensation
146     pzlev(ngrid,nlay+1), &! Altitude at the boundaries of the layers
147     nuice(ngrid,nlay),   &! Estimated effective variance of the size distribution
148     co2ice(ngrid)         ! Amount of co2 ice at the surface
149!----------------------------------------------------------------------------------------------------------------------!
150! Output arguments:
151!------------------
152  real, intent(out) :: &
153     rdust(ngrid,nlay),         & ! Dust geometric mean radius (m)
154     rice(ngrid,nlay),          & ! Water Ice mass mean radius (m)
155     rsedcloud(ngrid,nlay),     & ! Water Cloud sedimentation radius
156     rhocloud(ngrid,nlay),      & ! Water Cloud density (kg.m-3)
157     pdqs_sedco2(ngrid),        & ! CO2 flux at the surface
158     pdqcloudco2(ngrid,nlay,nq),& ! Tendency due to CO2 condensation (kg/kg.s-1)
159     pcondicea(ngrid,nlay),     & ! Rate of condensation/sublimation of co2 ice in layers
160     pdtcloudco2(ngrid,nlay),   & ! Tendency on temperature due to latent heat
161     riceco2(ngrid,nlay)          ! Ice mass mean radius (m) r_c in Montmessin et al. (2004)
162!----------------------------------------------------------------------------------------------------------------------!
163! Local:
164!-------
165!-----1) Parameters:
166!-------------------
167  integer, parameter :: &
168     uQext = 555,         &! file_qext unit ID
169     var_dim_qext = 10000  ! Exact dimension of radv and qextv1mic from file_qext
170
171  real, parameter :: &
172     mincloud = 0.1,  &! Minimum cloud fraction
173     beta = 0.85,     &! correction for the shape of the particles (see Murphy et al. JGR 1990 vol.95):
174                       !   beta = 1    for spheres
175                       !   beta = 0.85 for irregular particles
176                       !   beta = 0.5  for disk shaped particles
177     threshold = 1e-30 ! limit value
178
179  double precision, parameter :: &
180     rmin_cld = 1.e-9,  &! Minimum radius (m)
181     rmax_cld = 5.e-6,  &! Maximum radius (m)
182     rbmin_cld = 1.e-10,&! Minimum boundary radius (m)
183     rbmax_cld = 2.e-4, &! Maximum boundary radius (m)
184     Fo = 7.5e-7,       &! for sat index  (J.m-3)
185     lambdaH = 150.e3    ! for sat index (km)
186
187  character(len=23), parameter :: &
188     file_qext = 'optprop_co2ice_1mic.dat' ! File extinction coefficients of CO2 particles
189!----------------------------------------------------------------------------------------------------------------------!
190!-----2) Saved:
191!--------------
192  integer, save :: &
193     imicroco2 ! Time subsampling for coupled water microphysics sedimentation microtimestep timeloop for microphysics:
194               !   if imicroco2 = 1, subpdt is the same as pdt
195  real, save :: &
196     sigma_iceco2, &! Variance of the ice and CCN distributions
197     microtimestep  ! Integration timestep for coupled water microphysics & sedimentation
198
199  double precision, save :: &
200     dev2,                   &! 1. / ( sqrt(2.) * sigma_iceco2 )
201     Qext1bins(nbinco2_cld), &! Extinction coefficients for rb_cldco2 radius of CO2 ice particles
202     rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m)
203
204  logical, save :: &
205     firstcall = .true. ! Used to compute saved variables
206!----------------------------------------------------------------------------------------------------------------------!
207!-----3) Variables:
208!------------------
209  integer :: &
210     iq,       &! loop on tracers
211     ig,       &! loop on grid points
212     l,        &! loop on layers
213     i,        &! loop on nbinco2_cld
214     nelem,    &! number of point between lebon1 and lebon2 => interpolation
215     lebon1,   &! bound limit for the interpolation
216     lebon2,   &! bound limit for the interpolation
217     microstep  ! Time subsampling step variable
218
219  real :: &
220! ---Tendency given by clouds inside the micro loop
221     subpdqcloudco2(ngrid,nlay,nq), &! On tracers, cf. pdqcloud
222     subpdtcloudco2(ngrid,nlay),    &! On temperature, cf. pdtcloud
223!  ---Global tendency (clouds+physics)
224     sum_subpdq(ngrid,nlay,nq),     &! On tracers, cf. pdqcloud
225     sum_subpdt(ngrid,nlay),        &! On temperature, cf. pdtcloud
226!  ---Sedimentation
227     rsedcloudco2(ngrid,nlay),      &! Cloud sedimentation radius
228     ztsed(ngrid,nlay),             &! Temperature with real-time value in microtimeloop
229     zqsed(ngrid,nlay,nq),          &! Tracers with real-time value in µloop
230     zqsed0(ngrid,nlay,nq),         &! For sedimentation tendancy
231     subpdqsed(ngrid,nlay,nq),      &! Tendancy due to sedimentation
232     sum_subpdqs_sedco2(ngrid),     &! CO2 flux at the surface
233!  ---For sub grid T distribution
234     zt(ngrid,nlay),                &! Local value of temperature
235     zq_co2vap(ngrid, nlay),        &! Local value of CO2 vap
236     rhocloudco2t(ngrid, nlay),     &! Cloud density (kg.m-3)
237!  ---For Saturation Index computation
238     zdelt,                         &! Delta T for the temperature distribution
239     co2cloudfrac(ngrid,nlay),      &! Cloud fraction used only with CLFvarying is true
240!  ---Misc
241     rhocloudco2(ngrid, nlay),      &! Cloud density (kg.m-3)
242     Nccnco2,                       &! buffer: number of ccn used for co2 condensation
243     Qccnco2,                       &! buffer: mass of ccn used for co2 condensation
244     Niceco2,                       &! buffer: mmr co2 ice
245     epaisseur(ngrid,nlay),         &! Layer thickness (m)
246     masse(ngrid,nlay),             &! Layer  mass (kg.m-2)
247     pteff(ngrid, nlay),            &! Effective temperature in the cloud
248     pqeff(ngrid, nlay, nq),        &! Effective tracers quantities in the cloud
249     wq(ngrid,nlay+1),              &! Displaced tracer mass (kg.m-2) during microtimestep
250     satuco2(ngrid,nlay),           &! CO2 satu ratio for output diagfi
251     zqsatco2(ngrid,nlay)            ! Saturation co2
252
253  double precision :: &
254! ---Extinction coefficients at 1 micron of CO2 particles
255     vrat_cld,                      &! Volume ratio
256     n_derf,                        &! derf( (rb_cldco2(1)-dlog(riceco2(ig,l))) *dev2)
257     Qtemp,                         &! mean value in the interval during the interpolation
258     ltemp1(var_dim_qext),          &! abs(radv(:)-rb_cldco2(i))
259     ltemp2(var_dim_qext),          &! abs(radv(:)-rb_cldco2(i+1))
260     n_aer(nbinco2_cld),            &! -0.5 * Nccnco2*tauscaling(ig) * n_derf
261     tau1mic(ngrid),                &! CO2 ice column opacity at 1µm
262     Qext1bins2(ngrid,nlay),        &! CO2 ice opacities
263     radv(var_dim_qext),            &! radius of CO2 ice at 1 µm (read from file_qext)
264     Qextv1mic(var_dim_qext),       &! extinction coefficient of CO2 ice at 1 µm (read from file_qext)
265! ---For Saturation Index computation
266     rho,                           &! background density
267     zu,                            &! absolute value of zonal wind field
268     NN,                            &! N^2 static stability
269     gradT,                         &! thermal gradient
270     SatIndex(ngrid,nlay),          &! Saturation index S in Spiga 2012 paper, assuming like in the paper GW phase speed
271                                     !  (stationary waves): c = 0 m.s-1, lambdaH = 150 km, Fo = 7.5e-7 J.m-3
272     SatIndexmap(ngrid),            &! maxval(SatIndex(ig,12:26))
273! ---Misc
274     myT,                           &! Temperature scalar for co2 density computation
275     tcond(ngrid,nlay),             &! CO2 condensation temperature
276     rho_ice_co2T(ngrid,nlay)        ! T-dependant CO2 ice density
277
278  logical :: &
279     file_qext_ok ! Check if file_qext exists
280!======================================================================================================================!
281! BEGIN ===============================================================================================================!
282!======================================================================================================================!
283! 0. Firstcall
284!----------------------------------------------------------------------------------------------------------------------!
285  if (firstcall) then
286    firstcall=.false.
287!----------------------------------------------------------------------------------------------------------------------!
288! 0.1. Initialization of microtimestep from imicroco2
289!----------------------------------------------------------------------------------------------------------------------!
290#ifdef MESOSCALE
291    imicroco2 = 2
292#else
293    imicroco2 = 30
294#endif
295    call getin("imicroco2", imicroco2)
296    microtimestep = ptimestep/real(imicroco2)
297    sigma_iceco2 = sqrt(log(1.+nuiceco2_sed))
298    dev2 = 1. / ( sqrt(2.) * sigma_iceco2 )
299!----------------------------------------------------------------------------------------------------------------------!
300! 0.2. Compute the radius grid of CO2 ice particles (rb_cldco2)
301!        > the grid spacing is computed assuming a constant volume ratio between two consecutive bins; i.e. vrat_cld.
302!   - rad_cldco2 is the primary radius grid used for microphysics computation.
303!   - The grid spacing is computed assuming a constant volume ratio between two consecutive bins; i.e. vrat_cld.
304!   - vrat_cld is determined from the boundary values of the size grid: rmin_cld and rmax_cld.
305!   - The rb_cldco2 array contains the boundary values of each rad_cldco2 bin.
306!----------------------------------------------------------------------------------------------------------------------!
307    ! vrat_cld is determined from the boundary values of the size grid: rmin_cld and rmax_cld.
308    vrat_cld = exp(log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) * 3.)
309
310    ! rad_cldco2 is the primary radius grid used for microphysics computation.
311    rad_cldco2(1) = rmin_cld
312    do i = 1, nbinco2_cld-1
313      rad_cldco2(i+1) = rad_cldco2(i) * vrat_cld**(1./3.)
314    end do
315
316    ! rb_cldco2 array contains the boundary values of each rad_cldco2 bin.
317    rb_cldco2(1) = rbmin_cld
318    do i = 1, nbinco2_cld
319      rb_cldco2(i+1) = ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) * rad_cldco2(i)
320    end do
321    rb_cldco2(nbinco2_cld+1) = rbmax_cld
322!----------------------------------------------------------------------------------------------------------------------!
323! 0.3. Read file 'optprop_co2ice_1mic.dat' to extract optical properties of CO2 ice at 1 micron (Qext)
324!----------------------------------------------------------------------------------------------------------------------!
325    ! get information about file_qext
326    inquire(file=trim(datadir)//'/'//file_qext, exist=file_qext_ok)
327
328    ! if file_qext is missing then stop
329    if (.not. file_qext_ok) then
330      write(*,*)'file'//file_qext//'should be in ', trim(datadir)
331      call abort_physic('co2cloud', 'file missing', 1)
332    end if
333
334    ! read file_qext
335    open(unit=uQext,file=trim(datadir)//'/'//file_qext, form='formatted')
336
337    ! skip 1 line
338    read(uQext,*)
339
340    ! extract radv
341    do i = 1, var_dim_qext
342      read(uQext,'(E12.5)')radv(i)
343    end do
344
345    ! skip 1 line
346    read(uQext,*)
347
348    ! Qextv1mic
349    do i = 1 , var_dim_qext
350      read(uQext,'(E12.5)')Qextv1mic(i)
351    end do
352
353    ! close file_qext
354    close(uQext)
355!----------------------------------------------------------------------------------------------------------------------!
356! 0.4. Interpole the radius grid (rb_cldco2) to get the corresponding exctinction coefficients (Qext1bins), using
357!       file_qext values (radv, Qextv1mic)
358!----------------------------------------------------------------------------------------------------------------------!
359    do i = 1, nbinco2_cld
360      ltemp1 = abs(radv(:)-rb_cldco2(i))
361      ltemp2 = abs(radv(:)-rb_cldco2(i+1))
362      lebon1 = minloc(ltemp1,DIM=1)
363      lebon2 = min(minloc(ltemp2,DIM=1), var_dim_qext)
364      nelem = lebon2 - lebon1 + 1.
365
366      ! mean value in the interval
367      Qtemp = 0d0
368      do l = 0, nelem
369        Qtemp = Qtemp + Qextv1mic(min(lebon1+l, var_dim_qext))
370      end do
371
372      Qext1bins(i) = Qtemp / nelem
373    end do
374
375    Qext1bins(:) = Qext1bins(:) * pi * (rad_cldco2(:)**2)
376
377    ! print result of the interpolation
378    write(*,*)'--------------------------------------------'
379    write(*,*)'Microphysics co2: size bin-Qext information:'
380    write(*,*)'   i, rad_cldco2(i), Qext1bins(i)'
381    do i = 1, nbinco2_cld
382      write(*,'(i3,3x,3(e13.6,4x))')i, rad_cldco2(i), Qext1bins(i)
383    end do
384    write(*,*)'--------------------------------------------'
385!----------------------------------------------------------------------------------------------------------------------!
386! 0.5. Save the radius grid of CO2 particles (rb_cldco2)
387!----------------------------------------------------------------------------------------------------------------------!
388    do i = 1, nbinco2_cld+1
389      rb_cldco2(i) = log(rb_cldco2(i))
390    end do
391  end if ! of IF (firstcall)
392!----------------------------------------------------------------------------------------------------------------------!
393! 1. Initialization
394!----------------------------------------------------------------------------------------------------------------------!
395  sum_subpdq(1:ngrid,1:nlay,1:nq) = 0.
396  sum_subpdt(1:ngrid,1:nlay) = 0.
397
398  subpdqcloudco2(1:ngrid,1:nlay,1:nq) = 0.
399  subpdtcloudco2(1:ngrid,1:nlay) = 0.
400
401  pdqcloudco2(1:ngrid,1:nlay,1:nq) = 0.
402  pdtcloudco2(1:ngrid,1:nlay) = 0.
403
404  ! default value if no ice
405  rhocloudco2(1:ngrid,1:nlay) = rho_dust
406  rhocloudco2t(1:ngrid,1:nlay) = rho_dust
407  epaisseur(1:ngrid,1:nlay) = 0.
408  masse(1:ngrid,1:nlay) = 0.
409
410  zqsed0(1:ngrid,1:nlay,1:nq) = 0.
411  sum_subpdqs_sedco2(1:ngrid) = 0.
412  subpdqsed(1:ngrid,1:nlay,1:nq) = 0.
413!----------------------------------------------------------------------------------------------------------------------!
414! 2. Compute mass and thickness layers
415!----------------------------------------------------------------------------------------------------------------------!
416  do l = 1, nlay
417    do ig = 1, ngrid
418      masse(ig,l) = (pplev(ig,l) - pplev(ig,l+1)) / g
419      epaisseur(ig,l) = pzlev(ig,l+1) - pzlev(ig,l)
420    end do
421  end do
422!----------------------------------------------------------------------------------------------------------------------!
423! 3. Define the sub-grid cloud (CLFvaryingCO2)
424!----------------------------------------------------------------------------------------------------------------------!
425! 3.1. Representation of sub-grid CO2 ice clouds (CLFvaryingCO2 = True)
426!----------------------------------------------------------------------------------------------------------------------!
427  if (CLFvaryingCO2) then
428    ! effective temperature
429    pteff(:,:) = pt(:,:)
430
431    ! min co2cloudfrac when there is ice
432    co2cloudfrac(:,:) = mincloud
433
434    ! temperature
435    do l=1,nlay
436      do ig=1,ngrid
437        zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
438      end do
439    end do
440
441    ! Quantities of traceurs
442    if (igcm_co2 /= 0) then
443      do l = 1, nlay
444        do ig = 1, ngrid
445          zq_co2vap(ig,l) = pq(ig,l,igcm_co2) + pdq(ig,l,igcm_co2)*ptimestep
446        end do
447      end do
448    end if
449!----------------------------------------------------------------------------------------------------------------------!
450! 3.1.a. Saturation index CO2
451!----------------------------------------------------------------------------------------------------------------------!
452    ! if saturation index co2 is true
453    if (satindexco2) then
454      ! layers 12 --> 26 ~ 12->85 km
455      do l = 12, 26
456        do ig = 1, ngrid
457          ! compute N^2 static stability
458          gradT = (zt(ig,l+1)-zt(ig,l))/(pzlev(ig,l+1)-pzlev(ig,l))
459          NN = sqrt(g/zt(iq,l) * (max(gradT,-g/cpp) + g/cpp))
460
461          ! compute absolute value of zonal wind field
462          zu = abs(pu(ig,l) + pdu(ig,l)*ptimestep)
463
464          ! compute background density
465          rho = pplay(ig,l) / (rnew(ig,l)*zt(ig,l))
466
467          ! saturation index: Modulate the DeltaT by GW propagation index:
468          ! --------------------------------------------------------------
469          SatIndex(ig,l) = sqrt(Fo*lambdaH/(2.*pi)*NN / (rho*zu**3) )
470        end do
471      end do
472
473      ! Then compute Satindex map in layers 12 --> 26 ~ 12->85 km
474      do ig = 1, ngrid
475        SatIndexmap(ig) = maxval(SatIndex(ig,12:26))
476      end do
477
478      ! Write outputs in diagfi.nc
479      call WRITEDIAGFI(ngrid, "SatIndex", "SatIndex", " ", 3, SatIndex)
480
481      call WRITEDIAGFI(ngrid, "SatIndexmap", "SatIndexmap", "km", 2, SatIndexmap)
482    !------------------------------------------------------------------------------------------------------------------!
483    ! if saturation index co2 is false, set saturation index to 0.05
484    !------------------------------------------------------------------------------------------------------------------!
485    else
486      do ig = 1, ngrid
487        SatIndexmap(ig)=0.05
488      end do
489    end if ! of if (satindexco2)
490!----------------------------------------------------------------------------------------------------------------------!
491! 3.1.b. Compute tcond
492!----------------------------------------------------------------------------------------------------------------------!
493    call tcondco2(ngrid,nlay,pplay,zq_co2vap,tcond)
494!----------------------------------------------------------------------------------------------------------------------!
495! 3.1.c. Compute cloud fraction in cells
496!----------------------------------------------------------------------------------------------------------------------!
497    do ig = 1, ngrid
498      if (SatIndexmap(ig) <= 0.1) then
499        do l = 1, nlay-1
500
501          ! The entire fraction is saturated
502          if (tcond(ig,l) >= (zt(ig,l)+zdelt) .or. tcond(ig,l) <= 0.) then
503            pteff(ig,l) = zt(ig,l)
504            co2cloudfrac(ig,l) = 1.
505
506          ! No saturation at all
507          else if (tcond(ig,l) <= (zt(ig,l)-zdelt)) then
508            pteff(ig,l) = zt(ig,l) - zdelt
509            co2cloudfrac(ig,l) = mincloud
510
511          ! Mean temperature of the cloud fraction
512          else
513            pteff(ig,l) = (tcond(ig,l)+zt(ig,l)-zdelt) / 2.
514            co2cloudfrac(ig,l) = (tcond(ig,l)-zt(ig,l)+zdelt) / (2.0*zdelt)
515          end if
516
517          pteff(ig,l) = pteff(ig,l) - pdt(ig,l)*ptimestep
518
519          ! check boundary values of co2cloudfrac
520          if (co2cloudfrac(ig,l) <= mincloud) then
521            co2cloudfrac(ig,l) = mincloud
522          else if (co2cloudfrac(ig,l)> 1) then
523            co2cloudfrac(ig,l) = 1.
524          end if
525        end do
526
527      ! SatIndex not favorable for GW: leave pt untouched
528      else
529        pteff(ig,l) = pt(ig,l)
530        co2cloudfrac(ig,l) = mincloud
531      end if ! of if (SatIndexmap <= 0.1)
532    end do ! of ngrid
533    ! TODO: Totalcloud frac of the column missing here
534!----------------------------------------------------------------------------------------------------------------------!
535! 3.2. No sub-grid cloud representation (CLFvarying = False)
536!----------------------------------------------------------------------------------------------------------------------!
537  else
538    do l = 1, nlay
539      do ig = 1, ngrid
540        pteff(ig,l) = pt(ig,l)
541      end do
542    end do
543  end if ! end if (CLFvaryingco2)
544!----------------------------------------------------------------------------------------------------------------------!
545! 4. Microphysics of CO2 cloud formation
546!----------------------------------------------------------------------------------------------------------------------!
547  do microstep = 1, imicroco2
548!----------------------------------------------------------------------------------------------------------------------!
549! 4.1. Stepped entry for tendancies: At each micro timestep we add pdt in order to have a stepped entry
550!----------------------------------------------------------------------------------------------------------------------!
551    do l = 1, nlay
552      do ig = 1, ngrid
553        ! on temperature
554        sum_subpdt(ig,l) = sum_subpdt(ig,l) + pdt(ig,l)
555
556        ! on tracers
557        sum_subpdq(ig,l,igcm_dust_mass) = sum_subpdq(ig,l,igcm_dust_mass) + pdq(ig,l,igcm_dust_mass)
558
559        sum_subpdq(ig,l,igcm_dust_number) = sum_subpdq(ig,l,igcm_dust_number) + pdq(ig,l,igcm_dust_number)
560
561        sum_subpdq(ig,l,igcm_ccnco2_mass) = sum_subpdq(ig,l,igcm_ccnco2_mass) + pdq(ig,l,igcm_ccnco2_mass)
562
563        sum_subpdq(ig,l,igcm_ccnco2_number) = sum_subpdq(ig,l,igcm_ccnco2_number) + pdq(ig,l,igcm_ccnco2_number)
564
565        sum_subpdq(ig,l,igcm_co2_ice) = sum_subpdq(ig,l,igcm_co2_ice) + pdq(ig,l,igcm_co2_ice)
566
567        sum_subpdq(ig,l,igcm_co2) = sum_subpdq(ig,l,igcm_co2) + pdq(ig,l,igcm_co2)
568
569        if (co2useh2o) then
570          sum_subpdq(ig,l,igcm_h2o_ice) = sum_subpdq(ig,l,igcm_h2o_ice) + pdq(ig,l,igcm_h2o_ice)
571
572          sum_subpdq(ig,l,igcm_ccn_mass) = sum_subpdq(ig,l,igcm_ccn_mass) + pdq(ig,l,igcm_ccn_mass)
573
574          sum_subpdq(ig,l,igcm_ccn_number) = sum_subpdq(ig,l,igcm_ccn_number) + pdq(ig,l,igcm_ccn_number)
575        end if
576      end do ! ngrid
577    end do ! nlay
578!----------------------------------------------------------------------------------------------------------------------!
579! 4.2. Effective tracers quantities in the cloud
580!----------------------------------------------------------------------------------------------------------------------!
581    pqeff(:,:,:) = pq(:,:,:)
582
583    if (CLFvaryingCO2) then
584      pqeff(:,:,igcm_ccnco2_mass) = pq(:,:,igcm_ccnco2_mass) / co2cloudfrac(:,:)
585
586      pqeff(:,:,igcm_ccnco2_number) = pq(:,:,igcm_ccnco2_number) / co2cloudfrac(:,:)
587
588      pqeff(:,:,igcm_co2_ice) = pq(:,:,igcm_co2_ice) / co2cloudfrac(:,:)
589    end if
590!----------------------------------------------------------------------------------------------------------------------!
591! 4.3. Gravitational sedimentation
592!----------------------------------------------------------------------------------------------------------------------!
593    if (sedimentation) then
594!----------------------------------------------------------------------------------------------------------------------!
595! 4.3.a. Compute cloud density
596!----------------------------------------------------------------------------------------------------------------------!
597      do l = 1, nlay
598        do ig = 1, ngrid
599          ! temperature during the sedimentation process
600          ztsed(ig,l) = pteff(ig,l) + sum_subpdt(ig,l) * microtimestep
601
602          ! quantities tracers during the sedimentation process
603          zqsed(ig,l,:) = pqeff(ig,l,:) + sum_subpdq(ig,l,:) * microtimestep
604
605          ! density of co2 ice
606          call density_co2_ice(ztsed(ig,l), rho_ice_co2T(ig,l))
607
608          ! assure positive value of co2_ice mmr, ccnco2 number, ccnco2 mass
609          Niceco2 = max(zqsed(ig,l,igcm_co2_ice), threshold)
610          Nccnco2 = max(zqsed(ig,l,igcm_ccnco2_number), threshold)
611          Qccnco2 = max(zqsed(ig,l,igcm_ccnco2_mass), threshold)
612
613          ! Get density cloud and co2 ice particle radius
614          call updaterice_microco2(Niceco2, Qccnco2, Nccnco2, tauscaling(ig), riceco2(ig,l), rhocloudco2t(ig,l))
615
616          ! conflict with updaterice_microco2 C.M.
617          if (Niceco2 <= threshold .or. Nccnco2*tauscaling(ig) <= 1) then
618            riceco2(ig,l) = 1.e-9
619          end if
620
621          ! done in updaterice_microco2 C.M.
622          rhocloudco2t(ig,l) = min( max(rhocloudco2t(ig,l),rho_ice_co2T(ig,l)), rho_dust )
623        end do ! ngrid
624      end do ! nlay
625!----------------------------------------------------------------------------------------------------------------------!
626! 4.3.b. Save actualized tracer values to compute sedimentation tendancies
627!----------------------------------------------------------------------------------------------------------------------!
628      zqsed0(:,:,igcm_co2_ice) = zqsed(:,:,igcm_co2_ice)
629      zqsed0(:,:,igcm_ccnco2_mass) = zqsed(:,:,igcm_ccnco2_mass)
630      zqsed0(:,:,igcm_ccnco2_number) = zqsed(:,:,igcm_ccnco2_number)
631!----------------------------------------------------------------------------------------------------------------------!
632! 4.3.c. Sedimentation of co2 ice
633!----------------------------------------------------------------------------------------------------------------------!
634      rsedcloudco2(1:ngrid, 1:nlay) = 5e-6
635      wq(:,:) = 0.
636      call newsedim(ngrid, nlay, ngrid*nlay, ngrid*nlay, microtimestep, pplev, masse, epaisseur, ztsed, &
637                   rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_co2_ice), wq, beta)
638
639      do ig = 1, ngrid
640        sum_subpdqs_sedco2(ig) = sum_subpdqs_sedco2(ig) + wq(ig,1) / microtimestep
641      end do
642!----------------------------------------------------------------------------------------------------------------------!
643! 4.3.d. Sedimentation for other tracers
644!----------------------------------------------------------------------------------------------------------------------!
645      wq(:,:) = 0.
646      ! for ccnco2_mass
647      call newsedim(ngrid, nlay, ngrid*nlay, ngrid*nlay, microtimestep, pplev, masse, epaisseur, ztsed, &
648                   rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_mass), wq, beta)
649      !TODO: ajouter le calcule de la tendance a la surface comme co2ice
650
651      wq(:,:) = 0.
652      ! for ccnco2_number
653      call newsedim(ngrid, nlay, ngrid*nlay, ngrid*nlay,microtimestep, pplev, masse, epaisseur, ztsed, &
654                   rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_number), wq, beta)
655
656      !TODO: ajouter le calcule de la tendance a la surface comme co2ice
657!----------------------------------------------------------------------------------------------------------------------!
658! 4.3.e. Compute tendencies due to the sedimation process
659!----------------------------------------------------------------------------------------------------------------------!
660      do l = 1, nlay
661        do ig = 1, ngrid
662          subpdqsed(ig,l,igcm_ccnco2_mass) = ( zqsed(ig,l,igcm_ccnco2_mass) - zqsed0(ig,l,igcm_ccnco2_mass)  ) &
663                                             / microtimestep
664
665          subpdqsed(ig,l,igcm_ccnco2_number) = ( zqsed(ig,l,igcm_ccnco2_number) - zqsed0(ig,l,igcm_ccnco2_number) )&
666                                               / microtimestep
667
668          subpdqsed(ig,l,igcm_co2_ice) = ( zqsed(ig,l,igcm_co2_ice) - zqsed0(ig,l,igcm_co2_ice) ) / microtimestep
669        end do
670      end do
671
672      ! update subtimestep tendencies with sedimentation input
673      do l = 1, nlay
674        do ig = 1, ngrid
675          sum_subpdq(ig,l,igcm_ccnco2_mass) = sum_subpdq(ig,l,igcm_ccnco2_mass) + subpdqsed(ig,l,igcm_ccnco2_mass)
676
677          sum_subpdq(ig,l,igcm_ccnco2_number) = sum_subpdq(ig,l,igcm_ccnco2_number) + subpdqsed(ig,l,igcm_ccnco2_number)
678
679          sum_subpdq(ig,l,igcm_co2_ice) = sum_subpdq(ig,l,igcm_co2_ice) + subpdqsed(ig,l,igcm_co2_ice)
680        end do
681      end do
682    end if !(end if sedimentation)
683!----------------------------------------------------------------------------------------------------------------------!
684! 4.4. Main call to the cloud scheme
685!----------------------------------------------------------------------------------------------------------------------!
686    call improvedco2clouds(ngrid, nlay, microtimestep, pplay, pplev, pteff, sum_subpdt, pqeff, sum_subpdq, &
687                           subpdqcloudco2, subpdtcloudco2, nq, tauscaling, mem_Mccn_co2, mem_Mh2o_co2, mem_Nccn_co2, &
688                           rb_cldco2, sigma_iceco2, dev2)
689!----------------------------------------------------------------------------------------------------------------------!
690! 4.5. Updating tendencies after cloud scheme
691!----------------------------------------------------------------------------------------------------------------------!
692    do l = 1, nlay
693      do ig = 1, ngrid
694        sum_subpdt(ig,l) = sum_subpdt(ig,l) + subpdtcloudco2(ig,l)
695
696        sum_subpdq(ig,l,igcm_dust_mass) = sum_subpdq(ig,l,igcm_dust_mass) + subpdqcloudco2(ig,l,igcm_dust_mass)
697
698        sum_subpdq(ig,l,igcm_dust_number) = sum_subpdq(ig,l,igcm_dust_number) + subpdqcloudco2(ig,l,igcm_dust_number)
699
700        sum_subpdq(ig,l,igcm_ccnco2_mass) = sum_subpdq(ig,l,igcm_ccnco2_mass) + subpdqcloudco2(ig,l,igcm_ccnco2_mass)
701
702        sum_subpdq(ig,l,igcm_ccnco2_number) = sum_subpdq(ig,l,igcm_ccnco2_number) + &
703                                              subpdqcloudco2(ig,l,igcm_ccnco2_number)
704
705        sum_subpdq(ig,l,igcm_co2_ice) = sum_subpdq(ig,l,igcm_co2_ice) + subpdqcloudco2(ig,l,igcm_co2_ice)
706
707        sum_subpdq(ig,l,igcm_co2) = sum_subpdq(ig,l,igcm_co2) + subpdqcloudco2(ig,l,igcm_co2)
708
709        if (co2useh2o) then
710          sum_subpdq(ig,l,igcm_h2o_ice) = sum_subpdq(ig,l,igcm_h2o_ice) + subpdqcloudco2(ig,l,igcm_h2o_ice)
711
712          sum_subpdq(ig,l,igcm_ccn_mass) = sum_subpdq(ig,l,igcm_ccn_mass) + subpdqcloudco2(ig,l,igcm_ccn_mass)
713
714          sum_subpdq(ig,l,igcm_ccn_number) = sum_subpdq(ig,l,igcm_ccn_number) + subpdqcloudco2(ig,l,igcm_ccn_number)
715        end if
716      end do ! ngrid
717    end do ! nlay
718  end do ! of do microstep = 1, imicroco2
719!----------------------------------------------------------------------------------------------------------------------!
720! 5. Compute final tendencies after time loop
721!----------------------------------------------------------------------------------------------------------------------!
722! condensation/sublimation rate in the atmosphere
723  do l = nlay, 1, -1
724    do ig = 1, ngrid
725      pcondicea(ig,l) = sum_subpdq(ig,l,igcm_co2_ice) / real(imicroco2)
726    end do
727  end do
728
729  ! CO2 flux at surface (kg.m-2.s-1)
730  do ig = 1, ngrid
731    pdqs_sedco2(ig) = sum_subpdqs_sedco2(ig) / real(imicroco2)
732  end do
733
734  ! temperature tendency
735  do l = 1, nlay
736    do ig = 1, ngrid
737      pdtcloudco2(ig,l) = ( sum_subpdt(ig,l)/real(imicroco2) ) - pdt(ig,l)
738    end do
739  end do
740
741  ! tracers tendencies
742  do l = 1, nlay
743    do ig = 1, ngrid
744      pdqcloudco2(ig,l,igcm_co2) = 0. ! here is the trick, this tendency is computed in co2condens_mod4micro
745
746      pdqcloudco2(ig,l,igcm_co2_ice) = ( sum_subpdq(ig,l,igcm_co2_ice) / real(imicroco2) ) - pdq(ig,l,igcm_co2_ice)
747
748      pdqcloudco2(ig,l,igcm_ccnco2_mass) = ( sum_subpdq(ig,l,igcm_ccnco2_mass)/real(imicroco2) ) - &
749                                           pdq(ig,l,igcm_ccnco2_mass)
750
751      pdqcloudco2(ig,l,igcm_ccnco2_number) = ( sum_subpdq(ig,l,igcm_ccnco2_number) / real(imicroco2) ) - &
752                                             pdq(ig,l,igcm_ccnco2_number)
753
754      pdqcloudco2(ig,l,igcm_dust_mass) = ( sum_subpdq(ig,l,igcm_dust_mass) / real(imicroco2) ) - &
755                                         pdq(ig,l,igcm_dust_mass)
756
757      pdqcloudco2(ig,l,igcm_dust_number) = ( sum_subpdq(ig,l,igcm_dust_number)/real(imicroco2) ) - &
758                                           pdq(ig,l,igcm_dust_number)
759
760      if (co2useh2o) then
761        pdqcloudco2(ig,l,igcm_h2o_ice) = ( sum_subpdq(ig,l,igcm_h2o_ice) / real(imicroco2) ) -&
762                                         pdq(ig,l,igcm_h2o_ice)
763
764        pdqcloudco2(ig,l,igcm_ccn_mass) = ( sum_subpdq(ig,l,igcm_ccn_mass) / real(imicroco2) ) - &
765                                          pdq(ig,l,igcm_ccn_mass)
766
767        pdqcloudco2(ig,l,igcm_ccn_number) = ( sum_subpdq(ig,l,igcm_ccn_number) / real(imicroco2) ) - &
768                                            pdq(ig,l,igcm_ccn_number)
769      end if
770    end do ! ngrid
771  end do ! nlay
772!----------------------------------------------------------------------------------------------------------------------!
773! 6. Update clouds physical values in the cloud (for output)
774!----------------------------------------------------------------------------------------------------------------------!
775! 6.1. Update density of co2 ice, riceco2 and opacity
776!----------------------------------------------------------------------------------------------------------------------!
777  do l = 1, nlay
778    do ig = 1, ngrid
779      Niceco2 = pqeff(ig,l,igcm_co2_ice) + ( pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice) ) * ptimestep
780
781      Nccnco2 = max( (pqeff(ig,l,igcm_ccnco2_number) + (pdq(ig,l,igcm_ccnco2_number) + &
782                       pdqcloudco2(ig,l, igcm_ccnco2_number)) * ptimestep) &
783                    , threshold)
784
785      Qccnco2 = max( (pqeff(ig,l,igcm_ccnco2_mass) + (pdq(ig,l,igcm_ccnco2_mass) + &
786                       pdqcloudco2(ig,l, igcm_ccnco2_mass)) * ptimestep)&
787                    , threshold)
788
789      myT = pteff(ig,l) + (pdt(ig,l)+pdtcloudco2(ig,l))*ptimestep
790
791      ! compute density of co2 ice
792      call density_co2_ice(myT, rho_ice_co2T(ig,l))
793
794      rho_ice_co2 = rho_ice_co2T(ig,l) ! rho_ice_co2 is shared by tracer_mod and used in updaterice
795
796      ! Compute particle size
797      call updaterice_microco2(Niceco2, Qccnco2, Nccnco2, tauscaling(ig), riceco2(ig,l), rhocloudco2(ig,l))
798
799     ! Compute opacities
800      if ( (Niceco2 <= threshold .or. Nccnco2*tauscaling(ig) <= 1.) ) then
801        riceco2(ig,l) = 0.
802        Qext1bins2(ig,l) = 0.
803      else
804        n_derf = derf( (rb_cldco2(1)-dlog(riceco2(ig,l))) *dev2)
805        Qext1bins2(ig,l) = 0.
806
807        do i = 1, nbinco2_cld
808          n_aer(i) = -0.5 * Nccnco2*tauscaling(ig) * n_derf
809
810          n_derf = derf((rb_cldco2(i+1)-dlog(riceco2(ig,l))) *dev2)
811          n_aer(i) = n_aer(i) + (0.5 * Nccnco2*tauscaling(ig) * n_derf)
812
813          Qext1bins2(ig,l) = Qext1bins2(ig,l) + Qext1bins(i) * n_aer(i)
814        end do
815      end if
816!----------------------------------------------------------------------------------------------------------------------!
817! 6.2. Update rice and rdust
818!----------------------------------------------------------------------------------------------------------------------!
819      ! update rice water only if co2 use h2o ice as CCN
820      if (co2useh2o) then
821        call updaterice_micro( &
822              pqeff(ig,l,igcm_h2o_ice) + (pdq(ig,l,igcm_h2o_ice) + pdqcloudco2(ig,l,igcm_h2o_ice))*ptimestep, &
823              pqeff(ig,l,igcm_ccn_mass) + (pdq(ig,l,igcm_ccn_mass) + pdqcloudco2(ig,l,igcm_ccn_mass))*ptimestep, &
824              pqeff(ig,l,igcm_ccn_number) + (pdq(ig,l,igcm_ccn_number) + pdqcloudco2(ig,l,igcm_ccn_number))*ptimestep, &
825              tauscaling(ig),rice(ig,l),rhocloud(ig,l))
826      end if
827
828      ! update rdust
829      call updaterdust( &
830           pqeff(ig,l,igcm_dust_mass) + (pdq(ig,l,igcm_dust_mass) + pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep, &
831           pqeff(ig,l,igcm_dust_number) + (pdq(ig,l,igcm_dust_number) + pdqcloudco2(ig,l,igcm_dust_number))*ptimestep, &
832           rdust(ig,l))
833    end do ! ngrid
834  end do ! nlay
835!----------------------------------------------------------------------------------------------------------------------!
836! 7. A correction if a lot of subliming CO2 fills the 1st layer FF (04/2005). Then that should not affect the ice
837!      particle radius
838!----------------------------------------------------------------------------------------------------------------------!
839  do ig = 1, ngrid
840    if ( pdpsrf(ig)*ptimestep > 0.9*(pplev(ig,1)-pplev(ig,2))) then
841
842      if ( pdpsrf(ig)*ptimestep > 0.9*(pplev(ig,1)-pplev(ig,3)) ) then
843        riceco2(ig,2) = riceco2(ig,3)
844      end if
845
846      riceco2(ig,1) = riceco2(ig,2)
847    end if
848  end do
849!----------------------------------------------------------------------------------------------------------------------!
850! 8. Compute water cloud sedimentation radius
851!----------------------------------------------------------------------------------------------------------------------!
852  do l = 1, nlay
853    do ig= 1, ngrid
854      rsedcloud(ig,l) = max( rice(ig,l)*(1.+nuice_sed)* (1.+nuice_sed)*(1.+nuice_sed), rdust(ig,l) )
855    end do
856  end do
857!----------------------------------------------------------------------------------------------------------------------!
858! 9. CO2 saturated quantities
859!----------------------------------------------------------------------------------------------------------------------!
860! 9.1 Compute CO2 saturation in layers
861!----------------------------------------------------------------------------------------------------------------------!
862  call co2sat(ngrid*nlay, pteff+(pdt+pdtcloudco2)*ptimestep, zqsatco2)
863!----------------------------------------------------------------------------------------------------------------------!
864! 9.2 Compute CO2 saturated quantities in layers
865!----------------------------------------------------------------------------------------------------------------------!
866  do l = 1, nlay
867    do ig = 1, ngrid
868      satuco2(ig,l) = ( pqeff(ig,l,igcm_co2)  + (pdq(ig,l,igcm_co2) + pdqcloudco2(ig,l,igcm_co2))*ptimestep ) *  &
869                      (mmean(ig,l)/(mco2*1e3)) * pplay(ig,l) / zqsatco2(ig,l)
870    end do
871  end do
872!----------------------------------------------------------------------------------------------------------------------!
873! 10. Everything modified by CO2 microphysics must be wrt co2cloudfrac
874!----------------------------------------------------------------------------------------------------------------------!
875  if (CLFvaryingCO2) then
876    do l = 1, nlay
877      do ig = 1, ngrid
878        pdqcloudco2(ig,l,igcm_ccnco2_mass) = pdqcloudco2(ig,l,igcm_ccnco2_mass) * co2cloudfrac(ig,l)
879
880        pdqcloudco2(ig,l,igcm_ccnco2_number) = pdqcloudco2(ig,l,igcm_ccnco2_number) * co2cloudfrac(ig,l)
881
882        pdqcloudco2(ig,l,igcm_dust_mass) = pdqcloudco2(ig,l,igcm_dust_mass) * co2cloudfrac(ig,l)
883
884        pdqcloudco2(ig,l,igcm_dust_number) = pdqcloudco2(ig,l,igcm_dust_number) * co2cloudfrac(ig,l)
885
886        pdqcloudco2(ig,l,igcm_co2_ice) = pdqcloudco2(ig,l,igcm_co2_ice) * co2cloudfrac(ig,l)
887
888        pdqcloudco2(ig,l,igcm_co2) = pdqcloudco2(ig,l,igcm_co2) * co2cloudfrac(ig,l)
889
890        pdtcloudco2(ig,l) = pdtcloudco2(ig,l) * co2cloudfrac(ig,l)
891
892        Qext1bins2(ig,l) = Qext1bins2(ig,l) * co2cloudfrac(ig,l)
893
894        if (co2useh2o) then
895          pdqcloudco2(ig,l,igcm_h2o_ice) = pdqcloudco2(ig,l,igcm_h2o_ice) * co2cloudfrac(ig,l)
896
897          pdqcloudco2(ig,l,igcm_ccn_mass) = pdqcloudco2(ig,l,igcm_ccn_mass) * co2cloudfrac(ig,l)
898
899          pdqcloudco2(ig,l,igcm_ccn_number) = pdqcloudco2(ig,l,igcm_ccn_number) * co2cloudfrac(ig,l)
900        end if
901      end do ! ngrid
902    end do ! nlay
903  end if ! if CLFvaryingCO2 is true
904!----------------------------------------------------------------------------------------------------------------------!
905! 11. Compute opacity at 1 micron: Opacity in mesh ig is the sum over l of Qext1bins2. Is this true ?
906!----------------------------------------------------------------------------------------------------------------------!
907  tau1mic(:)=0.
908  do l = 1, nlay
909    do ig = 1, ngrid
910      tau1mic(ig) = tau1mic(ig) + Qext1bins2(ig,l)
911    end do
912  end do
913!----------------------------------------------------------------------------------------------------------------------!
914! 12. Write outputs in diagfi.nc
915!----------------------------------------------------------------------------------------------------------------------!
916  call WRITEDIAGFI(ngrid, "satuco2", "vap in satu", " ", 3, satuco2)
917
918  call WRITEDIAGFI(ngrid, "co2cloudfrac", "co2 cloud fraction", " ", 3, co2cloudfrac)
919
920  call writediagfi(ngrid, "riceco2", "ice radius", "m", 3, riceco2)
921
922  call WRITEDIAGFI(ngrid, "Tau3D1mic", " co2 ice opacities", " ", 3, Qext1bins2)
923
924  call WRITEDIAGFI(ngrid, "tau1mic", "co2 ice opacity 1 micron", " ", 2, tau1mic)
925
926  call WRITEDIAGFI(ngrid, "mem_Nccn_co2", "CCN number used by CO2", "kg/kg", 3, mem_Nccn_co2)
927
928  call WRITEDIAGFI(ngrid, "mem_Mccn_co2", "CCN mass used by CO2",  "kg/kg", 3, mem_Mccn_co2)
929
930  call WRITEDIAGFI(ngrid, "mem_Mh2o_co2", "H2O mass in CO2 crystal", "kg/kg", 3, mem_Mh2o_co2)
931!======================================================================================================================!
932! END =================================================================================================================!
933!======================================================================================================================!
934  end subroutine co2cloud
935
936
937!**********************************************************************************************************************!
938!**********************************************************************************************************************!
939
940
941!======================================================================================================================!
942! SUBROUTINE: ini_co2cloud ============================================================================================!
943!======================================================================================================================!
944! Subject:
945!---------
946!   Allocate arrays used for co2 microphysics
947!======================================================================================================================!
948  subroutine ini_co2cloud(ngrid,nlayer)
949
950  implicit none
951
952  integer, intent(in) :: &
953     ngrid, &! number of atmospheric columns
954     nlayer  ! number of atmospheric layers
955
956  allocate(mem_Nccn_co2(ngrid,nlayer))
957  allocate(mem_Mccn_co2(ngrid,nlayer))
958  allocate(mem_Mh2o_co2(ngrid,nlayer))
959
960  end subroutine ini_co2cloud
961
962!**********************************************************************************************************************!
963!**********************************************************************************************************************!
964
965!======================================================================================================================!
966! SUBROUTINE: end_co2cloud ============================================================================================!
967!======================================================================================================================!
968! Subject:
969!---------
970!   Deallocate arrays used for co2 microphysics
971!======================================================================================================================!
972  subroutine end_co2cloud
973
974  implicit none
975
976  if (allocated(mem_Nccn_co2)) deallocate(mem_Nccn_co2)
977  if (allocated(mem_Mccn_co2)) deallocate(mem_Mccn_co2)
978  if (allocated(mem_Mh2o_co2)) deallocate(mem_Mh2o_co2)
979
980  end subroutine end_co2cloud
981
982end module co2cloud_mod
Note: See TracBrowser for help on using the repository browser.