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

Last change on this file since 3026 was 3008, checked in by emillour, 18 months ago

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