source: trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90 @ 726

Last change on this file since 726 was 726, checked in by jleconte, 13 years ago

17/07/2012 == JL for LK

  • Generalization of aerosol scheme:
    • any number of aerosols can be used and id numbers are determined consistently by the code. Aerosol order not important anymore.
    • addition of a module with the id numbers for aerosols (aerosol_mod.F90).
    • initialization of aerosols id numbers in iniaerosol.F90
    • compile with -s x where x *must* be equal to the number of aerosols turned on in callphys.def (either by a flag or by dusttau>0 for dust). => may have to erase object files when compiling with s option for the first time.
  • For no aerosols, run with aeroco2=.true. and aerofixco2=.true (the default distribution for fixed co2

aerosols is 1.e-9; can be changed in aeropacity).

  • If starting from an old start file, recreate start file with the q=0 option in newstart.e.
  • update callphys.def with aeroXXX and aerofixXXX options (only XXX=co2,h2o supported for

now). Dust is activated by setting dusttau>0. See the early mars case in deftank.

  • To add other aerosols, see Laura Kerber.
File size: 13.8 KB
Line 
1      Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pq, &
2         aerosol,reffrad,QREFvis3d,QREFir3d,tau_col,cloudfrac,totcloudfrac,clearsky)
3
4       use radinc_h, only : L_TAUMAX,naerkind
5       use aerosol_mod
6                 
7       implicit none
8
9!==================================================================
10!     
11!     Purpose
12!     -------
13!     Compute aerosol optical depth in each gridbox.
14!     
15!     Authors
16!     -------
17!     F. Forget
18!     F. Montmessin (water ice scheme)
19!     update J.-B. Madeleine (2008)
20!     dust removal, simplification by Robin Wordsworth (2009)
21!
22!     Input
23!     -----
24!     ngrid             Number of horizontal gridpoints
25!     nlayer            Number of layers
26!     nq                Number of tracers
27!     pplev             Pressure (Pa) at each layer boundary
28!     pq                Aerosol mixing ratio
29!     reffrad(ngrid,nlayer,naerkind)         Aerosol effective radius
30!     QREFvis3d(ngridmx,nlayermx,naerkind) \ 3d extinction coefficients
31!     QREFir3d(ngridmx,nlayermx,naerkind)  / at reference wavelengths
32!
33!     Output
34!     ------
35!     aerosol            Aerosol optical depth in layer l, grid point ig
36!     tau_col            Total column optical depth at grid point ig
37!
38!=======================================================================
39
40#include "dimensions.h"
41#include "dimphys.h"
42#include "callkeys.h"
43#include "comcstfi.h"
44#include "comgeomfi.h"
45#include "tracer.h"
46#include "comvert.h"
47
48      INTEGER ngrid,nlayer,nq
49
50      REAL pplay(ngrid,nlayer)
51      REAL pplev(ngrid,nlayer+1)
52      REAL pq(ngrid,nlayer,nq)
53      REAL aerosol(ngrid,nlayer,naerkind)
54      REAL reffrad(ngrid,nlayer,naerkind)
55      REAL QREFvis3d(ngridmx,nlayermx,naerkind)
56      REAL QREFir3d(ngridmx,nlayermx,naerkind)
57
58      REAL tau_col(ngrid)
59!      REAL tauref(ngrid), tau_col(ngrid)
60
61      real cloudfrac(ngridmx,nlayermx)
62      real aerosol0
63
64      INTEGER l,ig,iq,iaer
65
66      LOGICAL firstcall
67      DATA firstcall/.true./
68      SAVE firstcall
69
70      REAL CBRT
71      EXTERNAL CBRT
72
73      INTEGER,SAVE :: i_co2ice=0      ! co2 ice
74      INTEGER,SAVE :: i_h2oice=0      ! water ice
75      CHARACTER(LEN=20) :: tracername ! to temporarily store text
76
77      ! for fixed dust profiles
78      real topdust, expfactor, zp
79      REAL taudusttmp(ngridmx) ! Temporary dust opacity used before scaling
80      REAL tauh2so4tmp(ngridmx) ! Temporary h2so4 opacity used before scaling
81      ! BENJAMIN MODIFS
82      real CLFtot,CLF1,CLF2
83      real totcloudfrac(ngridmx)
84      logical clearsky
85
86      ! identify tracers
87      IF (firstcall) THEN
88
89         ! are these tests of any real use ?
90        IF(ngrid.NE.ngridmx) THEN
91            PRINT*,'STOP in aeropacity!'
92            PRINT*,'problem of dimensions:'
93            PRINT*,'ngrid  =',ngrid
94            PRINT*,'ngridmx  =',ngridmx
95            STOP
96        ENDIF
97
98        if (nq.gt.nqmx) then
99           write(*,*) 'STOP in aeropacity: (nq .gt. nqmx)!'
100           write(*,*) 'nq=',nq,' nqmx=',nqmx
101           stop
102        endif
103
104        write(*,*) "Tracers found in aeropacity:"
105        do iq=1,nqmx
106          tracername=noms(iq)
107          if (tracername.eq."co2_ice") then
108            i_co2ice=iq
109          write(*,*) "i_co2ice=",i_co2ice
110
111          endif
112          if (tracername.eq."h2o_ice") then
113            i_h2oice=iq
114            write(*,*) "i_h2oice=",i_h2oice
115          endif
116        enddo
117
118        if (naerkind.ne.0) then
119          print*, "If you would like to use aerosols, make sure any old"
120          print*, "start files are updated in newstart using the option"
121          print*, "q=0"
122          write(*,*) "Aerosols found in aeropacity:"
123        endif
124
125        if (iaero_co2.ne.0) then
126          print*, 'iaero_co2=  ',iaero_co2
127        endif
128        if (iaero_h2o.ne.0) then
129          print*,'iaero_h2o=  ',iaero_h2o   
130        endif
131        if (iaero_dust.ne.0) then
132          print*,'iaero_dust= ',iaero_dust
133        endif
134        if (iaero_h2so4.ne.0) then
135          print*,'iaero_h2so4= ',iaero_h2so4
136        endif
137
138        firstcall=.false.
139      ENDIF ! of IF (firstcall)
140
141
142!     ---------------------------------------------------------
143!==================================================================
144!    CO2 ice aerosols
145!==================================================================
146
147        if (iaero_co2.ne.0) then
148           iaer=iaero_co2
149!       1. Initialization
150            aerosol(:,:,iaer)=0.0
151!       2. Opacity calculation
152            if (aerofixco2.or.(i_co2ice.eq.0)) then !  CO2 ice cloud prescribed
153               do ig=1, ngrid
154                  do l=1,nlayer
155                     aerosol(ig,l,iaer)=1.e-9
156                  end do
157                  !aerosol(ig,12,iaer)=4.0 ! single cloud layer option
158               end do
159            else
160               DO ig=1, ngrid
161                  DO l=1,nlayer-1 ! to stop the rad tran bug
162
163                     aerosol0 =                         &
164                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
165                          ( rho_co2 * reffrad(ig,l,iaer) )  ) *   &
166                          ( pq(ig,l,i_co2ice) + 1.E-9 ) *         &
167                          ( pplev(ig,l) - pplev(ig,l+1) ) / g
168                     aerosol0           = max(aerosol0,1.e-9)
169                     aerosol0           = min(aerosol0,L_TAUMAX)
170                     aerosol(ig,l,iaer) = aerosol0
171!                     aerosol(ig,l,iaer) = 0.0
172!                     print*, aerosol(ig,l,iaer)
173!        using cloud fraction
174!                     aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF))
175!                     aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),L_TAUMAX)
176
177
178                  ENDDO
179               ENDDO
180            end if ! if fixed or varying
181        end if ! if CO2 aerosols   
182!==================================================================
183!     Water ice / liquid
184!==================================================================
185
186        if (iaero_h2o.ne.0) then
187           iaer=iaero_h2o
188!       1. Initialization
189            aerosol(:,:,iaer)=0.0
190!       2. Opacity calculation
191            if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then
192               do ig=1,ngrid ! to stop the rad tran bug
193                  do l=1,nlayer
194                     aerosol(ig,l,iaer) =1.e-9
195                  end do
196               end do
197
198               ! put cloud at cloudlvl
199               if(kastprof.and.(cloudlvl.ne.0.0))then
200                  ig=1
201                  do l=1,nlayer
202                     if(int(cloudlvl).eq.l)then
203                     !if(cloudlvl.gt.(pplay(ig,l)/pplev(ig,1)))then
204                        print*,'Inserting cloud at level ',l
205                        !aerosol(ig,l,iaer)=10.0
206
207                        rho_ice=920.0
208
209                        ! the Kasting approximation
210                        aerosol(ig,l,iaer) =                      &
211                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
212                          ( rho_ice * reffrad(ig,l,iaer) )  ) *   &
213                          !( pq(ig,l,i_h2oice) + 1.E-9 ) *         &
214                          ( 4.0e-4 + 1.E-9 ) *         &
215                          ( pplev(ig,l) - pplev(ig,l+1) ) / g
216
217
218                        open(115,file='clouds.out',form='formatted')
219                        write(115,*) l,aerosol(ig,l,iaer)
220                        close(115)
221
222                        return
223                     endif
224                  end do
225
226                  call abort
227               endif
228
229            else
230               do ig=1, ngrid
231                  do l=1,nlayer-1 ! to stop the rad tran bug
232
233                     if(CLFvarying)then
234                        CLF1    = max(cloudfrac(ig,l),0.01)
235                        CLFtot  = max(totcloudfrac(ig),0.01)
236                        CLF2    = CLF1/CLFtot
237                        CLF2    = min(1.,CLF2)
238                     else
239                        CLF1=1.0
240                        CLF2=CLFfixval
241                     endif
242                     
243                     aerosol0 =                                   &
244                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
245                          ( rho_ice * reffrad(ig,l,iaer) )  ) *   &
246                          !  pq(ig,l,i_h2oice) *                   & !JL I dropped the +1e-9 here to have the same
247                          !( pplev(ig,l) - pplev(ig,l+1) ) / g /   & !   opacity in the clearsky=true and the
248                          !  CLF1                                   !   clear=false/pq=0 case
249                          ( pq(ig,l,i_h2oice) + 1.E-9 ) *         & ! Doing this makes the code unstable, so I have restored it (RW)
250                          ( pplev(ig,l) - pplev(ig,l+1) ) / g /   & 
251                            CLF1
252
253                     aerosol(ig,l,iaer) = -log(1 - CLF2 + CLF2*exp(-aerosol0))
254                     ! why no division in exponential?
255                     ! because its already done in aerosol0
256                     
257                     aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
258                     aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),aerosol0)                     
259                     
260                  enddo
261               enddo
262            end if
263        end if ! End if h2o aerosol
264
265!==================================================================
266!             Dust
267!==================================================================
268        if (iaero_dust.ne.0) then
269          iaer=iaero_dust
270
271!         1. Initialization
272            aerosol(:,:,iaer)=0.0
273         
274            topdust=30.0 ! km  (used to be 10.0 km) LK
275
276!       2. Opacity calculation
277
278!           expfactor=0.
279           DO l=1,nlayer-1
280             DO ig=1,ngrid
281!             Typical mixing ratio profile
282
283                 zp=(preff/pplay(ig,l))**(70./topdust)
284                 expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
285
286!             Vertical scaling function
287              aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) &
288               *expfactor
289
290
291             ENDDO
292           ENDDO
293
294!          Rescaling each layer to reproduce the choosen (or assimilated)
295!          dust extinction opacity at visible reference wavelength, which
296!          is scaled to the "preff" reference surface pressure available in comvert.h
297!          and stored in startfi.nc
298
299            taudusttmp(1:ngrid)=0.
300              DO l=1,nlayer
301                DO ig=1,ngrid
302                   taudusttmp(ig) = taudusttmp(ig) &
303                          +  aerosol(ig,l,iaer)
304                ENDDO
305              ENDDO
306            DO l=1,nlayer-1
307               DO ig=1,ngrid
308                  aerosol(ig,l,iaer) = max(1E-20, &
309                          dusttau &
310                       *  pplev(ig,1) / preff &
311                       *  aerosol(ig,l,iaer) &
312                       /  taudusttmp(ig))
313
314              ENDDO
315            ENDDO
316        end if ! If dust aerosol   
317
318!==================================================================
319!           H2SO4
320!==================================================================
321! added by LK
322      if (iaero_h2so4.ne.0) then
323          iaer=iaero_h2so4
324
325!       1. Initialization
326            aerosol(:,:,iaer)=0.0
327
328
329!       2. Opacity calculation
330
331!           expfactor=0.
332           DO l=1,nlayer-1
333             DO ig=1,ngrid
334!             Typical mixing ratio profile
335
336                 zp=(preff/pplay(ig,l))**(70./30) !emulating topdust
337                 expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
338
339!             Vertical scaling function
340              aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) &
341               *expfactor
342
343
344             ENDDO
345           ENDDO
346            tauh2so4tmp(1:ngrid)=0.
347              DO l=1,nlayer
348                DO ig=1,ngrid
349                   tauh2so4tmp(ig) = tauh2so4tmp(ig) &
350                          +  aerosol(ig,l,iaer)
351                ENDDO
352              ENDDO
353            DO l=1,nlayer-1
354               DO ig=1,ngrid
355                  aerosol(ig,l,iaer) = max(1E-20, &
356                          1 &
357                       *  pplev(ig,1) / preff &
358                       *  aerosol(ig,l,iaer) &
359                       /  tauh2so4tmp(ig))
360
361              ENDDO
362            ENDDO
363
364! 1/700. is assuming a "sulfurtau" of 1
365! Sulfur aerosol routine to be improved.
366!                     aerosol0 =                         &
367!                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
368!                          ( rho_h2so4 * reffrad(ig,l,iaer) )  ) *   &
369!                          ( pq(ig,l,i_h2so4) + 1.E-9 ) *         &
370!                          ( pplev(ig,l) - pplev(ig,l+1) ) / g
371!                     aerosol0           = max(aerosol0,1.e-9)
372!                     aerosol0           = min(aerosol0,L_TAUMAX)
373!                     aerosol(ig,l,iaer) = aerosol0
374
375!                  ENDDO
376!               ENDDO
377      end if
378 
379           
380
381! --------------------------------------------------------------------------
382! Column integrated visible optical depth in each point (used for diagnostic)
383
384      tau_col(:)=0.0
385      do iaer = 1, naerkind
386         do l=1,nlayer
387            do ig=1,ngrid
388               tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
389            end do
390         end do
391      end do
392
393      do ig=1, ngrid
394         do l=1,nlayer
395            do iaer = 1, naerkind
396               if(aerosol(ig,l,iaer).gt.1.e3)then
397                  print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
398                  print*,'at ig=',ig,',  l=',l,', iaer=',iaer
399                  print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
400                  print*,'reffrad=',reffrad(ig,l,iaer)
401               endif
402            end do
403         end do
404      end do
405
406      do ig=1, ngrid
407         if(tau_col(ig).gt.1.e3)then
408            print*,'WARNING: tau_col=',tau_col(ig)
409            print*,'at ig=',ig
410            print*,'aerosol=',aerosol(ig,:,:)
411            print*,'QREFvis3d=',QREFvis3d(ig,:,:)
412            print*,'reffrad=',reffrad(ig,:,:)
413         endif
414      end do
415      return
416    end subroutine aeropacity
417     
Note: See TracBrowser for help on using the repository browser.