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

Last change on this file since 2613 was 2340, checked in by jvatant, 5 years ago

In addition to r2297, for the n-layer aerosol scheme, enables to set the particle size effective variance with aeronlay_nueff in callphys.def.
--JVO
File size: 12.8 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     
12      real, save ::  rad_h2o
13      real, save ::  rad_h2o_ice
14      real, save ::  Nmix_h2o
15      real, save ::  Nmix_h2o_ice
16!$OMP THREADPRIVATE(rad_h2o,rad_h2o_ice,Nmix_h2o,Nmix_h2o_ice)
17      real, parameter ::  coef_chaud=0.13
18      real, parameter ::  coef_froid=0.09
19
20
21contains
22
23
24!==================================================================
25   subroutine su_aer_radii(ngrid,nlayer,reffrad,nueffrad)
26!==================================================================
27!     Purpose
28!     -------
29!     Compute the effective radii of liquid and icy water particles
30!
31!     Authors
32!     -------
33!     Jeremy Leconte (2012)
34!
35!==================================================================
36      use ioipsl_getin_p_mod, only: getin_p
37      use radinc_h, only: naerkind
38      use aerosol_mod, only: iaero_back2lay, iaero_co2, iaero_dust, &
39                             iaero_h2o, iaero_h2so4, iaero_nh3, iaero_nlay, iaero_aurora
40      use callkeys_mod, only: size_nh3_cloud, nlayaero, aeronlay_size, aeronlay_nueff
41
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, ia   
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!                |-> Done in th n-layer aerosol case (JVO 20)
61
62            if(iaer.eq.iaero_co2)then ! CO2 ice
63               reffrad(1:ngrid,1:nlayer,iaer) = 1.e-4
64               nueffrad(1:ngrid,1:nlayer,iaer) = 0.1
65            endif
66
67            if(iaer.eq.iaero_h2o)then ! H2O ice
68               reffrad(1:ngrid,1:nlayer,iaer) = 1.e-5
69               nueffrad(1:ngrid,1:nlayer,iaer) = 0.1
70            endif
71
72            if(iaer.eq.iaero_dust)then ! dust
73               reffrad(1:ngrid,1:nlayer,iaer) = 1.e-5
74               nueffrad(1:ngrid,1:nlayer,iaer) = 0.1
75            endif
76 
77            if(iaer.eq.iaero_h2so4)then ! H2O ice
78               reffrad(1:ngrid,1:nlayer,iaer) = 1.e-6
79               nueffrad(1:ngrid,1:nlayer,iaer) = 0.1
80            endif
81           
82            if(iaer.eq.iaero_back2lay)then ! Two-layer aerosols
83               reffrad(1:ngrid,1:nlayer,iaer) = 2.e-6
84               nueffrad(1:ngrid,1:nlayer,iaer) = 0.1
85            endif
86
87
88            if(iaer.eq.iaero_nh3)then ! Nh3 cloud
89               reffrad(1:ngrid,1:nlayer,iaer) = size_nh3_cloud
90               nueffrad(1:ngrid,1:nlayer,iaer) = 0.1
91            endif
92
93            do ia=1,nlayaero
94              if(iaer.eq.iaero_nlay(ia))then ! N-layer aerosols
95                 reffrad(1:ngrid,1:nlayer,iaer) = aeronlay_size(ia)
96                 nueffrad(1:ngrid,1:nlayer,iaer) = aeronlay_nueff(ia)
97              endif
98            enddo
99
100            if(iaer.eq.iaero_aurora)then ! Auroral aerosols
101               reffrad(1:ngrid,1:nlayer,iaer) = 3.e-7
102               nueffrad(1:ngrid,1:nlayer,iaer) = 0.1
103            endif
104
105
106            if(iaer.gt.5)then
107               print*,'Error in callcorrk, naerkind is too high (>5).'
108               print*,'The code still needs generalisation to arbitrary'
109               print*,'aerosol kinds and number.'
110               call abort
111            endif
112
113         enddo
114
115
116         if (radfixed) then
117
118            write(*,*)"radius of H2O water particles:"
119            rad_h2o=13. ! default value
120            call getin_p("rad_h2o",rad_h2o)
121            write(*,*)" rad_h2o = ",rad_h2o
122
123            write(*,*)"radius of H2O ice particles:"
124            rad_h2o_ice=35. ! default value
125            call getin_p("rad_h2o_ice",rad_h2o_ice)
126            write(*,*)" rad_h2o_ice = ",rad_h2o_ice
127
128         else
129
130            write(*,*)"Number mixing ratio of H2O water particles:"
131            Nmix_h2o=1.e6 ! default value
132            call getin_p("Nmix_h2o",Nmix_h2o)
133            write(*,*)" Nmix_h2o = ",Nmix_h2o
134
135            write(*,*)"Number mixing ratio of H2O ice particles:"
136            Nmix_h2o_ice=Nmix_h2o ! default value
137            call getin_p("Nmix_h2o_ice",Nmix_h2o_ice)
138            write(*,*)" Nmix_h2o_ice = ",Nmix_h2o_ice
139         endif
140
141      print*,'exit su_aer_radii'
142
143   end subroutine su_aer_radii
144!==================================================================
145
146
147!==================================================================
148   subroutine h2o_reffrad(ngrid,nlayer,pq,pt,reffrad,nueffrad)
149!==================================================================
150!     Purpose
151!     -------
152!     Compute the effective radii of liquid and icy water particles
153!
154!     Authors
155!     -------
156!     Jeremy Leconte (2012)
157!
158!==================================================================
159      use watercommon_h, Only: T_h2O_ice_liq,T_h2O_ice_clouds,rhowater,rhowaterice
160      use comcstfi_mod, only: pi
161      Implicit none
162
163      integer,intent(in) :: ngrid
164      integer,intent(in) :: nlayer
165
166      real, intent(in) :: pq(ngrid,nlayer) !water ice mixing ratios (kg/kg)
167      real, intent(in) :: pt(ngrid,nlayer) !temperature (K)
168      real, intent(out) :: reffrad(ngrid,nlayer)      !aerosol radii
169      real, intent(out) :: nueffrad(ngrid,nlayer) ! dispersion     
170
171      integer :: ig,l
172      real zfice ,zrad,zrad_liq,zrad_ice
173      real,external :: CBRT           
174     
175
176      if (radfixed) then
177         do l=1,nlayer
178            do ig=1,ngrid
179               zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
180               zfice = MIN(MAX(zfice,0.0),1.0)
181               reffrad(ig,l)= rad_h2o * (1.-zfice) + rad_h2o_ice * zfice
182               nueffrad(ig,l) = coef_chaud * (1.-zfice) + coef_froid * zfice
183            enddo
184         enddo
185      else
186         do l=1,nlayer
187            do ig=1,ngrid
188               zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
189               zfice = MIN(MAX(zfice,0.0),1.0)
190               zrad_liq  = CBRT( 3*pq(ig,l)/(4*Nmix_h2o*pi*rhowater) )
191               zrad_ice  = CBRT( 3*pq(ig,l)/(4*Nmix_h2o_ice*pi*rhowaterice) )
192               nueffrad(ig,l) = coef_chaud * (1.-zfice) + coef_froid * zfice
193               zrad = zrad_liq * (1.-zfice) + zrad_ice * zfice
194
195               reffrad(ig,l) = min(max(zrad,1.e-6),1000.e-6)
196               enddo
197            enddo     
198      end if
199
200   end subroutine h2o_reffrad
201!==================================================================
202
203
204!==================================================================
205   subroutine h2o_cloudrad(ngrid,nlayer,pql,reffliq,reffice)
206!==================================================================
207!     Purpose
208!     -------
209!     Compute the effective radii of liquid and icy water particles
210!
211!     Authors
212!     -------
213!     Jeremy Leconte (2012)
214!
215!==================================================================
216      use watercommon_h, Only: rhowater,rhowaterice
217      use comcstfi_mod, only: pi
218      Implicit none
219
220      integer,intent(in) :: ngrid
221      integer,intent(in) :: nlayer
222
223      real, intent(in) :: pql(ngrid,nlayer) !condensed water mixing ratios (kg/kg)
224      real, intent(out) :: reffliq(ngrid,nlayer),reffice(ngrid,nlayer)     !liquid and ice water particle radii (m)
225
226      real,external :: CBRT           
227      integer :: i,k
228
229      if (radfixed) then
230         reffliq(1:ngrid,1:nlayer)= rad_h2o
231         reffice(1:ngrid,1:nlayer)= rad_h2o_ice
232      else
233         do k=1,nlayer
234           do i=1,ngrid
235             reffliq(i,k) = CBRT(3*pql(i,k)/(4*Nmix_h2o*pi*rhowater))
236             reffliq(i,k) = min(max(reffliq(i,k),1.e-6),1000.e-6)
237           
238             reffice(i,k) = CBRT(3*pql(i,k)/(4*Nmix_h2o_ice*pi*rhowaterice))
239             reffice(i,k) = min(max(reffice(i,k),1.e-6),1000.e-6)
240           enddo
241         enddo
242      endif
243
244   end subroutine h2o_cloudrad
245!==================================================================
246
247
248
249!==================================================================
250   subroutine co2_reffrad(ngrid,nlayer,nq,pq,reffrad)
251!==================================================================
252!     Purpose
253!     -------
254!     Compute the effective radii of co2 ice particles
255!
256!     Authors
257!     -------
258!     Jeremy Leconte (2012)
259!
260!==================================================================
261      USE tracer_h, only:igcm_co2_ice,rho_co2
262      use comcstfi_mod, only: pi
263      Implicit none
264
265      integer,intent(in) :: ngrid,nlayer,nq
266
267      real, intent(in) :: pq(ngrid,nlayer,nq) !tracer mixing ratios (kg/kg)
268      real, intent(out) :: reffrad(ngrid,nlayer)      !co2 ice particles radii (m)
269
270      integer :: ig,l
271      real :: zrad   
272      real,external :: CBRT           
273           
274     
275
276      if (radfixed) then
277         reffrad(1:ngrid,1:nlayer) = 5.e-5 ! CO2 ice
278      else
279         do l=1,nlayer
280            do ig=1,ngrid
281               zrad = CBRT( 3*pq(ig,l,igcm_co2_ice)/(4*Nmix_co2*pi*rho_co2) )
282               reffrad(ig,l) = min(max(zrad,1.e-6),100.e-6)
283            enddo
284         enddo     
285      end if
286
287   end subroutine co2_reffrad
288!==================================================================
289
290
291
292!==================================================================
293   subroutine dust_reffrad(ngrid,nlayer,reffrad)
294!==================================================================
295!     Purpose
296!     -------
297!     Compute the effective radii of dust particles
298!
299!     Authors
300!     -------
301!     Jeremy Leconte (2012)
302!
303!==================================================================
304      Implicit none
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      integer,intent(in) :: ngrid
332      integer,intent(in) :: nlayer
333
334      real, intent(out) :: reffrad(ngrid,nlayer)      !h2so4 particle radii (m)
335               
336      reffrad(1:ngrid,1:nlayer) = 1.e-6 ! h2so4
337
338   end subroutine h2so4_reffrad
339!==================================================================
340
341!==================================================================
342   subroutine back2lay_reffrad(ngrid,reffrad,nlayer,pplev)
343!==================================================================
344!     Purpose
345!     -------
346!     Compute the effective radii of particles in a 2-layer model
347!
348!     Authors
349!     -------
350!     Sandrine Guerlet (2013)
351!
352!==================================================================
353      use callkeys_mod, only: pres_bottom_tropo,pres_top_tropo,size_tropo,  &
354                              pres_bottom_strato,size_strato
355 
356      Implicit none
357
358      integer,intent(in) :: ngrid
359
360      real, intent(out) :: reffrad(ngrid,nlayer)      ! particle radii (m)
361      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
362      INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
363      REAL :: expfactor
364      INTEGER l,ig
365           
366      reffrad(:,:)=1e-6  !!initialization, not important
367          DO ig=1,ngrid
368            DO l=1,nlayer-1
369              IF (pplev(ig,l) .le. pres_bottom_tropo .and. pplev(ig,l) .ge. pres_top_tropo) THEN
370                reffrad(ig,l) = size_tropo
371              ELSEIF (pplev(ig,l) .lt. pres_top_tropo .and. pplev(ig,l) .gt. pres_bottom_strato) THEN
372                expfactor=log(size_strato/size_tropo) / log(pres_bottom_strato/pres_top_tropo)
373                reffrad(ig,l)= size_tropo*((pplev(ig,l)/pres_top_tropo)**expfactor)
374              ELSEIF (pplev(ig,l) .le. pres_bottom_strato) then
375                reffrad(ig,l) = size_strato
376              ENDIF
377            ENDDO
378          ENDDO
379
380   end subroutine back2lay_reffrad
381!==================================================================
382
383
384end module radii_mod
385!==================================================================
Note: See TracBrowser for help on using the repository browser.