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

Last change on this file since 1016 was 858, checked in by emillour, 13 years ago

Generic GCM:

  • Fixed sedimentation issue: ensure in callsedim that the correct radii are provided to newsedim and also that the updated temperature and tracer mixing ratios are used to compute sedimentation.
  • Updated the way aerosol radii are considered and used; routines in radii_mod (h2o_reffrad, co2_reffrad, etc.) only handle a single aerosol. The idea here is that these can be called from anywhere and that the caller doesn't need to have the full (naerkind size) array of aerosol radii.
  • cleanup (addition of intent(..) to routine arguments) in various routines

EM

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