source: trunk/LMDZ.MARS/libf/phymars/growthrate.F @ 524

Last change on this file since 524 was 358, checked in by aslmd, 13 years ago

LMDZ.MARS:

A FIRST VERSION WITH SCAVENGING
TRANSPARENT TO CASUAL USER
FOR WATER CYCLE FOLKS, SEE BELOW

[TN and AS solved conflicts due to r330 >> r357]

07/11/11 == JBM

Changed watercloud.F to call two separate routines,

simpleclouds.F or improvedclouds.F, which are a simplified and
full microphysical scheme for cloud formation, respectively.
Removed the tag called "improved" in watercloud.F, and added
another tag called "microphys" which is defined in callphys.F
instead. Changed routines: callkeys, inifis, physiq, watercloud.

Reimplemented the use of typical profiles of dust particle sizes

and CCNs in simpleclouds.F. Corrected the previously used
analytical CCN profile. Moved ccn_factor to simpleclouds.F,
which won't be used in the new microphysical scheme. Changed
routines: aeropacity, physiq, simpleclouds, watercloud.

Computed rdust at the end of callsedim instead of updatereffrad,

except at firstcall. Renamed rfall into rsedcloud and moved it
in simpleclouds. Moved nuice_sed to initracer.F and added it to
tracer.h. Changed routines: callsedim, physiq, tracer.h,
watercloud, initracer, simpleclouds, updatereffrad.

Added two tracers for the CCNs. Added the different tests in

initracer.F to make sure that, for example, microphys is not
used without doubleq, etc. Corrected an inconsistency in
callsedim.F, and changed the way r0 is updated. Changes
routines: callsedim, inifis, initracer, physiq, testphys1d,
tracer.h.

Added the ability to have a spatially variable density in

newsedim (same method as that used for the radius of
sedimentation). Required the addition of one input to newsedim,
which is the size of the density variable "rho". Changed
routines: callsedim, newsedim.

Added an output to aeropacity called "tauscaling", which is a

factor to convert dust_mass and dust_number to true quantities,
based on tauvis. Removed ccn and qdust from aeropacity, which
are obsolete now.

Wrote improvedcloud.F which includes all the microphysical

scheme, and connected it to the sedimentation scheme. Added and
changed routines: callsedim, physiq, growthrate, nuclea,
improvedclouds, initracer, watercloud, microphys.h.

07/11/11 == TN

Added CCN & dust tracers updates in physiq.F

Corrected a bug that can give negative CCN numbers, due to the
use of single precision values (only 6 decimals) whereas up to 10E+10
CCN can be formed or disappears...

Corrected physical bug that causes h2o_vap or h2o_ice

to be negative in improvedclouds.F.

Corrected physical bug that causes CCN & dust tracers

to be negative in improvedclouds.F when all ice sublimates and
releases dust

Added parameter contact mteta in callphys.def

Default value is still 0.95, see inifis.F

Changed tendancy computation for dust_number in improvedclouds.F

that was not the right one. Indeed, the scavenged dust_number tracer
is computed from the dust_mass one, and its tendancy before scavenging
must be taken into account to compute its scavenging's tendancy.
File size: 3.1 KB
Line 
1      subroutine growthrate(timestep,temp,pmid,ph2o,psat,
2     &                      seq,rcrystal,growth)
3
4      IMPLICIT NONE
5
6c=======================================================================
7c
8c     Determination of the water ice crystal growth rate
9c
10c     Authors: F. Montmessin
11c       Adapted for the LMD/GCM by J.-B. Madeleine (October 2011)
12c     
13c=======================================================================
14
15c-----------------------------------------------------------------------
16c   declarations:
17c   -------------
18
19#include "dimensions.h"
20#include "dimphys.h"
21#include "comcstfi.h"
22#include "tracer.h"
23#include "microphys.h"
24
25c
26c   arguments:
27c   ----------
28
29      REAL timestep
30      REAL temp     ! temperature in the middle of the layer (K)
31      REAL pmid     ! pressure in the middle of the layer (K)
32      REAL*8 ph2o   ! water vapor partial pressure (Pa)
33      REAL*8 psat   ! water vapor saturation pressure (Pa)
34      REAL rcrystal ! crystal radius before condensation (m)
35      REAL*8 seq    ! Equilibrium saturation ratio
36      REAL*8 growth
37
38c   local:
39c   ------
40
41      REAL*8 k,Lv                 
42      REAL*8 knudsen           ! Knudsen number (gas mean free path/particle radius)
43      REAL*8 afactor,Dv,lambda       ! Intermediate computations for growth rate
44      REAL*8 Rk,Rd
45
46c-----------------------------------------------------------------------
47c      Ice particle growth rate by diffusion/impegement of water molecules
48c                r.dr/dt = (S-Seq) / (Seq*Rk+Rd)
49c        with r the crystal radius, Rk and Rd the resistances due to
50c        latent heat release and to vapor diffusion respectively
51c-----------------------------------------------------------------------
52
53c     - Equilibrium saturation accounting for KeLvin Effect
54c      seq=exp(2*sigh2o*mh2o/(rho_ice*rgp*t*r))
55c      (already computed in improvedcloud.F)
56
57c     - Thermal conductibility of CO2
58      k  = (0.17913 * temp - 13.9789) * 4.184e-4
59c     - Latent heat of h2o (J.kg-1)
60      Lv = (2834.3 - 0.28 * (temp-To) - 0.004 * (temp-To)**2 ) * 1.e+3
61
62c     - Constant to compute gas mean free path
63c     l= (T/P)*a, with a = (  0.707*8.31/(4*pi*molrad**2 * avogadro))
64      afactor = 0.707*rgp/(4 * pi* molco2**2  * nav)
65
66c     - Compute Dv, water vapor diffusion coefficient
67c       accounting for both kinetic and continuum regime of diffusion,
68c       the nature of which depending on the Knudsen number.
69
70      Dv = 1./3. * sqrt( 8*kbz*temp/(pi*mh2o/nav) )* kbz * temp /
71     &   ( pi * pmid * (molco2+molh2o)**2 * sqrt(1.+mh2o/mco2) )
72
73      knudsen = temp / pmid * afactor / rcrystal
74      lambda  = (1.333+0.71/knudsen) / (1.+1./knudsen)
75      Dv      = Dv / (1. + lambda * knudsen)
76
77c     - Compute Rk
78      Rk = Lv**2 * rho_ice * mh2o / (k*rgp*temp**2.)
79c     - Compute Rd
80      Rd = rgp * temp *rho_ice / (Dv*psat*mh2o)
81
82c     - Compute growth=rdr/dt, then r(t+1)= sqrt(r(t)**2.+2.*growth*dt)
83      growth = 1. / (seq*Rk+Rd)
84c       growth = (ph2o/psat-seq) / (seq*Rk+Rd)
85c      rf   = sqrt( max( r**2.+2.*growth*timestep , 0. ) )
86c      dr   = rf-r
87
88      RETURN
89      END
90
Note: See TracBrowser for help on using the repository browser.