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

Last change on this file since 1477 was 1397, checked in by milmd, 10 years ago

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