source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/aeropacity.F90 @ 263

Last change on this file since 263 was 227, checked in by milmd, 10 years ago

Last LMDZ version (1315) with OpenMP directives and other stuff

File size: 16.1 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       USE comgeomfi_h
7       USE tracer_h, only: noms,rho_co2,rho_ice
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
32!     QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients
33!     QREFir3d(ngrid,nlayer,naerkind)  / at reference wavelengths
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"
46!#include "comvert.h"
47
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,nlayer,naerkind) ! extinction coefficient in the visible
57      REAL,INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind)
58      REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth
59      ! BENJAMIN MODIFS
60      real,intent(in) :: cloudfrac(ngrid,nlayer) ! cloud fraction
61      real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
62      logical,intent(in) :: clearsky
63
64      real aerosol0
65
66      INTEGER l,ig,iq,iaer
67
68      LOGICAL,SAVE :: firstcall=.true.
69!$OMP THREADPRIVATE(firstcall)
70      REAL CBRT
71      EXTERNAL CBRT
72
73      INTEGER,SAVE :: i_co2ice=0      ! co2 ice
74      INTEGER,SAVE :: i_h2oice=0      ! water ice
75!$OMP THREADPRIVATE(i_co2ice,i_h2oice)
76      CHARACTER(LEN=20) :: tracername ! to temporarily store text
77
78      ! for fixed dust profiles
79      real topdust, expfactor, zp
80      REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling
81      REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling
82
83      real CLFtot
84
85      ! identify tracers
86      IF (firstcall) THEN
87
88        write(*,*) "Tracers found in aeropacity:"
89        do iq=1,nq
90          tracername=noms(iq)
91          if (tracername.eq."co2_ice") then
92            i_co2ice=iq
93          write(*,*) "i_co2ice=",i_co2ice
94
95          endif
96          if (tracername.eq."h2o_ice") then
97            i_h2oice=iq
98            write(*,*) "i_h2oice=",i_h2oice
99          endif
100        enddo
101
102        if (noaero) then
103          print*, "No active aerosols found in aeropacity"
104        else
105          print*, "If you would like to use aerosols, make sure any old"
106          print*, "start files are updated in newstart using the option"
107          print*, "q=0"
108          write(*,*) "Active aerosols found in aeropacity:"
109        endif
110
111        if ((iaero_co2.ne.0).and.(.not.noaero)) then
112          print*, 'iaero_co2=  ',iaero_co2
113        endif
114        if (iaero_h2o.ne.0) then
115          print*,'iaero_h2o=  ',iaero_h2o   
116        endif
117        if (iaero_dust.ne.0) then
118          print*,'iaero_dust= ',iaero_dust
119        endif
120        if (iaero_h2so4.ne.0) then
121          print*,'iaero_h2so4= ',iaero_h2so4
122        endif
123        if (iaero_back2lay.ne.0) then
124          print*,'iaero_back2lay= ',iaero_back2lay
125        endif
126
127        firstcall=.false.
128      ENDIF ! of IF (firstcall)
129
130
131!     ---------------------------------------------------------
132!==================================================================
133!    CO2 ice aerosols
134!==================================================================
135
136      if (iaero_co2.ne.0) then
137           iaer=iaero_co2
138!       1. Initialization
139            aerosol(1:ngrid,1:nlayer,iaer)=0.0
140!       2. Opacity calculation
141            if (noaero) then ! aerosol set to zero
142             aerosol(1:ngrid,1:nlayer,iaer)=0.0
143            elseif (aerofixco2.or.(i_co2ice.eq.0)) then !  CO2 ice cloud prescribed
144               aerosol(1:ngrid,1:nlayer,iaer)=1.e-9
145               !aerosol(1:ngrid,12,iaer)=4.0 ! single cloud layer option
146            else
147               DO ig=1, ngrid
148                  DO l=1,nlayer-1 ! to stop the rad tran bug
149
150                     aerosol0 =                         &
151                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
152                          ( rho_co2 * reffrad(ig,l,iaer) )  ) *   &
153                          ( pq(ig,l,i_co2ice) + 1.E-9 ) *         &
154                          ( pplev(ig,l) - pplev(ig,l+1) ) / g
155                     aerosol0           = max(aerosol0,1.e-9)
156                     aerosol0           = min(aerosol0,L_TAUMAX)
157                     aerosol(ig,l,iaer) = aerosol0
158!                     aerosol(ig,l,iaer) = 0.0
159!                     print*, aerosol(ig,l,iaer)
160!        using cloud fraction
161!                     aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF))
162!                     aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),L_TAUMAX)
163
164
165                  ENDDO
166               ENDDO
167            end if ! if fixed or varying
168      end if ! if CO2 aerosols   
169!==================================================================
170!     Water ice / liquid
171!==================================================================
172
173      if (iaero_h2o.ne.0) then
174           iaer=iaero_h2o
175!       1. Initialization
176            aerosol(1:ngrid,1:nlayer,iaer)=0.0
177!       2. Opacity calculation
178            if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then
179               aerosol(1:ngrid,1:nlayer,iaer) =1.e-9
180
181               ! put cloud at cloudlvl
182               if(kastprof.and.(cloudlvl.ne.0.0))then
183                  ig=1
184                  do l=1,nlayer
185                     if(int(cloudlvl).eq.l)then
186                     !if(cloudlvl.gt.(pplay(ig,l)/pplev(ig,1)))then
187                        print*,'Inserting cloud at level ',l
188                        !aerosol(ig,l,iaer)=10.0
189
190                        rho_ice=920.0
191
192                        ! the Kasting approximation
193                        aerosol(ig,l,iaer) =                      &
194                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
195                          ( rho_ice * reffrad(ig,l,iaer) )  ) *   &
196                          !( pq(ig,l,i_h2oice) + 1.E-9 ) *         &
197                          ( 4.0e-4 + 1.E-9 ) *         &
198                          ( pplev(ig,l) - pplev(ig,l+1) ) / g
199
200
201                        open(115,file='clouds.out',form='formatted')
202                        write(115,*) l,aerosol(ig,l,iaer)
203                        close(115)
204
205                        return
206                     endif
207                  end do
208
209                  call abort
210               endif
211
212            else
213
214               do ig=1, ngrid
215                  do l=1,nlayer-1 ! to stop the rad tran bug
216
217                     aerosol(ig,l,iaer) =                                    & !modification by BC
218                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
219                          ( rho_ice * reffrad(ig,l,iaer) )  ) *   &
220                          !  pq(ig,l,i_h2oice) *                   & !JL I dropped the +1e-9 here to have the same
221                          !( pplev(ig,l) - pplev(ig,l+1) ) / g       !   opacity in the clearsky=true and the
222                                                                     !   clear=false/pq=0 case
223                          ( pq(ig,l,i_h2oice) + 1.E-9 ) *         & ! Doing this makes the code unstable, so I have restored it (RW)
224                          ( pplev(ig,l) - pplev(ig,l+1) ) / g 
225
226                  enddo
227               enddo
228
229               if(CLFvarying)then
230                  call totalcloudfrac(ngrid,nlayer,nq,cloudfrac,totcloudfrac,pplev,pq,aerosol(1,1,iaer))
231                  do ig=1, ngrid
232                     do l=1,nlayer-1 ! to stop the rad tran bug
233                        CLFtot  = max(totcloudfrac(ig),0.01)
234                        aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
235                        aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
236                     enddo
237                  enddo
238               else
239                  do ig=1, ngrid
240                     do l=1,nlayer-1 ! to stop the rad tran bug
241                        CLFtot  = CLFfixval
242                        aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
243                        aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
244                     enddo
245                  enddo
246              end if!(CLFvarying)
247            endif !(aerofixed.or.(i_h2oice.eq.0).or.clearsky)
248             
249      end if ! End if h2o aerosol
250
251!==================================================================
252!             Dust
253!==================================================================
254      if (iaero_dust.ne.0) then
255          iaer=iaero_dust
256!         1. Initialization
257          aerosol(1:ngrid,1:nlayer,iaer)=0.0
258         
259          topdust=30.0 ! km  (used to be 10.0 km) LK
260
261!       2. Opacity calculation
262
263!           expfactor=0.
264           DO l=1,nlayer-1
265             DO ig=1,ngrid
266!             Typical mixing ratio profile
267
268                 zp=(pplev(ig,1)/pplay(ig,l))**(70./topdust)
269                 expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
270
271!             Vertical scaling function
272              aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) &
273               *expfactor
274
275
276             ENDDO
277           ENDDO
278
279!          Rescaling each layer to reproduce the choosen (or assimilated)
280!          dust extinction opacity at visible reference wavelength, which
281!          is scaled to the surface pressure pplev(ig,1)
282
283            taudusttmp(1:ngrid)=0.
284              DO l=1,nlayer
285                DO ig=1,ngrid
286                   taudusttmp(ig) = taudusttmp(ig) &
287                          +  aerosol(ig,l,iaer)
288                ENDDO
289              ENDDO
290            DO l=1,nlayer-1
291               DO ig=1,ngrid
292                  aerosol(ig,l,iaer) = max(1E-20, &
293                          dusttau &
294                       *  pplev(ig,1) / pplev(ig,1) &
295                       *  aerosol(ig,l,iaer) &
296                       /  taudusttmp(ig))
297
298              ENDDO
299            ENDDO
300      end if ! If dust aerosol   
301
302!==================================================================
303!           H2SO4
304!==================================================================
305! added by LK
306      if (iaero_h2so4.ne.0) then
307         iaer=iaero_h2so4
308
309!       1. Initialization
310         aerosol(1:ngrid,1:nlayer,iaer)=0.0
311
312
313!       2. Opacity calculation
314
315!           expfactor=0.
316         DO l=1,nlayer-1
317            DO ig=1,ngrid
318!              Typical mixing ratio profile
319
320               zp=(pplev(ig,1)/pplay(ig,l))**(70./30) !emulating topdust
321               expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
322
323!             Vertical scaling function
324               aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1))*expfactor
325
326            ENDDO
327         ENDDO
328         tauh2so4tmp(1:ngrid)=0.
329         DO l=1,nlayer
330            DO ig=1,ngrid
331               tauh2so4tmp(ig) = tauh2so4tmp(ig) + aerosol(ig,l,iaer)
332            ENDDO
333         ENDDO
334         DO l=1,nlayer-1
335            DO ig=1,ngrid
336               aerosol(ig,l,iaer) = max(1E-20, &
337                          1 &
338                       *  pplev(ig,1) / pplev(ig,1) &
339                       *  aerosol(ig,l,iaer) &
340                       /  tauh2so4tmp(ig))
341
342            ENDDO
343         ENDDO
344
345! 1/700. is assuming a "sulfurtau" of 1
346! Sulfur aerosol routine to be improved.
347!                     aerosol0 =                         &
348!                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
349!                          ( rho_h2so4 * reffrad(ig,l,iaer) )  ) *   &
350!                          ( pq(ig,l,i_h2so4) + 1.E-9 ) *         &
351!                          ( pplev(ig,l) - pplev(ig,l+1) ) / g
352!                     aerosol0           = max(aerosol0,1.e-9)
353!                     aerosol0           = min(aerosol0,L_TAUMAX)
354!                     aerosol(ig,l,iaer) = aerosol0
355
356!                  ENDDO
357!               ENDDO
358      end if
359 
360           
361!     ---------------------------------------------------------
362!==================================================================
363!    Two-layer aerosols (unknown composition)
364!    S. Guerlet (2013)
365!==================================================================
366
367      if (iaero_back2lay .ne.0) then
368           iaer=iaero_back2lay
369!       1. Initialization
370            aerosol(1:ngrid,1:nlayer,iaer)=0.0
371!       2. Opacity calculation
372          DO ig=1,ngrid
373           DO l=1,nlayer-1
374             aerosol(ig,l,iaer) = ( pplev(ig,l) - pplev(ig,l+1) )
375             !! 1. below tropospheric layer: no aerosols
376             IF (pplev(ig,l) .gt. pres_bottom_tropo) THEN
377               aerosol(ig,l,iaer) = 0.*aerosol(ig,l,iaer)
378             !! 2. tropo layer
379             ELSEIF (pplev(ig,l) .le. pres_bottom_tropo .and. pplev(ig,l) .ge. pres_top_tropo) THEN
380               aerosol(ig,l,iaer) = obs_tau_col_tropo*aerosol(ig,l,iaer)
381             !! 3. linear transition
382             ELSEIF (pplev(ig,l) .lt. pres_top_tropo .and. pplev(ig,l) .gt. pres_bottom_strato) THEN
383               expfactor=log(obs_tau_col_strato/obs_tau_col_tropo)/log(pres_bottom_strato/pres_top_tropo)
384               aerosol(ig,l,iaer)= obs_tau_col_tropo*((pplev(ig,l)/pres_top_tropo)**expfactor)*aerosol(ig,l,iaer)/1.5
385             !! 4. strato layer
386             ELSEIF (pplev(ig,l) .le. pres_bottom_strato .and. pplev(ig,l) .gt. pres_top_strato) THEN
387               aerosol(ig,l,iaer)= obs_tau_col_strato*aerosol(ig,l,iaer)
388             !! 5. above strato layer: no aerosols
389             ELSEIF (pplev(ig,l) .lt. pres_top_strato) THEN
390               aerosol(ig,l,iaer) = 0.*aerosol(ig,l,iaer)
391             ENDIF
392           ENDDO
393          ENDDO
394
395 !       3. Re-normalize to observed total column
396         tau_col(:)=0.0
397         DO l=1,nlayer
398          DO ig=1,ngrid
399               tau_col(ig) = tau_col(ig) &
400                     + aerosol(ig,l,iaer)/(obs_tau_col_tropo+obs_tau_col_strato)
401            ENDDO
402         ENDDO
403
404         DO ig=1,ngrid
405           DO l=1,nlayer-1
406                aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/tau_col(ig)
407           ENDDO
408         ENDDO
409
410
411      end if ! if Two-layer aerosols 
412
413
414! --------------------------------------------------------------------------
415! Column integrated visible optical depth in each point (used for diagnostic)
416
417      tau_col(:)=0.0
418      do iaer = 1, naerkind
419         do l=1,nlayer
420            do ig=1,ngrid
421               tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
422            end do
423         end do
424      end do
425
426      do ig=1,ngrid
427         do l=1,nlayer
428            do iaer = 1, naerkind
429               if(aerosol(ig,l,iaer).gt.1.e3)then
430                  print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
431                  print*,'at ig=',ig,',  l=',l,', iaer=',iaer
432                  print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
433                  print*,'reffrad=',reffrad(ig,l,iaer)
434               endif
435            end do
436         end do
437      end do
438
439      do ig=1,ngrid
440         if(tau_col(ig).gt.1.e3)then
441            print*,'WARNING: tau_col=',tau_col(ig)
442            print*,'at ig=',ig
443            print*,'aerosol=',aerosol(ig,:,:)
444            print*,'QREFvis3d=',QREFvis3d(ig,:,:)
445            print*,'reffrad=',reffrad(ig,:,:)
446         endif
447      end do
448      return
449    end subroutine aeropacity
450     
Note: See TracBrowser for help on using the repository browser.