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

Last change on this file since 837 was 804, checked in by lkerber, 12 years ago

A bug was found in gfluxi.F. If co2_ice was the only aerosol, the single scattering albedo was occasionally equal to 1. This led to NaNs? in gfluxi.F because we divide by a (1-W0) statement. Now there is a test in gfluxi.F which puts any W0=1 to W0=99999. This change removed the necessity in aeropacity.F90 for aerosols which are supposed to be zero to be put at 1.E-9 (they are now set to zero). In the case of zero aerosols, a dummy co2 aerosol is created and set to zero abundance everywhere. An obsolete test was removed from inifis.F. Some extraneous print statements were removed from physiq.F90.

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