source: trunk/LMDZ.MARS/libf/phymars/updaterad.F90 @ 740

Last change on this file since 740 was 740, checked in by tnavarro, 12 years ago

module for ice and radius radius computation

File size: 7.8 KB
Line 
1module updaterad
2
3
4
5! This module intents to group together all ice and dust radius computations done in the GCM,
6! so that it stays coherent through the code and make (numerous) radius bugs easier to debug.
7! All different thresholds values are defined below.
8
9! Radius computation do not always occur on the whole grid (cf. improvedcloud).
10! So, subroutines are designed for scalar values instead of tables
11
12
13! T. Navarro, June 2012
14
15
16
17! For instance, if R^3 is lower than r3icemin, then R is set to ricemin.
18! So, ricemin should be the cubic root of r3icemin, but not necessarily ...
19
20real, parameter :: r3icemin = 1.e-30   ! ie ricemin  = 0.0001 microns
21real, parameter :: ricemin  = 1.e-10
22
23real, parameter :: r3icemax = 125.e-12 ! ie ricemax  = 500 microns
24real, parameter :: ricemax  = 500.e-6
25
26real, parameter :: qice_threshold  = 1.e-15 ! 1.e-10
27
28
29
30real, parameter :: nccn_threshold  =  1.
31real, parameter :: qccn_threshold  =  1.e-20
32
33real, parameter :: r3ccnmin = 1.e-21    ! ie rccnmin = 0.1 microns
34real, parameter :: rccnmin  = 0.1e-6
35
36real, parameter :: r3ccnmax = 125.e-12  ! ie rccnmax  = 500 microns
37real, parameter :: rccnmax  = 500.e-6 
38
39
40
41real, parameter :: ndust_threshold  = 1.
42real, parameter :: qdust_threshold  = 1.e-20
43
44real, parameter :: r3dustmin = 1.e-21  ! ie rdustmin = 0.1 microns
45real, parameter :: rdustmin  = 0.1e-6
46
47real, parameter :: r3dustmax = 125.e-12 ! ie rdustmax  = 500 microns
48real, parameter :: rdustmax  = 500.e-6 
49
50
51real, parameter :: rdust0 = 0.8e-7
52
53
54     
55
56contains
57
58
59!============================================================================
60!============================================================================
61!============================================================================
62! Update ice radius if microphys == true
63subroutine updaterice_micro(qice,qccn,nccn,coeff,rice,rhocloud)
64implicit none
65
66#include "dimensions.h"
67#include "dimphys.h"
68#include "comcstfi.h"
69#include "tracer.h"
70
71real, intent(in)  :: qice,qccn,nccn
72real, intent(in)  :: coeff         ! this coeff is tauscaling if microphy = T (possibly ccn_factor^-1 otherwise)
73real, intent(out) :: rice,rhocloud ! rhocloud is needed for sedimentation and is also a good diagnostic variable
74real nccn_true,qccn_true ! radius of the ice crystal core to the power of 3
75
76   
77nccn_true = max(nccn * coeff, 1.e-30)
78qccn_true = max(qccn * coeff, 1.e-30)
79
80
81!! Nota: It is very dangerous to apply a threshold on qccn or nccn to force rice to be ricemin.
82!! Indeed, one can obtain ricemin for small but non negligible qice values, and therefore hugely opaque clouds.
83 
84if (qice .ge. qice_threshold) then
85
86  rhocloud = (qice*rho_ice + qccn_true*rho_dust) / (qice + qccn_true)
87  rhocloud = min(max(rhocloud,rho_ice),rho_dust)
88
89  rice = (qice + qccn_true) * 0.75 / pi / rhocloud / nccn_true
90 
91  if (rice .le. r3icemin) then
92    rice = ricemin
93  else if (rice .ge. r3icemax) then
94    rice = ricemax
95  else
96    rice = rice**(1./3.) ! here rice is always positive
97  endif
98
99else
100
101  rice = ricemin
102  rhocloud = rho_dust
103
104endif
105
106end subroutine updaterice_micro
107!============================================================================
108!============================================================================
109!============================================================================
110
111
112
113
114
115!============================================================================
116!============================================================================
117!============================================================================
118! Update ice radius from a typical profile if microphys == false
119subroutine updaterice_typ(qice,tau,pzlay,rice)
120implicit none
121
122#include "dimensions.h"
123#include "dimphys.h"
124#include "comcstfi.h"
125#include "tracer.h"
126
127real, intent(in)  :: qice
128real, intent(in)  :: tau   ! tau for dust
129real, intent(in)  :: pzlay ! altitude at the middle of the layers
130real, intent(out) :: rice
131real rccn,nccn ! radius  and number of ice crystals
132
133!   Typical CCN profile following Montmessin et al. 2004
134!   (N0=2e6 m-3 has been converted into N0=1.3e8 kg-1, otherwise the equation for rice is not homogeneous...)
135     nccn  = 1.3e+8*max(tau,0.001)/0.1*exp(-pzlay/10000.)
136!   The previously used profile was not correct:
137!   Nccn=( epaisseur/masse ) * 2.e+6/0.1*max(tau,0.001)*exp(-pzlay/10000.)
138     
139if (nccn .le. 1) then
140
141  rice = ricemin
142
143else
144
145! Typical dust radius profile:
146  rccn = max(rdust0*exp(-pzlay/18000.),1.e-9)
147  rice  = qice * 0.75 / pi / rho_ice / nccn + rccn*rccn*rccn
148       
149  if (rice .le. r3icemin) then
150    rice = ricemin
151  else if (rice .ge. r3icemax) then
152    rice = ricemax
153  else
154    rice = rice**(1./3.) ! here rice is always positive
155  endif
156
157endif
158 
159end subroutine updaterice_typ
160!============================================================================
161!============================================================================
162!============================================================================
163
164
165
166
167
168
169!============================================================================
170!============================================================================
171!============================================================================
172! This subroutine computes the geometric mean radius(or number median radius)
173! For a lognormal distribution :
174! geometric mean radius = mass mean radius x exp(-1.5 sigma0^2)
175! To be used with doubleq == true. otherwise, rdust is constant !!!
176subroutine updaterdust(qdust,ndust,rdust,tauscaling)
177implicit none
178
179#include "dimensions.h"
180#include "dimphys.h"
181#include "comcstfi.h"
182#include "tracer.h"
183
184real, intent(in) :: qdust,ndust ! needed if doubleq
185real, intent(in), optional :: tauscaling ! useful for realistic thresholds
186real, intent(out) :: rdust
187real coeff
188
189if (present(tauscaling)) then
190  coeff = tauscaling ! thresholds on realistic values
191else
192  coeff = 1. ! thresholds on virtual values
193endif
194
195if ((ndust .le. ndust_threshold/coeff) .or. (qdust .le. qdust_threshold/coeff)) then
196
197   rdust = rdustmin
198
199else
200
201   rdust  = r3n_q * qdust / ndust
202
203   if (rdust .le. r3dustmin) then
204       rdust = rdustmin
205   else if (rdust .ge. r3dustmax) then
206       rdust = rdustmax
207   else
208      rdust = rdust**(1./3.)
209   endif
210
211endif
212     
213
214end subroutine updaterdust
215!============================================================================
216!============================================================================
217!============================================================================
218
219
220
221
222
223
224!============================================================================
225!============================================================================
226!============================================================================
227! This subroutine computes the mass mean radius,
228! used for heterogenous nucleation on CCNs in microphysics.
229! For a lognormal distribution :
230! geometric mean radius = mass mean radius x exp(-1.5 sigma0^2)
231subroutine updaterccn(qccn,nccn,rccn,tauscaling)
232implicit none
233
234#include "dimensions.h"
235#include "dimphys.h"
236#include "comcstfi.h"
237#include "tracer.h"
238
239real, intent(in) :: qccn,nccn ! needed if doubleq
240real, intent(in), optional :: tauscaling ! useful for realistic thresholds
241
242real, intent(out) :: rccn
243
244real coeff
245
246
247if (present(tauscaling)) then
248  coeff = tauscaling ! threshold on realistic values
249else
250  coeff = 1. ! threshold on virtual values
251endif
252
253if ((nccn .le. nccn_threshold/coeff) .or. (qccn .le. qccn_threshold/coeff)) then
254
255   rccn = rccnmin
256
257else
258
259   rccn  =  qccn * 0.75 / pi / rho_dust / nccn
260
261   if (rccn .le. r3ccnmin) then
262       rccn = rccnmin
263   else if (rccn .ge. r3ccnmax) then
264       rccn = rccnmax
265   else
266       rccn = rccn**(1./3.)
267   endif
268
269endif
270     
271end subroutine updaterccn
272!============================================================================
273!============================================================================
274!============================================================================
275
276
277
278end module updaterad
279
Note: See TracBrowser for help on using the repository browser.