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

Last change on this file since 2561 was 2561, checked in by jnaar, 3 years ago

Water cycle flag update & addition :

  • "cap_albedo" renamed "cst_cap_albedo" (default false) : if true, water ice cap albedo remains constant even when frost with higher albedo condensates on it
  • "refill_watercap" added (default false) : turns h2o_ice_s into watercap when above a given threshold
  • "frost_metam_threshold" added (default 0.05 km.m-2) : threshold used by "refill_watercap"

JN
File size: 6.0 KB
Line 
1*******************************************************
2*                                                     *
3      subroutine nuclea(ph2o,temp,sat,n_ccn,nucrate)
4      USE comcstfi_h
5      implicit none
6*                                                     *
7*   This subroutine computes the nucleation rate      *
8*   as given in Pruppacher & Klett (1978) in the      *
9*   case of water ice forming on a solid substrate.   *
10*     Definition refined by Keese (jgr,1989)          *
11*   Authors: F. Montmessin                            *
12*     Adapted for the LMD/GCM by J.-B. Madeleine      *
13*     (October 2011)                                  *
14*     Optimisation by A. Spiga (February 2012)        * 
15*******************************************************
16
17#include "microphys.h"
18      include "callkeys.h"
19
20c     Inputs
21      DOUBLE PRECISION ph2o,sat
22      DOUBLE PRECISION n_ccn(nbin_cld)
23      REAL temp
24
25c     Output
26   !   DOUBLE PRECISION nucrate(nbin_cld)
27      REAL 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
41c     double precision xratio
42     
43      double precision mtetalocal ! local mteta in double precision
44
45      double precision fshapesimple,zefshape
46
47
48      integer i
49     
50      LOGICAL firstcall
51      DATA firstcall/.true./
52      SAVE firstcall
53
54c     *************************************************
55
56      mtetalocal = mteta  !! use mtetalocal for better performance
57
58      IF (temp_dependant_m) THEN
59c     Simple linear parametrisation from Maattaanen 2014
60c     Smectite sample
61c     Maxed out at 0.97 for physical realism
62         mtetalocal = min(0.0044*temp + 0.1831,0.97)
63      ENDIF ! (temp_dependant_m) THEN
64cccccccccccccccccccccccccccccccccccccccccccccccccc
65ccccccccccc ESSAIS TN MTETA = F (T) cccccccccccccc
66c      if (temp .gt. 200) then
67c         mtetalocal = mtetalocal
68c      else if (temp .lt. 190) then
69c         mtetalocal = mtetalocal-0.05
70c      else
71c         mtetalocal = mtetalocal - (190-temp)*0.005
72c      endif
73c----------------exp law, see Trainer 2008, J. Phys. Chem. C 2009, 113, 2036\u20132040
74       !mtetalocal = max(mtetalocal - 6005*exp(-0.065*temp),0.1)
75       !mtetalocal = max(mtetalocal - 6005*exp(-0.068*temp),0.1)
76               !print*, mtetalocal, temp
77cccccccccccccccccccccccccccccccccccccccccccccccccc
78cccccccccccccccccccccccccccccccccccccccccccccccccc
79      IF (firstcall.and.temp_dependant_m) THEN
80          print*, ' ' 
81          print*, 'dear user, please keep in mind that'
82          print*, 'contact parameter IS NOT constant ;'
83          print*, 'Using the following linear fit from'
84          print*, 'Maattanen et al. 2014 (SM linear fit) :'
85          print*, 'min(0.0044*temp + 0.1831,0.97)'
86          print*, ' ' 
87         firstcall=.false.
88      ELSE IF (firstcall.and.(.not.(temp_dependant_m))) THEN
89          print*, ' ' 
90          print*, 'dear user, please keep in mind that'
91          print*, 'contact parameter IS constant'
92          print*, ' ' 
93         firstcall=.false.
94      END IF
95cccccccccccccccccccccccccccccccccccccccccccccccccc
96cccccccccccccccccccccccccccccccccccccccccccccccccc
97   
98
99      if (sat .gt. 1.) then    ! minimum condition to activate nucleation
100
101        nh2o   = ph2o / kbz / temp
102        rstar  = 2. * sig(temp) * vo1 / (rgp*temp*log(sat))
103        gstar  = 4. * nav * pi * (rstar * rstar * rstar) / (3.*vo1)
104       
105        fshapesimple = (2.+mtetalocal)*(1.-mtetalocal)*(1.-mtetalocal)
106     &                   / 4.
107
108c       Loop over size bins
109        do 200 i=1,nbin_cld
110
111          if ( n_ccn(i) .lt. 1e-10 ) then
112c           no dust, no need to compute nucleation!
113            nucrate(i)=0.
114            goto 200
115          endif
116
117          if (rad_cld(i).gt.3000.*rstar) then
118            zefshape = fshapesimple
119          else
120            zefshape = fshape(mtetalocal,rad_cld(i)/rstar)
121          endif
122
123          fistar = (4./3.*pi) * sig(temp) * (rstar * rstar) *
124     &             zefshape
125          deltaf = (2.*desorp-surfdif-fistar)/
126     &             (kbz*temp)
127          deltaf = min( max(deltaf, -100.d0), 100.d0)
128
129          if (deltaf.eq.-100.) then
130            nucrate(i) = 0.
131          else
132            nucrate(i)= real(sqrt ( fistar /
133     &               (3.*pi*kbz*temp*(gstar*gstar)) )
134     &                  * kbz * temp * rstar
135     &                  * rstar * 4. * pi
136     &                  * ( nh2o*rad_cld(i) )
137     &                  * ( nh2o*rad_cld(i) )
138     &                  / ( zefshape * nus * m0 )
139     &                  * exp (deltaf))
140          endif
141
142200     continue
143
144      else
145
146        do i=1,nbin_cld
147          nucrate(i) = 0.
148        enddo
149
150      endif
151
152      return
153      end
154
155*********************************************************
156      double precision function fshape(cost,rap)
157      implicit none
158*        function computing the f(m,x) factor           *
159* related to energy required to form a critical embryo  *
160*********************************************************
161
162      double precision cost,rap
163      double precision yeah
164
165          !! PHI
166          yeah = sqrt( 1. - 2.*cost*rap + rap*rap )
167          !! FSHAPE = TERM A
168          fshape = (1.-cost*rap) / yeah
169          fshape = fshape * fshape * fshape
170          fshape = 1. + fshape
171          !! ... + TERM B
172          yeah = (rap-cost)/yeah
173          fshape = fshape +
174     & rap*rap*rap*(2.-3.*yeah+yeah*yeah*yeah)
175          !! ... + TERM C
176          fshape = fshape + 3. * cost * rap * rap * (yeah-1.)
177          !! FACTOR 1/2
178          fshape = 0.5*fshape
179
180      return
181      end
Note: See TracBrowser for help on using the repository browser.