source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/radii_mod.F90 @ 227

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

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