source: trunk/LMDZ.MARS/libf/phymars/nucleaco2.F @ 2156

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

changed upper-case CO2 in lower-case co2 to be consistent everywhere. also Fortran is not case-sensitive so using a mix of upper case and lower case is not usually good practice.
  • Property svn:executable set to *
File size: 6.3 KB
Line 
1*******************************************************
2*                                                     *
3      subroutine nucleaco2(pco2,temp,sat,n_ccn,nucrate,
4     &           n_ccn_h2oice,rad_h2oice,nucrate_h2oice,
5     &           vo2co2)
6      USE comcstfi_h
7
8      implicit none
9*                                                     *
10*   This subroutine computes the nucleation rate      *
11*   as given in Pruppacher & Klett (1978) in the      *
12*   case of water ice forming on a solid substrate.   *
13*     Definition refined by Keese (jgr,1989)          *
14*   Authors: F. Montmessin                            *
15*     Adapted for the LMD/GCM by J.-B. Madeleine      *
16*     (October 2011)                                  *
17*     Optimisation by A. Spiga (February 2012)        *
18*   CO2 nucleation routine dev. by Constantino        *
19*     Listowski and Joachim Audouard (2016-2017),     *
20*     adapted from the water ice nucleation     
21* It computes two different nucleation rates : one
22* on the dust CCN distribution and the other one on
23* the water ice particles distribution
24*******************************************************
25 ! nucrate = output
26 ! nucrate_h2o en sortie aussi :
27!nucleation sur dust et h2o separement ici
28
29      include "microphys.h"
30      include "callkeys.h"
31
32c     Inputs
33      DOUBLE PRECISION pco2,sat,vo2co2
34      DOUBLE PRECISION n_ccn(nbinco2_cld), n_ccn_h2oice(nbinco2_cld)
35      REAL temp !temperature
36c     Output
37      DOUBLE PRECISION nucrate(nbinco2_cld)
38      DOUBLE PRECISION nucrate_h2oice(nbinco2_cld) ! h2o as substrate
39      double precision rad_h2oice(nbinco2_cld) ! h2o ice grid (as substrate)
40
41c     Local variables
42      DOUBLE PRECISION nco2
43      DOUBLE PRECISION rstar    ! Radius of the critical germ (m)
44      DOUBLE PRECISION gstar    ! # of molecules forming a critical embryo
45      DOUBLE PRECISION fistar   ! Activation energy required to form a critical embryo (J)
46      DOUBLE PRECISION fshapeco2   ! function defined at the end of the file
47      DOUBLE PRECISION deltaf     
48      double precision mtetalocal,mtetalocalh ! local mteta in double precision
49      double precision fshapeco2simple,zefshapeco2
50      integer i
51c     *************************************************
52
53      mtetalocal = dble(mtetaco2)  !! use mtetalocal for better performance
54      mtetalocalh=dble(mteta)
55
56
57      IF (sat .gt. 1.) THEN    ! minimum condition to activate nucleation
58
59        nco2   = pco2 / kbz / temp
60        rstar  = 2. * sigco2 * vo2co2 / (kbz*temp*dlog(sat))
61        gstar  = 4. * pi * (rstar * rstar * rstar) / (3.*vo2co2)
62       
63       fshapeco2simple = (2.+mtetalocal)*(1.-mtetalocal)*(1.-mtetalocal)
64     &                   / 4.
65
66c       Loop over size bins
67        do i=1,nbinco2_cld
68c            write(*,*) "IN NUCLEA, i, RAD_CLDCO2(i) = ",i, rad_cldco2(i),
69c     &          n_ccn(i)
70
71          if ( n_ccn(i) .lt. 1e-10 ) then
72c           no dust, no need to compute nucleation!
73            nucrate(i)=0.
74c            goto 210
75c          endif
76          else
77            if (rad_cldco2(i).gt.3000.*rstar) then
78              zefshapeco2 = fshapeco2simple
79            else
80             zefshapeco2 = fshapeco2(mtetalocal,rad_cldco2(i)/rstar)
81            endif
82
83            fistar = (4./3.*pi) * sigco2 * (rstar * rstar) *
84     &             zefshapeco2
85            deltaf = (2.*desorpco2-surfdifco2-fistar)/
86     &             (kbz*temp)
87            deltaf = min( max(deltaf, -100.d0), 100.d0)
88
89            if (deltaf.eq.-100.) then
90                nucrate(i) = 0.
91            else
92                nucrate(i)= dble(sqrt ( fistar /
93     &               (3.*pi*kbz*temp*(gstar*gstar)) )
94     &                  * kbz * temp * rstar
95     &                  * rstar * 4. * pi
96     &                  * ( nco2*rad_cldco2(i) )
97     &                  * ( nco2*rad_cldco2(i) )
98     &                  / ( zefshapeco2 * nusco2 * m0co2 )
99     &                  * dexp (deltaf))
100
101           
102            endif
103          endif ! if n_ccn(i) .lt. 1e-10
104
105          if (co2useh2o) then
106
107            if ( n_ccn_h2oice(i) .lt. 1e-10 ) then
108c           no H2O ice, no need to compute nucleation!
109               nucrate_h2oice(i)=0.
110            else   
111               if (rad_h2oice(i).gt.3000.*rstar) then
112                 zefshapeco2 = (2.+mtetalocalh)*(1.-mtetalocalh)*
113     &                (1.-mtetalocalh) / 4.
114               else  ! same m for dust/h2o ice
115                 zefshapeco2 = fshapeco2(mtetalocalh,
116     &                            (rad_h2oice(i)/rstar))
117               endif
118
119               fistar = (4./3.*pi) * sigco2 * (rstar * rstar) *
120     &             zefshapeco2
121              deltaf = (2.*desorpco2-surfdifco2-fistar)/
122     &             (kbz*temp)
123              deltaf = min( max(deltaf, -100.d0), 100.d0)
124
125              if (deltaf.eq.-100.) then
126                  nucrate_h2oice(i) = 0.
127              else
128                  nucrate_h2oice(i)= dble(sqrt ( fistar /
129     &               (3.*pi*kbz*temp*(gstar*gstar)) )
130     &                  * kbz * temp * rstar
131     &                  * rstar * 4. * pi
132     &                  * ( nco2*rad_h2oice(i) )
133     &                  * ( nco2*rad_h2oice(i) )
134     &                  / ( zefshapeco2 * nusco2 * m0co2 )
135     &                  * dexp (deltaf))
136              endif
137            endif
138          endif
139        enddo
140
141      ELSE ! parcelle d'air non saturée
142
143        do i=1,nbinco2_cld
144          nucrate(i) = 0.
145          nucrate_h2oice(i) = 0.
146        enddo
147
148      ENDIF ! if (sat .gt. 1.)
149
150      end
151
152*********************************************************
153      double precision function fshapeco2(cost,rap)
154      implicit none
155*        function computing the f(m,x) factor           *
156* related to energy required to form a critical embryo  *
157*********************************************************
158
159      double precision cost,rap
160      double precision yeah
161
162          !! PHI
163          yeah = sqrt( 1. - 2.*cost*rap + rap*rap )
164          !! FSHAPECO2 = TERM A
165          fshapeco2 = (1.-cost*rap) / yeah
166          fshapeco2 = fshapeco2 * fshapeco2 * fshapeco2
167          fshapeco2 = 1. + fshapeco2
168          !! ... + TERM B
169          yeah = (rap-cost)/yeah
170          fshapeco2 = fshapeco2 +
171     & rap*rap*rap*(2.-3.*yeah+yeah*yeah*yeah)
172          !! ... + TERM C
173          fshapeco2 = fshapeco2 + 3. * cost * rap * rap * (yeah-1.)
174          !! FACTOR 1/2
175          fshapeco2 = 0.5*fshapeco2
176
177      end
Note: See TracBrowser for help on using the repository browser.