source: trunk/LMDZ.MARS/libf/phymars/nuclea.F @ 522

Last change on this file since 522 was 520, checked in by tnavarro, 14 years ago

10/02/12 == TN

Major update on watercycle: a smaller integration timestep is now used

in watercloud.F, sedimentation of clouds is done in watercloud instead of
callsedim.F

Temperature-dependant contact parameter in nuclea.F
No dust lifting if CO2 ice in vdif.c
Ice integrated column opacity is written in diagfi from physiq.F, instead

of aeropacity.F. Mandatory if iradia is not 1.

New definition of permanent ice in surfini.F and possibility to have an ice

cap in it in 1d.

Update in deftank: callphys.def.outliers,run.def.1d; added traceur.def.scavenging
File size: 4.9 KB
RevLine 
[358]1*******************************************************
2*                                                     *
3      subroutine nuclea(ph2o,temp,sat,n_ccn,nucrate)
4      implicit none
5*                                                     *
6*   This subroutine computes the nucleation rate      *
7*   as given in Pruppacher & Klett (1978) in the      *
8*   case of water ice forming on a solid substrate.   *
9*     Definition refined by Keese (jgr,1989)          *
10*   Authors: F. Montmessin                            *
11*     Adapted for the LMD/GCM by J.-B. Madeleine      *
12*     (October 2011)                                  *
13*******************************************************
14
15#include "dimensions.h"
16#include "dimphys.h"
17#include "comcstfi.h"
18#include "tracer.h"
19#include "microphys.h"
20
21c     Inputs
22      DOUBLE PRECISION ph2o,sat
23      DOUBLE PRECISION n_ccn(nbin_cld)
24      REAL temp
25
26c     Output
27      DOUBLE PRECISION nucrate(nbin_cld)
28
29c     Local variables
30      DOUBLE PRECISION nh2o
31      DOUBLE PRECISION sig      ! Water-ice/air surface tension  (N.m)
32      external sig
33      DOUBLE PRECISION rstar    ! Radius of the critical germ (m)
34      DOUBLE PRECISION gstar    ! # of molecules forming a critical embryo
35      DOUBLE PRECISION fistar   ! Activation energy required to form a critical embryo (J)
36      DOUBLE PRECISION zeldov   ! Zeldovitch factor (no dim)
37      DOUBLE PRECISION fshape   ! function defined at the end of the file
38      DOUBLE PRECISION deltaf
39
40c     Ratio rstar/radius of the nucleating dust particle
41      double precision xratio
42     
43      double precision mtetalocal ! local mteta in double precision
44
45      integer i
[520]46     
47      LOGICAL firstcall
48      DATA firstcall/.true./
49      SAVE firstcall
[358]50
51c     *************************************************
52
53      mtetalocal = mteta
54
[520]55cccccccccccccccccccccccccccccccccccccccccccccccccc
56ccccccccccc ESSAIS TN MTETA = F (T) cccccccccccccc
57c      if (temp .gt. 200) then
58c         mtetalocal = mtetalocal
59c      else if (temp .lt. 190) then
60c         mtetalocal = mtetalocal-0.05
61c      else
62c         mtetalocal = mtetalocal - (190-temp)*0.005
63c      endif
64c----------------exp law, see Trainer 2008, J. Phys. Chem. C 2009, 113, 2036\u20132040
65       !mtetalocal = max(mtetalocal - 6005*exp(-0.065*temp),0.1)
66       !mtetalocal = max(mtetalocal - 6005*exp(-0.068*temp),0.1)
67               !print*, mtetalocal, temp
68cccccccccccccccccccccccccccccccccccccccccccccccccc
69cccccccccccccccccccccccccccccccccccccccccccccccccc
70      IF (firstcall) THEN
71          print*, ' ' 
72          print*, 'dear user, please keep in mind that'
73          print*, 'contact parameter IS constant'
74          !print*, 'contact parameter IS NOT constant:'
75          !print*, 'max(mteta - 6005*exp(-0.065*temp),0.1)'
76          !print*, 'max(mteta - 6005*exp(-0.068*temp),0.1)'
77          print*, ' ' 
78         firstcall=.false.
79      END IF
80cccccccccccccccccccccccccccccccccccccccccccccccccc
81cccccccccccccccccccccccccccccccccccccccccccccccccc
82   
83
[358]84      if (sat .gt. 1.) then    ! minimum condition to activate nucleation
85
86        nh2o   = ph2o / kbz / temp
87        rstar  = 2. * sig(temp) * vo1 / (rgp*temp*dlog(sat))
88        gstar  = 4. * nav * pi * (rstar**3) / (3.*vo1)
89
90c       Loop over size bins
91        do 200 i=1,nbin_cld
92
93          if ( n_ccn(i) .lt. 1e-10 ) then
94c           no dust, no need to compute nucleation!
95            nucrate(i)=0.
96            goto 200
97          endif
98
99          xratio   = rad_cld(i) / rstar
100          fistar = (4./3.*pi) * sig(temp) * (rstar**2.) *
101     &             fshape(mtetalocal,xratio)
102          deltaf = (2.*desorp-surfdif-fistar)/
103     &             (kbz*temp)
104          deltaf = min( max(deltaf, -100.), 100.)
105
106          if (deltaf.eq.-100.) then
107            nucrate(i) = 0.
108          else
109            zeldov = sqrt ( fistar /
110     &               (3.*pi*kbz*temp*(gstar**2.)) )
111            nucrate(i)= zeldov * kbz * temp * rstar
112     &                  * rstar * 4. * pi
113     &                  * ( nh2o*rad_cld(i) )**2.
114     &                  / ( fshape(mtetalocal,xratio) * nus * m0 )
115     &                  * dexp (deltaf)
116          endif
117
118200     continue
119
120      else
121
122        do i=1,nbin_cld
123          nucrate(i) = 0.
124        enddo
125
126      endif
127
128      return
129      end
130
131*********************************************************
132      double precision function fshape(cost,rap)
133      implicit none
134*        function computing the f(m,x) factor           *
135* related to energy required to form a critical embryo  *
136*********************************************************
137
138      double precision cost,rap
139      double precision phi
140      double precision a,b,c
141
142      phi = sqrt( 1. - 2.*cost*rap + rap**2 )
143      a = 1. + ( (1.-cost*rap)/phi )**3
144      b = (rap**3) * (2.-3.*(rap-cost)/phi+((rap-cost)/phi)**3)
145      c = 3. * cost * (rap**2) * ((rap-cost)/phi-1.)
146
147      fshape = 0.5*(a+b+c)
148
149      if (rap.gt.3000.) fshape = ((2.+cost)*(1.-cost)**2)/4.
150
151      return
152      end
Note: See TracBrowser for help on using the repository browser.