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

Last change on this file since 3813 was 3802, checked in by debatzbr, 7 months ago

Pluto PCM: add some sanity checks to ensure haze stability
BBT

File size: 7.5 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
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)
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
77      INTEGER l,ig,iq,iaer,ia
78
79      LOGICAL,SAVE :: firstcall=.true.
80!$OMP THREADPRIVATE(firstcall)
81      REAL CBRT
82      EXTERNAL CBRT
83
84      INTEGER,SAVE :: i_n2ice=0      ! n2 ice
85!$OMP THREADPRIVATE(i_n2ice)
86      CHARACTER(LEN=20) :: tracername ! to temporarily store text
87
88      real CLFtot
89      integer igen_ice,igen_gas ! to store the index of generic tracer
90      logical dummy_bool ! dummy boolean just in case we need one
91      !  for venus clouds
92      real      :: p_bot,p_top,h_bot,h_top,mode_dens,h_lay
93
94      ! identify tracers
95      IF (firstcall) THEN
96        ia =0
97        write(*,*) "Tracers found in aeropacity:"
98        do iq=1,nq
99          tracername=noms(iq)
100          if (tracername.eq."n2_ice") then
101            i_n2ice=iq
102          write(*,*) "i_n2ice=",i_n2ice
103
104          endif
105        enddo
106
107        firstcall=.false.
108      ENDIF ! of IF (firstcall)
109
110
111!     ---------------------------------------------------------
112!==================================================================
113!    Haze aerosols
114!==================================================================
115
116      if (iaero_haze.ne.0) then
117         if (callmufi) then
118            ! Convert intensive microphysical tracers to extensive [m-2]
119            m0as(:,:) = pq(:,:,micro_indx(1)) * (pplev(:,1:nlayer) - pplev(:,2:nlayer+1)) / g
120            m0af(:,:) = pq(:,:,micro_indx(3)) * (pplev(:,1:nlayer) - pplev(:,2:nlayer+1)) / g
121           
122            ! Spherical aerosols
123            sig = 0.2
124            where ((m0as(:,:) >= 1e-8).and.(mp2m_rc_sph(:,:) >= 5e-9))
125               dtau_aer(:,:,1) = m0as(:,:) * QREFvis3d(:,:,1) * pi * mp2m_rc_sph(:,:)**2 * exp(2*sig**2)
126            elsewhere
127               dtau_aer(:,:,1) = 0d0
128            endwhere
129            ! Fractal aerosols
130            sig = 0.35
131            where ((m0af(:,:) >= 1e-8).and.(mp2m_rc_fra(:,:) >= 1e-8))
132               dtau_aer(:,:,2) = m0af(:,:) * QREFvis3d(:,:,2) * pi * mp2m_rc_fra(:,:)**2 * exp(2*sig**2)
133            elsewhere
134               dtau_aer(:,:,2) = 0d0
135            endwhere
136           
137            ! write(*,*) 'dtau_as  :', MINVAL(dtau_aer(:,:,1)), '-', MAXVAL(dtau_aer(:,:,1))
138            ! write(*,*) 'dtau_af  :', MINVAL(dtau_aer(:,:,2)), '-', MAXVAL(dtau_aer(:,:,2))
139         
140         else
141            do iaer = 1, naerkind
142               ! 1. Initialization
143               dtau_aer(1:ngrid,1:nlayer,iaer)=0.0
144               ! 2. Opacity calculation
145               DO ig = 1, ngrid
146                  DO l = 1, nlayer-1 ! to stop the rad tran bug
147                     ! If fractal, radius doit etre equivalent sphere radius
148                     ! Eq. 2.37 - Madeleine's PhD (2011).
149                     aerosol0 =                         &
150                        (  0.75 * QREFvis3d(ig,l,iaer) /        &
151                        ( rho_q(i_haze) * reffrad(ig,l,iaer) )  ) *   &
152                        ( pq(ig,l,i_haze) + 1.E-10 ) *         &
153                        ( pplev(ig,l) - pplev(ig,l+1) ) / g
154                     aerosol0            = max(aerosol0,1.e-10)
155                     aerosol0            = min(aerosol0,L_TAUMAX)
156                     dtau_aer(ig,l,iaer) = aerosol0
157                  ENDDO
158               ENDDO
159               !QREF est le meme dans toute la colonne par def si size uniforme
160               !print*, 'TB17: QREFvis3d=',QREFvis3d(1,:,1)
161               !print*, 'TB17: rho_q=',rho_q(i_haze)
162               !print*, 'TB17: reffrad=',reffrad(1,:,1)
163               !print*, 'TB17: pq=',pq(1,:,i_haze)
164               !print*, 'TB17: deltap=',pplev(1,1) - pplev(1,nlayer)
165            enddo ! end iaer = 1, naerkind
166         endif ! end callmufi
167      endif ! if haze aerosols
168
169! --------------------------------------------------------------------------
170! Column integrated visible optical depth in each point (used for diagnostic)
171
172   tau_col(:)=0.0
173   do iaer = 1, naerkind
174      do l=1,nlayer
175         do ig=1,ngrid
176            tau_col(ig) = tau_col(ig) + dtau_aer(ig,l,iaer)
177         end do
178      end do
179   end do
180
181   do ig=1,ngrid
182      do l=1,nlayer
183         do iaer = 1, naerkind
184            if(dtau_aer(ig,l,iaer).gt.1.e3)then
185               print*,'WARNING: dtau_aer=',dtau_aer(ig,l,iaer)
186               print*,'at ig=',ig,',  l=',l,', iaer=',iaer
187               print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
188               print*,'reffrad=',reffrad(ig,l,iaer)
189            endif
190         end do
191      end do
192   end do
193
194   do ig=1,ngrid
195      if(tau_col(ig).gt.1.e3)then
196         print*,'WARNING: tau_col=',tau_col(ig)
197         print*,'at ig=',ig
198         print*,'dtau_aer=',dtau_aer(ig,:,:)
199         print*,'QREFvis3d=',QREFvis3d(ig,:,:)
200         print*,'reffrad=',reffrad(ig,:,:)
201      endif
202   end do
203  end subroutine aeropacity
204
205end module aeropacity_mod
Note: See TracBrowser for help on using the repository browser.