source: trunk/LMDZ.PLUTO/libf/phypluto/aeropacity.F90

Last change on this file was 3959, checked in by debatzbr, 6 weeks ago

Pluto PCM: Take clouds into account in radiative transfer.
BBT

File size: 8.2 KB
Line 
1module aeropacity_mod
2
3implicit none
4
5contains
6
7      Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,zzlev,pt,pq, &
8         dtau_aer,reffrad,nueffrad, QREFvis3d,QREFir3d,tau_col)
9
10       use radinc_h, only : L_TAUMAX,naerkind
11       use aerosol_mod, only: iaero_haze, i_haze
12       USE tracer_h, only: noms,rho_n2,rho_ice,rho_q,mmol,micro_indx
13       use comcstfi_mod, only: g, pi, mugaz, avocado
14       use geometry_mod, only: latitude
15       use callkeys_mod, only: kastprof, callmufi, callmuclouds
16       use mp2m_diagnostics
17       implicit none
18
19!==================================================================
20!
21!     Purpose
22!     -------
23!     Compute aerosol optical depth in each gridbox.
24!
25!     Authors
26!     -------
27!     F. Forget
28!     F. Montmessin (water ice scheme)
29!     update J.-B. Madeleine (2008)
30!     dust removal, simplification by Robin Wordsworth (2009)
31!
32!     Input
33!     -----
34!     ngrid                            Number of horizontal gridpoints
35!     nlayer                           Number of layers
36!     nq                               Number of tracers
37!     pplev                            Pressure (Pa) at each layer boundary
38!     pq                               Aerosol mixing ratio
39!     reffrad(ngrid,nlayer,naerkind)   Aerosol effective radius
40!     QREFvis3d(ngrid,nlayer,naerkind) \ 3D extinction coefficients
41!     QREFir3d(ngrid,nlayer,naerkind)  / at reference wavelengths
42!
43!     Output
44!     ------
45!     dtau_aer           Aerosol optical depth in layer l, grid point ig
46!     tau_col            Total column optical depth at grid point ig
47!
48!=======================================================================
49
50      INTEGER,INTENT(IN) :: ngrid                            ! Number of atmospheric columns
51      INTEGER,INTENT(IN) :: nlayer                           ! Number of atmospheric layers
52      INTEGER,INTENT(IN) :: nq                               ! Number of tracers
53      REAL,INTENT(IN)    :: pplay(ngrid,nlayer)              ! Mid-layer pressure (Pa)
54      REAL,INTENT(IN)    :: pplev(ngrid,nlayer+1)            ! Inter-layer pressure (Pa)
55      REAL,INTENT(IN)    :: zzlev(ngrid,nlayer)              ! Altitude at the layer boundaries.
56      REAL,INTENT(IN)    :: pq(ngrid,nlayer,nq)              ! Tracers (.../kg_of_air)
57      REAL,INTENT(IN)    :: pt(ngrid,nlayer)                 ! Mid-layer temperature (K)
58      REAL,INTENT(OUT)   :: dtau_aer(ngrid,nlayer,naerkind)  ! Aerosol optical depth
59      REAL,INTENT(IN)    :: reffrad(ngrid,nlayer,naerkind)   ! Aerosol effective radius
60      REAL,INTENT(IN)    :: nueffrad(ngrid,nlayer,naerkind)  ! Aerosol effective variance
61      REAL,INTENT(IN)    :: QREFvis3d(ngrid,nlayer,naerkind) ! Extinction coefficient in the visible
62      REAL,INTENT(IN)    :: QREFir3d(ngrid,nlayer,naerkind)  ! Extinction coefficient in the infrared
63      REAL,INTENT(OUT)   :: tau_col(ngrid)                   ! Column integrated visible optical depth
64
65      real aerosol0, obs_tau_col_aurora, pm
66      real pcloud_deck, cloud_slope
67
68      real dp_strato(ngrid)
69      real dp_tropo(ngrid)
70      real dp_layer(ngrid)
71
72      ! Microphysical tracers
73      real sig
74      real m0as(ngrid,nlayer)
75      real m0af(ngrid,nlayer)
76      real m0ccn(ngrid,nlayer)
77
78      INTEGER l,ig,iq,iaer,ia
79
80      LOGICAL,SAVE :: firstcall=.true.
81!$OMP THREADPRIVATE(firstcall)
82      REAL CBRT
83      EXTERNAL CBRT
84
85      INTEGER,SAVE :: i_n2ice=0      ! n2 ice
86!$OMP THREADPRIVATE(i_n2ice)
87      CHARACTER(LEN=20) :: tracername ! to temporarily store text
88
89      real CLFtot
90      integer igen_ice,igen_gas ! to store the index of generic tracer
91      logical dummy_bool ! dummy boolean just in case we need one
92      !  for venus clouds
93      real      :: p_bot,p_top,h_bot,h_top,mode_dens,h_lay
94
95      ! identify tracers
96      IF (firstcall) THEN
97        ia =0
98        write(*,*) "Tracers found in aeropacity:"
99        do iq=1,nq
100          tracername=noms(iq)
101          if (tracername.eq."n2_ice") then
102            i_n2ice=iq
103          write(*,*) "i_n2ice=",i_n2ice
104
105          endif
106        enddo
107
108        firstcall=.false.
109      ENDIF ! of IF (firstcall)
110
111
112!     ---------------------------------------------------------
113!==================================================================
114!    Haze aerosols
115!==================================================================
116
117      if (iaero_haze.ne.0) then
118         if (callmufi) then
119            ! Convert intensive microphysical tracers to extensive [m-2]
120            m0as(:,:) = pq(:,:,micro_indx(1)) * (pplev(:,1:nlayer) - pplev(:,2:nlayer+1)) / g
121            m0af(:,:) = pq(:,:,micro_indx(3)) * (pplev(:,1:nlayer) - pplev(:,2:nlayer+1)) / g
122           
123            ! Spherical aerosols
124            sig = 0.2
125            where ((m0as(:,:) >= 1e-8).and.(mp2m_rc_sph(:,:) >= 5e-9))
126               dtau_aer(:,:,1) = m0as(:,:) * QREFvis3d(:,:,1) * pi * mp2m_rc_sph(:,:)**2 * exp(2*sig**2)
127            elsewhere
128               dtau_aer(:,:,1) = 0d0
129            endwhere
130            ! Fractal aerosols
131            sig = 0.35
132            where ((m0af(:,:) >= 1e-8).and.(mp2m_rc_fra(:,:) >= 1e-8))
133               dtau_aer(:,:,2) = m0af(:,:) * QREFvis3d(:,:,2) * pi * mp2m_rc_fra(:,:)**2 * exp(2*sig**2)
134            elsewhere
135               dtau_aer(:,:,2) = 0d0
136            endwhere
137           
138            ! Cloud drops
139            if (callmuclouds) then
140               m0ccn(:,:) = pq(:,:,micro_indx(5)) * (pplev(:,1:nlayer) - pplev(:,2:nlayer+1)) / g
141               sig = 0.2
142               where ((m0ccn(:,:) >= 1e-8).and.(mp2m_rc_cld(:,:) >= 5e-9))
143                  dtau_aer(:,:,3) = m0ccn(:,:) * QREFvis3d(:,:,3) * pi * mp2m_rc_cld(:,:)**2 * exp(2*sig**2)
144               elsewhere
145                  dtau_aer(:,:,3) = 0d0
146               endwhere
147            endif ! end callmuclouds
148         
149         else
150            do iaer = 1, naerkind
151               ! 1. Initialization
152               dtau_aer(1:ngrid,1:nlayer,iaer)=0.0
153               ! 2. Opacity calculation
154               DO ig = 1, ngrid
155                  DO l = 1, nlayer-1 ! to stop the rad tran bug
156                     ! If fractal, radius doit etre equivalent sphere radius
157                     ! Eq. 2.37 - Madeleine's PhD (2011).
158                     aerosol0 =                         &
159                        (  0.75 * QREFvis3d(ig,l,iaer) /        &
160                        ( rho_q(i_haze) * reffrad(ig,l,iaer) )  ) *   &
161                        ( pq(ig,l,i_haze) + 1.E-10 ) *         &
162                        ( pplev(ig,l) - pplev(ig,l+1) ) / g
163                     aerosol0            = max(aerosol0,1.e-10)
164                     aerosol0            = min(aerosol0,L_TAUMAX)
165                     dtau_aer(ig,l,iaer) = aerosol0
166                  ENDDO
167               ENDDO
168               !QREF est le meme dans toute la colonne par def si size uniforme
169               !print*, 'TB17: QREFvis3d=',QREFvis3d(1,:,1)
170               !print*, 'TB17: rho_q=',rho_q(i_haze)
171               !print*, 'TB17: reffrad=',reffrad(1,:,1)
172               !print*, 'TB17: pq=',pq(1,:,i_haze)
173               !print*, 'TB17: deltap=',pplev(1,1) - pplev(1,nlayer)
174            enddo ! end iaer = 1, naerkind
175         endif ! end callmufi
176      endif ! if haze aerosols
177
178! --------------------------------------------------------------------------
179! Column integrated visible optical depth in each point (used for diagnostic)
180
181   tau_col(:)=0.0
182   do iaer = 1, naerkind
183      do l=1,nlayer
184         do ig=1,ngrid
185            tau_col(ig) = tau_col(ig) + dtau_aer(ig,l,iaer)
186         end do
187      end do
188   end do
189
190   do ig=1,ngrid
191      do l=1,nlayer
192         do iaer = 1, naerkind
193            if(dtau_aer(ig,l,iaer).gt.1.e3)then
194               print*,'WARNING: dtau_aer=',dtau_aer(ig,l,iaer)
195               print*,'at ig=',ig,',  l=',l,', iaer=',iaer
196               print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
197               print*,'reffrad=',reffrad(ig,l,iaer)
198            endif
199         end do
200      end do
201   end do
202
203   do ig=1,ngrid
204      if(tau_col(ig).gt.1.e3)then
205         print*,'WARNING: tau_col=',tau_col(ig)
206         print*,'at ig=',ig
207         print*,'dtau_aer=',dtau_aer(ig,:,:)
208         print*,'QREFvis3d=',QREFvis3d(ig,:,:)
209         print*,'reffrad=',reffrad(ig,:,:)
210      endif
211   end do
212  end subroutine aeropacity
213
214end module aeropacity_mod
Note: See TracBrowser for help on using the repository browser.