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

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

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