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

Last change on this file since 3325 was 3064, checked in by abierjon, 14 months ago

Mars GCM:
Following r3061 and r3062 :

  • adapt files callphys.def.GCM6 and callphys.def.MCD6 in deftank/ to match with the new orthograph of temp_dependent_m
  • restore the possibility to run with a temperature-dependent contact parameter without using the cloud adaptative subtimestep (MCD6.1 configuration, mteta = linear fit of temp)

AB
File size: 6.9 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_dependent_m) THEN
62c     J.Naar - sep 2023 :
63        if (.not.cloud_adapt_ts) then
64          ! Simple linear parametrisation from Maattaanen 2014
65          ! Maxed out at 0.97 for physical realism
66          ! Used to tune the WC without adaptative ts of microphy (MCD6.1 configuration)
67          mtetalocal = min(0.0044*temp + 0.1831,0.97)
68        else
69          ! With adapt_ts, need to use this relation (same paper, tanh relation):
70          mtetalocal = 0.469+((0.972-0.469)*tanh((temp/158.282)**4.244))
71        endif ! (.not.cloud_adapt_ts)
72      ENDIF ! (temp_dependent_m)
73cccccccccccccccccccccccccccccccccccccccccccccccccc
74ccccccccccc ESSAIS TN MTETA = F (T) cccccccccccccc
75c      if (temp .gt. 200) then
76c         mtetalocal = mtetalocal
77c      else if (temp .lt. 190) then
78c         mtetalocal = mtetalocal-0.05
79c      else
80c         mtetalocal = mtetalocal - (190-temp)*0.005
81c      endif
82c----------------exp law, see Trainer 2008, J. Phys. Chem. C 2009, 113, 2036\u20132040
83       !mtetalocal = max(mtetalocal - 6005*exp(-0.065*temp),0.1)
84       !mtetalocal = max(mtetalocal - 6005*exp(-0.068*temp),0.1)
85               !print*, mtetalocal, temp
86cccccccccccccccccccccccccccccccccccccccccccccccccc
87cccccccccccccccccccccccccccccccccccccccccccccccccc
88      IF (firstcall.and.temp_dependent_m) THEN
89        if (.not.cloud_adapt_ts) then
90          print*, ' ' 
91          print*, 'dear user, please keep in mind that'
92          print*, 'contact parameter IS NOT constant ;'
93          print*, 'Using the following linear fit from'
94          print*, 'Maattanen et al. 2014 (SM linear fit) :'
95          print*, 'm=min(0.0044*temp + 0.1831,0.97)'
96          print*, ' '
97        else !cloud_adapt_ts=.true.
98          print*, ' ' 
99          print*, 'dear user, please keep in mind that'
100          print*, 'contact parameter IS NOT constant ;'
101          print*, 'Using the tanh fit n2 (jsc sample) from'
102          print*, 'Maattanen et al. 2014 (10.1016/j.grj.2014.09.002) :'
103          print*, 'm=0.469+((0.972-0.469)*tanh((temp/158.282)**4.244))'
104          print*, ' '
105        endif
106        firstcall=.false.
107      ELSE IF (firstcall.and.(.not.(temp_dependent_m))) THEN
108        print*, ' ' 
109        print*, 'dear user, please keep in mind that'
110        print*, 'contact parameter IS constant'
111        print*, ' ' 
112        firstcall=.false.
113      END IF
114cccccccccccccccccccccccccccccccccccccccccccccccccc
115cccccccccccccccccccccccccccccccccccccccccccccccccc
116   
117
118      if (sat .gt. 1.) then    ! minimum condition to activate nucleation
119
120        nh2o   = ph2o / kbz / temp
121        rstar  = 2. * sig(temp) * vo1 / (rgp*temp*log(sat))
122        gstar  = 4. * nav * pi * (rstar * rstar * rstar) / (3.*vo1)
123       
124        fshapesimple = (2.+mtetalocal)*(1.-mtetalocal)*(1.-mtetalocal)
125     &                   / 4.
126
127c       Loop over size bins
128        do i=1,nbin_cld
129
130          if ( n_ccn(i) .lt. 1e-10 ) then
131c           no dust, no need to compute nucleation!
132            nucrate(i)=0.
133            ! move on to next bin
134            cycle
135          endif
136
137          if (rad_cld(i).gt.3000.*rstar) then
138            zefshape = fshapesimple
139          else
140            zefshape = fshape(mtetalocal,rad_cld(i)/rstar)
141          endif
142
143          fistar = (4./3.*pi) * sig(temp) * (rstar * rstar) *
144     &             zefshape
145          deltaf = (2.*desorp-surfdif-fistar)/
146     &             (kbz*temp)
147          deltaf = min( max(deltaf, -100.d0), 100.d0)
148
149          if (deltaf.eq.-100.) then
150            nucrate(i) = 0.
151          else
152            nucrate(i)= real(sqrt ( fistar /
153     &               (3.*pi*kbz*temp*(gstar*gstar)) )
154     &                  * kbz * temp * rstar
155     &                  * rstar * 4. * pi
156     &                  * ( nh2o*rad_cld(i) )
157     &                  * ( nh2o*rad_cld(i) )
158     &                  / ( zefshape * nus * m0 )
159     &                  * exp (deltaf))
160          endif
161
162        enddo ! of do i=1,nbin_cld
163
164      else
165
166        do i=1,nbin_cld
167          nucrate(i) = 0.
168        enddo
169
170      endif ! of if (sat .gt. 1.)
171
172      end subroutine nuclea
173
174*********************************************************
175      double precision function fshape(cost,rap)
176      implicit none
177*        function computing the f(m,x) factor           *
178* related to energy required to form a critical embryo  *
179*********************************************************
180
181      double precision, intent(in) :: cost,rap
182      double precision yeah
183
184          !! PHI
185          yeah = sqrt( 1. - 2.*cost*rap + rap*rap )
186          !! FSHAPE = TERM A
187          fshape = (1.-cost*rap) / yeah
188          fshape = fshape * fshape * fshape
189          fshape = 1. + fshape
190          !! ... + TERM B
191          yeah = (rap-cost)/yeah
192          fshape = fshape +
193     & rap*rap*rap*(2.-3.*yeah+yeah*yeah*yeah)
194          !! ... + TERM C
195          fshape = fshape + 3. * cost * rap * rap * (yeah-1.)
196          !! FACTOR 1/2
197          fshape = 0.5*fshape
198
199      end function fshape
200
201
202      END MODULE nuclea_mod
Note: See TracBrowser for help on using the repository browser.