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

Last change on this file since 1974 was 1677, checked in by sglmd, 8 years ago

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