source: trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90 @ 1255

Last change on this file since 1255 was 1026, checked in by sglmd, 11 years ago

Added a flexible, 2-layer aerosol scenario (initially for Saturn, but can be simply tuned). Called aeroback2lay.
File size: 11.9 KB
Line 
1!==================================================================
2module radii_mod
3!==================================================================
4!  module to centralize the radii calculations for aerosols
5! OK for water but should be extended to other aerosols (CO2,...)
6!==================================================================
7     
8!     water cloud optical properties
9     
10      real, save ::  rad_h2o
11      real, save ::  rad_h2o_ice
12      real, save ::  Nmix_h2o
13      real, save ::  Nmix_h2o_ice
14      real, parameter ::  coef_chaud=0.13
15      real, parameter ::  coef_froid=0.09
16
17
18      contains
19
20
21!==================================================================
22   subroutine su_aer_radii(ngrid,reffrad,nueffrad)
23!==================================================================
24!     Purpose
25!     -------
26!     Compute the effective radii of liquid and icy water particles
27!
28!     Authors
29!     -------
30!     Jeremy Leconte (2012)
31!
32!==================================================================
33 ! to use  'getin'
34      use ioipsl_getincom
35      use radinc_h, only: naerkind
36      use aerosol_mod
37!      USE tracer_h
38      Implicit none
39
40#include "callkeys.h"
41#include "dimensions.h"
42#include "dimphys.h"
43
44      integer,intent(in) :: ngrid
45
46      real, intent(out) :: reffrad(ngrid,nlayermx,naerkind)      !aerosols radii (K)
47      real, intent(out) :: nueffrad(ngrid,nlayermx,naerkind)     !variance     
48
49      logical, save :: firstcall=.true.
50      integer :: iaer   
51     
52      print*,'enter su_aer_radii'
53          do iaer=1,naerkind
54!     these values will change once the microphysics gets to work
55!     UNLESS tracer=.false., in which case we should be working with
56!     a fixed aerosol layer, and be able to define reffrad in a
57!     .def file. To be improved!
58
59            if(iaer.eq.iaero_co2)then ! CO2 ice
60               reffrad(1:ngrid,1:nlayermx,iaer) = 1.e-4
61               nueffrad(1:ngrid,1:nlayermx,iaer) = 0.1
62            endif
63
64            if(iaer.eq.iaero_h2o)then ! H2O ice
65               reffrad(1:ngrid,1:nlayermx,iaer) = 1.e-5
66               nueffrad(1:ngrid,1:nlayermx,iaer) = 0.1
67            endif
68
69            if(iaer.eq.iaero_dust)then ! dust
70               reffrad(1:ngrid,1:nlayermx,iaer) = 1.e-5
71               nueffrad(1:ngrid,1:nlayermx,iaer) = 0.1
72            endif
73 
74            if(iaer.eq.iaero_h2so4)then ! H2O ice
75               reffrad(1:ngrid,1:nlayermx,iaer) = 1.e-6
76               nueffrad(1:ngrid,1:nlayermx,iaer) = 0.1
77            endif
78           
79            if(iaer.eq.iaero_back2lay)then ! Two-layer aerosols
80               reffrad(1:ngrid,1:nlayermx,iaer) = 2.e-6
81               nueffrad(1:ngrid,1:nlayermx,iaer) = 0.1
82            endif
83
84
85
86            if(iaer.gt.5)then
87               print*,'Error in callcorrk, naerkind is too high (>5).'
88               print*,'The code still needs generalisation to arbitrary'
89               print*,'aerosol kinds and number.'
90               call abort
91            endif
92
93         enddo
94
95
96         if (radfixed) then
97
98            write(*,*)"radius of H2O water particles:"
99            rad_h2o=13. ! default value
100            call getin("rad_h2o",rad_h2o)
101            write(*,*)" rad_h2o = ",rad_h2o
102
103            write(*,*)"radius of H2O ice particles:"
104            rad_h2o_ice=35. ! default value
105            call getin("rad_h2o_ice",rad_h2o_ice)
106            write(*,*)" rad_h2o_ice = ",rad_h2o_ice
107
108         else
109
110            write(*,*)"Number mixing ratio of H2O water particles:"
111            Nmix_h2o=1.e6 ! default value
112            call getin("Nmix_h2o",Nmix_h2o)
113            write(*,*)" Nmix_h2o = ",Nmix_h2o
114
115            write(*,*)"Number mixing ratio of H2O ice particles:"
116            Nmix_h2o_ice=Nmix_h2o ! default value
117            call getin("Nmix_h2o_ice",Nmix_h2o_ice)
118            write(*,*)" Nmix_h2o_ice = ",Nmix_h2o_ice
119         endif
120
121      print*,'exit su_aer_radii'
122
123   end subroutine su_aer_radii
124!==================================================================
125
126
127!==================================================================
128   subroutine h2o_reffrad(ngrid,pq,pt,reffrad,nueffrad)
129!==================================================================
130!     Purpose
131!     -------
132!     Compute the effective radii of liquid and icy water particles
133!
134!     Authors
135!     -------
136!     Jeremy Leconte (2012)
137!
138!==================================================================
139      use watercommon_h, Only: T_h2O_ice_liq,T_h2O_ice_clouds,rhowater,rhowaterice
140      Implicit none
141
142#include "callkeys.h"
143#include "dimensions.h"
144#include "dimphys.h"
145#include "comcstfi.h"
146
147      integer,intent(in) :: ngrid
148
149      real, intent(in) :: pq(ngrid,nlayermx) !water ice mixing ratios (kg/kg)
150      real, intent(in) :: pt(ngrid,nlayermx) !temperature (K)
151      real, intent(out) :: reffrad(ngrid,nlayermx)      !aerosol radii
152      real, intent(out) :: nueffrad(ngrid,nlayermx) ! dispersion     
153
154      integer :: ig,l
155      real zfice ,zrad,zrad_liq,zrad_ice
156      real,external :: CBRT           
157     
158
159      if (radfixed) then
160         do l=1,nlayermx
161            do ig=1,ngrid
162               zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
163               zfice = MIN(MAX(zfice,0.0),1.0)
164               reffrad(ig,l)= rad_h2o * (1.-zfice) + rad_h2o_ice * zfice
165               nueffrad(ig,l) = coef_chaud * (1.-zfice) + coef_froid * zfice
166            enddo
167         enddo
168      else
169         do l=1,nlayermx
170            do ig=1,ngrid
171               zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
172               zfice = MIN(MAX(zfice,0.0),1.0)
173               zrad_liq  = CBRT( 3*pq(ig,l)/(4*Nmix_h2o*pi*rhowater) )
174               zrad_ice  = CBRT( 3*pq(ig,l)/(4*Nmix_h2o_ice*pi*rhowaterice) )
175               nueffrad(ig,l) = coef_chaud * (1.-zfice) + coef_froid * zfice
176               zrad = zrad_liq * (1.-zfice) + zrad_ice * zfice
177
178               reffrad(ig,l) = min(max(zrad,1.e-6),1000.e-6)
179               enddo
180            enddo     
181      end if
182
183   end subroutine h2o_reffrad
184!==================================================================
185
186
187!==================================================================
188   subroutine h2o_cloudrad(ngrid,pql,reffliq,reffice)
189!==================================================================
190!     Purpose
191!     -------
192!     Compute the effective radii of liquid and icy water particles
193!
194!     Authors
195!     -------
196!     Jeremy Leconte (2012)
197!
198!==================================================================
199      use watercommon_h, Only: rhowater,rhowaterice
200      Implicit none
201
202#include "callkeys.h"
203#include "dimensions.h"
204#include "dimphys.h"
205#include "comcstfi.h"
206
207      integer,intent(in) :: ngrid
208
209      real, intent(in) :: pql(ngrid,nlayermx) !condensed water mixing ratios (kg/kg)
210      real, intent(out) :: reffliq(ngrid,nlayermx),reffice(ngrid,nlayermx)     !liquid and ice water particle radii (K)
211
212      real,external :: CBRT           
213     
214
215      if (radfixed) then
216         reffliq(1:ngrid,1:nlayermx)= rad_h2o
217         reffice(1:ngrid,1:nlayermx)= rad_h2o_ice
218      else
219         reffliq(1:ngrid,1:nlayermx)  = CBRT( 3*pql(1:ngrid,1:nlayermx)/(4*Nmix_h2o*pi*rhowater) )
220         reffliq(1:ngrid,1:nlayermx)  = min(max(reffliq(1:ngrid,1:nlayermx),1.e-6),1000.e-6)
221         reffice(1:ngrid,1:nlayermx)  = CBRT( 3*pql(1:ngrid,1:nlayermx)/(4*Nmix_h2o_ice*pi*rhowaterice) )
222         reffice(1:ngrid,1:nlayermx)  = min(max(reffice(1:ngrid,1:nlayermx),1.e-6),1000.e-6)
223      end if
224
225   end subroutine h2o_cloudrad
226!==================================================================
227
228
229
230!==================================================================
231   subroutine co2_reffrad(ngrid,nq,pq,reffrad)
232!==================================================================
233!     Purpose
234!     -------
235!     Compute the effective radii of co2 ice particles
236!
237!     Authors
238!     -------
239!     Jeremy Leconte (2012)
240!
241!==================================================================
242      USE tracer_h, only:igcm_co2_ice,rho_co2
243      Implicit none
244
245#include "callkeys.h"
246#include "dimensions.h"
247#include "dimphys.h"
248#include "comcstfi.h"
249
250      integer,intent(in) :: ngrid,nq
251
252      real, intent(in) :: pq(ngrid,nlayermx,nq) !tracer mixing ratios (kg/kg)
253      real, intent(out) :: reffrad(ngrid,nlayermx)      !co2 ice particles radii (K)
254
255      integer :: ig,l
256      real :: zrad   
257      real,external :: CBRT           
258           
259     
260
261      if (radfixed) then
262         reffrad(1:ngrid,1:nlayermx) = 5.e-5 ! CO2 ice
263      else
264         do l=1,nlayermx
265            do ig=1,ngrid
266               zrad = CBRT( 3*pq(ig,l,igcm_co2_ice)/(4*Nmix_co2*pi*rho_co2) )
267               reffrad(ig,l) = min(max(zrad,1.e-6),100.e-6)
268            enddo
269         enddo     
270      end if
271
272   end subroutine co2_reffrad
273!==================================================================
274
275
276
277!==================================================================
278   subroutine dust_reffrad(ngrid,reffrad)
279!==================================================================
280!     Purpose
281!     -------
282!     Compute the effective radii of dust particles
283!
284!     Authors
285!     -------
286!     Jeremy Leconte (2012)
287!
288!==================================================================
289      Implicit none
290
291#include "callkeys.h"
292#include "dimensions.h"
293#include "dimphys.h"
294
295      integer,intent(in) :: ngrid
296
297      real, intent(out) :: reffrad(ngrid,nlayermx)      !dust particles radii (K)
298           
299      reffrad(1:ngrid,1:nlayermx) = 2.e-6 ! dust
300
301   end subroutine dust_reffrad
302!==================================================================
303
304
305!==================================================================
306   subroutine h2so4_reffrad(ngrid,reffrad)
307!==================================================================
308!     Purpose
309!     -------
310!     Compute the effective radii of h2so4 particles
311!
312!     Authors
313!     -------
314!     Jeremy Leconte (2012)
315!
316!==================================================================
317      Implicit none
318
319#include "callkeys.h"
320#include "dimensions.h"
321#include "dimphys.h"
322
323      integer,intent(in) :: ngrid
324
325      real, intent(out) :: reffrad(ngrid,nlayermx)      !h2so4 particle radii (K)
326               
327      reffrad(1:ngrid,1:nlayermx) = 1.e-6 ! h2so4
328
329   end subroutine h2so4_reffrad
330!==================================================================
331
332!==================================================================
333   subroutine back2lay_reffrad(ngrid,reffrad,nlayer,pplev)
334!==================================================================
335!     Purpose
336!     -------
337!     Compute the effective radii of particles in a 2-layer model
338!
339!     Authors
340!     -------
341!     Sandrine Guerlet (2013)
342!
343!==================================================================
344 
345      use aerosol_mod   !! Particle sizes and boundaries of aerosol layers defined there
346     Implicit none
347
348#include "callkeys.h"
349#include "dimensions.h"
350#include "dimphys.h"
351
352      integer,intent(in) :: ngrid
353
354      real, intent(out) :: reffrad(ngrid,nlayermx)      ! particle radii
355      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
356      INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
357      REAL :: expfactor
358      INTEGER l,ig
359           
360      reffrad(:,:)=1e-6  !!initialization, not important
361          DO ig=1,ngrid
362            DO l=1,nlayer-1
363              IF (pplev(ig,l) .le. pres_bottom_tropo .and. pplev(ig,l) .ge. pres_top_tropo) THEN
364                reffrad(ig,l) = size_tropo
365              ELSEIF (pplev(ig,l) .lt. pres_top_tropo .and. pplev(ig,l) .gt. pres_bottom_strato) THEN
366                expfactor=log(size_strato/size_tropo) / log(pres_bottom_strato/pres_top_tropo)
367                reffrad(ig,l)= size_tropo*((pplev(ig,l)/pres_top_tropo)**expfactor)
368              ELSEIF (pplev(ig,l) .le. pres_bottom_strato) then
369                reffrad(ig,l) = size_strato
370              ENDIF
371            ENDDO
372          ENDDO
373
374   end subroutine back2lay_reffrad
375!==================================================================
376
377
378
379end module radii_mod
380!==================================================================
Note: See TracBrowser for help on using the repository browser.