source: trunk/LMDZ.PLUTO.old/libf/phypluto/aeropacity.F90 @ 3436

Last change on this file since 3436 was 3175, checked in by emillour, 11 months ago

Pluto PCM:
Add the old Pluto LMDZ for reference (required prior step to making
an LMDZ.PLUTO using the same framework as the other physics packages).
TB+EM

File size: 5.9 KB
RevLine 
[3175]1      Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pq, &
2         aerosol,reffrad,QREFvis3d,QREFir3d,tau_col)
3
4       use radinc_h, only : L_TAUMAX,naerkind
5       use aerosol_mod
6       use radii_mod
7                 
8       implicit none
9
10#include "dimensions.h"
11#include "dimphys.h"
12#include "comcstfi.h"
13#include "callkeys.h"
14#include "tracer.h"
15
16!==================================================================
17!     
18!     Purpose
19!     -------
20!     Compute aerosol optical depth in each gridbox.
21!     
22!     Authors
23!     -------
24!     F. Forget
25!     F. Montmessin (water ice scheme)
26!     update J.-B. Madeleine (2008)
27!     Adaptation Pluto : Tanguy Bertrand (2017)
28!
29!     Input
30!     -----
31!     ngrid             Number of horizontal gridpoints
32!     nlayer            Number of layers
33!     nq                Number of tracers
34!     pplev             Pressure (Pa) at each layer boundary
35!     pq                Aerosol mixing ratio
36!     reffrad(ngrid,nlayer,naerkind)         Aerosol effective radius
37!     QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients
38!     QREFir3d(ngrid,nlayer,naerkind)  / at reference wavelengths
39!
40!     Output
41!     ------
42!     aerosol            Aerosol optical depth in layer l, grid point ig
43!     tau_col            Total column optical depth at grid point ig
44!
45!=======================================================================
46
47      INTEGER,INTENT(IN) :: ngrid  ! number of atmospheric columns
48      INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
49      INTEGER,INTENT(IN) :: nq     ! number of tracers
50      REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
51      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
52      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
53      REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol optical depth
54      REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind) ! aerosol effective radius
55      REAL,INTENT(IN) :: QREFvis3d(ngrid,nlayer,naerkind) ! extinction coefficient in the visible
56      REAL,INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind)
57      REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth
58
59      real aerosol0
60
61      INTEGER l,ig,iq,iaer
62
63      LOGICAL,SAVE :: firstcall=.true.
64!$OMP THREADPRIVATE(firstcall)
65      REAL CBRT
66      EXTERNAL CBRT
67
68      !INTEGER,SAVE :: i_haze=0      !
69      CHARACTER(LEN=20) :: tracername ! to temporarily store text
70
71      ! for fixed dust profiles
72      real topdust, expfactor, zp
73      REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling
74      REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling
75
76      real CLFtot
77!      integer block
78
79      ! identify tracers
80      IF (firstcall) THEN
81
82        write(*,*) "Tracers found in aeropacity:"
83        do iq=1,nq
84          tracername=noms(iq)
85          write(*,*) iq,tracername
86        enddo
87
88        if (iaero_haze.eq.0) then
89          print*, "No active aerosols found in aeropacity"
90          print*, "Check traceur.def or haze option in callphys.def"
91          stop
92        else
93          write(*,*) "Active aerosols found in aeropacity:"
94        endif
95
96        if (iaero_haze.ne.0) then
97          print*, 'iaero_haze=  ',iaero_haze
98        endif
99        print*, 'naerkind aeropacity=',naerkind
100        firstcall=.false.
101      ENDIF ! of IF (firstcall)
102
103
104!     ---------------------------------------------------------
105!==================================================================
106!    Haze aerosols
107!==================================================================
108
109      if (iaero_haze.ne.0) then
110           iaer=iaero_haze
111!       1. Initialization
112            aerosol(1:ngrid,1:nlayer,iaer)=0.0
113!       2. Opacity calculation
114            DO ig=1, ngrid
115                  DO l=1,nlayer-1 ! to stop the rad tran bug
116                     ! if fractal, radius doit etre equivalent sphere radius
117                     aerosol0 =                         &
118                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
119                          ( rho_q(i_haze) * reffrad(ig,l,iaer) )  ) *   &
120                          ( pq(ig,l,i_haze) + 1.E-10 ) *         &
121                          ( pplev(ig,l) - pplev(ig,l+1) ) / g
122                     aerosol0           = max(aerosol0,1.e-10)
123                     aerosol0           = min(aerosol0,L_TAUMAX)
124                     aerosol(ig,l,iaer) = aerosol0
125                  ENDDO
126            ENDDO
127            !QREF est le meme dans toute la colonne par def si size uniforme
128            !print*, 'TB17: QREFvis3d=',QREFvis3d(1,:,1)           
129            !print*, 'TB17: rho_q=',rho_q(i_haze)           
130            !print*, 'TB17: reffrad=',reffrad(1,:,1)       
131            !print*, 'TB17: pq=',pq(1,:,i_haze)
132            !print*, 'TB17: deltap=',pplev(1,1) - pplev(1,nlayer)
133      end if ! if haze aerosols   
134           
135! --------------------------------------------------------------------------
136! Column integrated visible optical depth in each point (used for diagnostic)
137
138      tau_col(:)=0.0
139      do iaer = 1, naerkind
140         do l=1,nlayer
141            do ig=1,ngrid
142               tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
143            end do
144         end do
145      end do
146
147      do ig=1,ngrid
148         do l=1,nlayer
149            do iaer = 1, naerkind
150               if(aerosol(ig,l,iaer).gt.1.e3)then
151                  print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
152                  print*,'at ig=',ig,',  l=',l,', iaer=',iaer
153                  print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
154                  print*,'reffrad=',reffrad(ig,l,iaer)
155               endif
156            end do
157         end do
158      end do
159
160      do ig=1,ngrid
161         if(tau_col(ig).gt.1.e3)then
162            print*,'WARNING: tau_col=',tau_col(ig)
163            print*,'at ig=',ig
164            print*,'aerosol=',aerosol(ig,:,:)
165            print*,'QREFvis3d=',QREFvis3d(ig,:,:)
166            print*,'reffrad=',reffrad(ig,:,:)
167         endif
168      end do
169      return
170    end subroutine aeropacity
171     
Note: See TracBrowser for help on using the repository browser.