source: trunk/LMDZ.VENUS/libf/phyvenus/aeropacity.F90 @ 3094

Last change on this file since 3094 was 2622, checked in by slebonnois, 3 years ago

SL: VENUS update (i) bug correction (2 bugs, phytrac and physiq), affected meam molec mass computations... (ii) updates for VCD 2.0 (iii) aeropacity: for latitudinal variations of the cloud distribution

File size: 13.8 KB
Line 
1      Subroutine aeropacity(ngrid,nlayer,pplay,pplev,pt, &
2         aerosol,reffrad,nueffrad,QREFvis3d,tau_col)
3
4       use radinc_h, only : naerkind
5       use aerosol_mod
6                   
7       implicit none
8#include "YOMCST.h"
9
10!==================================================================
11!     
12!     Purpose
13!     -------
14!     Compute aerosol optical depth in each gridbox.
15!     
16!     Authors
17!     -------
18!     F. Forget
19!     F. Montmessin (water ice scheme)
20!     update J.-B. Madeleine (2008)
21!     dust removal, simplification by Robin Wordsworth (2009)
22!     Generic n-layer aerosol - J. Vatant d'Ollone (2020)
23!
24!     Input
25!     -----
26!     ngrid             Number of horizontal gridpoints
27!     nlayer            Number of layers
28!     pplev             Pressure (Pa) at each layer boundary
29!     reffrad(ngrid,nlayer,naerkind)         Aerosol effective radius
30!     QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients
31!     QREFir3d(ngrid,nlayer,naerkind)  / at reference wavelengths
32!
33!     Output
34!     ------
35!     aerosol            Aerosol optical depth in layer l, grid point ig
36!     tau_col            Total column optical depth at grid point ig
37!
38!=======================================================================
39
40      INTEGER,INTENT(IN) :: ngrid  ! number of atmospheric columns
41      INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
42      REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
43      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
44      REAL,INTENT(IN) :: pt(ngrid,nlayer) ! mid-layer temperatre (K)
45      REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol optical depth
46      REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind) ! aerosol effective radius
47      REAL,INTENT(IN) :: nueffrad(ngrid,nlayer,naerkind) ! aerosol effective variance
48      REAL,INTENT(IN) :: QREFvis3d(ngrid,nlayer,naerkind) ! extinction coefficient in the visible
49      REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth
50
51      INTEGER l,ig,iaer,ll
52
53      LOGICAL,SAVE :: firstcall=.true.
54!$OMP THREADPRIVATE(firstcall)
55      character*80 abort_message
56      character*20 modname
57
58      !  for venus clouds
59      real,ALLOCATABLE,SAVE :: pbot2(:,:),pbot(:,:),ptop(:,:),ptop2(:,:)
60      real,ALLOCATABLE,SAVE :: hbot2(:,:),hbot(:,:),htop(:,:),htop2(:,:)
61      real,ALLOCATABLE,SAVE :: nmode(:,:)
62!$OMP THREADPRIVATE(pbot2,pbot,ptop,ptop2)
63!$OMP THREADPRIVATE(hbot2,hbot,htop,htop2)
64!$OMP THREADPRIVATE(nmode)
65      real :: mode_dens,hlay
66
67! ---------------------------------------------------------
68! ---------------------------------------------------------
69! First call only
70      IF (firstcall) THEN
71     
72        modname = 'aeropacity'
73
74! verify tracers
75        write(*,*) "Active aerosols found in aeropacity, designed for Venus"
76
77        if (iaero_venus1.eq.1) then
78          print*,'iaero_venus1= ',iaero_venus1
79        else
80            abort_message='iaero_venus1 is not 1...'
81            call abort_physic(modname,abort_message,1)
82        endif
83        if (iaero_venus2.eq.2) then
84          print*,'iaero_venus2= ',iaero_venus2
85        else
86            abort_message='iaero_venus2 is not 2...'
87            call abort_physic(modname,abort_message,1)
88        endif
89        if (iaero_venus2p.eq.3) then
90          print*,'iaero_venus2p= ',iaero_venus2p
91        else
92            abort_message='iaero_venus2p is not 3...'
93            call abort_physic(modname,abort_message,1)
94        endif
95        if (iaero_venus3.eq.4) then
96          print*,'iaero_venus3= ',iaero_venus3
97        else
98            abort_message='iaero_venus3 is not 4...'
99            call abort_physic(modname,abort_message,1)
100        endif
101        if (iaero_venusUV.eq.5) then
102          print*,'iaero_venusUV= ',iaero_venusUV
103        else
104            abort_message='iaero_venus1UV is not 5...'
105            call abort_physic(modname,abort_message,1)
106        endif
107
108! cloud model
109        allocate(pbot2(ngrid,naerkind),pbot(ngrid,naerkind))
110        allocate(ptop(ngrid,naerkind),ptop2(ngrid,naerkind))
111        allocate(hbot2(ngrid,naerkind),hbot(ngrid,naerkind))
112        allocate(htop(ngrid,naerkind),htop2(ngrid,naerkind))
113        allocate(nmode(ngrid,naerkind))
114        call cloud_haus16(ngrid,pbot2,pbot,hbot2,hbot,ptop,ptop2,htop,htop2,nmode)
115
116        firstcall=.false.
117      ENDIF ! of IF (firstcall)
118! ---------------------------------------------------------
119! ---------------------------------------------------------
120
121    aerosol(:,:,:) = 0.
122        tau_col(:) = 0.
123
124      do ig=1,ngrid
125        do iaer=1,naerkind
126
127!       Opacity calculation
128
129! determine the approximate middle of the mode layer
130           ll=1
131           DO l=1,nlayer-1    ! high p to low p
132             IF (pplay(ig,l) .gt. (ptop(ig,iaer)+pbot(ig,iaer))/2) ll=l
133           ENDDO
134! from there go DOWN for profile N
135           mode_dens = nmode(ig,iaer)  ! m-3
136           DO l=ll,1,-1
137             hlay=RD*(pt(ig,l+1)+pt(ig,l))/2./RG
138             IF (pplay(ig,l) .le. pbot(ig,iaer)) THEN
139               mode_dens = nmode(ig,iaer)  ! m-3
140             ELSEIF (pplay(ig,l) .gt. pbot(ig,iaer) .and. pplay(ig,l) .le. pbot2(ig,iaer)) THEN
141               mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(hlay/hbot(ig,iaer))
142             ELSEIF (pplay(ig,l) .gt. pbot2(ig,iaer)) THEN
143               mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(hlay/hbot2(ig,iaer))
144             ENDIF
145             mode_dens=max(1.e-30,mode_dens)
146             if (iaer.lt.5) then
147              aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
148               RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* &
149               mode_dens*hlay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
150             else  ! for UV absorber
151! normalized to 0.35 microns (peak of absorption)
152              aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*mode_dens
153             endif
154           ENDDO
155! ... then UP for profile N
156           mode_dens = nmode(ig,iaer)  ! m-3
157           DO l=ll+1,nlayer-1
158             hlay=RD*(pt(ig,l-1)+pt(ig,l))/2./RG
159             IF (pplay(ig,l) .gt. ptop(ig,iaer)) THEN
160               mode_dens = nmode(ig,iaer)  ! m-3
161             ELSEIF (pplay(ig,l) .le. ptop(ig,iaer) .and. pplay(ig,l) .gt. ptop2(ig,iaer)) THEN
162               mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(hlay/htop(ig,iaer))
163             ELSEIF (pplay(ig,l) .le. ptop2(ig,iaer)) THEN
164               mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(hlay/htop2(ig,iaer))
165             ENDIF
166             mode_dens=max(1.e-30,mode_dens)
167             if (iaer.lt.5) then
168              aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
169               RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* &
170               mode_dens*hlay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
171             else  ! for UV absorber
172! normalized to 0.35 microns (peak of absorption)
173              aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*mode_dens
174             endif
175           ENDDO
176
177           if (iaer.eq.5) then  ! for UV absorber
178            DO l=1,nlayer
179               tau_col(ig) = tau_col(ig) &
180                     + aerosol(ig,l,iaer)
181            ENDDO
182            DO l=1,nlayer
183                aerosol(ig,l,iaer)=aerosol(ig,l,iaer)*0.205/tau_col(ig)
184            ENDDO
185           endif
186
187        enddo
188      enddo
189
190!==================================================================
191!     if (is_master) then
192!      ig=10
193!      do l=1,nlayer
194!         write(82,*) l,pplay(ig,l),aerosol(ig,l,1)
195!         write(83,*) l,pplay(ig,l),aerosol(ig,l,2)
196!         write(84,*) l,pplay(ig,l),aerosol(ig,l,3)
197!         write(85,*) l,pplay(ig,l),aerosol(ig,l,4)
198!         write(86,*) l,pplay(ig,l),aerosol(ig,l,5)
199!      enddo
200!     endif   
201!      stop         
202!==================================================================
203
204! --------------------------------------------------------------------------
205! Column integrated visible optical depth in each point (used for diagnostic)
206
207      tau_col(:)=0.0
208      do ig=1,ngrid
209         do l=1,nlayer
210            do iaer = 1, naerkind
211               tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
212            end do
213         end do
214      end do
215
216      do ig=1,ngrid
217         do l=1,nlayer
218            do iaer = 1, naerkind
219               if(aerosol(ig,l,iaer).gt.1.e3)then
220                  print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
221                  print*,'at ig=',ig,',  l=',l,', iaer=',iaer
222                  print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
223                  print*,'reffrad=',reffrad(ig,l,iaer)
224               endif
225            end do
226         end do
227      end do
228
229      do ig=1,ngrid
230         if(tau_col(ig).gt.1.e3)then
231            print*,'WARNING: tau_col=',tau_col(ig)
232            print*,'at ig=',ig
233            print*,'aerosol=',aerosol(ig,:,:)
234            print*,'QREFvis3d=',QREFvis3d(ig,:,:)
235            print*,'reffrad=',reffrad(ig,:,:)
236         endif
237      end do
238
239      return
240    end subroutine aeropacity
241
242! --------------------------------------------------------------------------
243! --------------------------------------------------------------------------
244subroutine cloud_haus16(ngrid,pbot2,pbot,hbot2,hbot,ptop,ptop2,htop,htop2,nmode)
245
246!==================================================================
247!         Venus clouds (4 modes)
248!   S. Lebonnois (jan 2016)
249!==================================================================
250! distributions from Haus et al, 2016
251! ------------------------------------
252!
253! mode             1      2      2p     3
254! r (microns)     0.30   1.05   1.40   3.65
255! sigma           1.56   1.29   1.23   1.28
256! reff (microns)  0.49   1.23   1.56   4.25
257! nueff           0.21   0.067  0.044  0.063
258! (nueff=exp(ln^2 sigma)-1)
259!
260! pbot <=> zb ; ptop <=> zb+zc ; hbot <=> Hlo ; htop <=> Hup
261! p<ptop: N(i+1)=N(i)*(p(i+1)/p(i))**(hlay(i)/htop)     
262!    hlay=R(T(i)+T(i+1))/2/g  (in m)
263!    R=8.31/mu (mu in kg/mol)
264! p>pbot: N(i-1)=N(i)*(p(i)/p(i-1))**(hlay(i)/hbot)     
265! N is in m-3
266!
267! values in each mode below from Table 1 and text (p.186)
268!
269! dTau = Qext*[pi*reff**2*exp(-3*ln(1+nueff))]*N*hlay*(-dp)/p
270!
271! Latitudinal dependence (values from tables 3 and 4):
272! mode factor MF12(lat) for modes 1 and 2
273! mode 2' unchanged
274! mode factor MF3(lat) for mode 3
275! + for mode 2 only : pbot(lat), ptop(lat), htop(lat), htop2(lat)
276!==================================================================
277
278      use radinc_h, only : naerkind
279      USE geometry_mod, ONLY: latitude_deg             
280      implicit none
281       
282      INTEGER,INTENT(IN) :: ngrid  ! number of atmospheric columns
283      real,INTENT(out) :: pbot2(ngrid,naerkind),pbot(ngrid,naerkind)
284      real,INTENT(out) :: ptop(ngrid,naerkind),ptop2(ngrid,naerkind)
285      real,INTENT(out) :: hbot2(ngrid,naerkind),hbot(ngrid,naerkind)
286      real,INTENT(out) :: htop(ngrid,naerkind),htop2(ngrid,naerkind)
287      real,INTENT(out) :: nmode(ngrid,naerkind)
288     
289! For latitudinal dependence
290      integer :: i,j,ilat
291      real :: latit,factlat
292      integer,parameter :: nlat_H16 = 16
293      real :: lat_H16(nlat_H16)
294      real :: MF12(nlat_H16),MF3(nlat_H16)
295      real :: pbotM2(nlat_H16),ptopM2(nlat_H16),htopM2(nlat_H16)
296
297!-------------------------
298!!   Equatorial model 
299!-------------------------
300! initialization
301      pbot2(:,:) = 1.e8
302      pbot(:,:)  = 1.e8
303      ptop2(:,:) = 0.
304      ptop(:,:)  = 0.
305      hbot2(:,:) = 1.
306      hbot(:,:)  = 1.
307      htop2(:,:) = 1.
308      htop(:,:)  = 1.
309      nmode(:,:) = 0.
310       
311! table of values(lat,mode)
312
313! Mode 1
314      pbot2(:,1) = 2.0e5 ! Pa
315      hbot2(:,1) = 5.0e3 ! m
316      pbot(:,1)  = 1.2e5
317      hbot(:,1)  = 1.0e3
318      nmode(:,1) = 1.935e8 ! m-3
319      ptop(:,1)  = 1.0e4
320      htop(:,1)  = 3.5e3
321      ptop2(:,1) = 4.5e2
322      htop2(:,1) = 2.0e3
323
324! Mode 2
325      pbot(:,2)  = 1.0e4
326      hbot(:,2)  = 3.0e3
327      nmode(:,2) = 1.00e8 ! m-3
328      ptop(:,2)  = 8.0e3
329      htop(:,2)  = 3.5e3
330      ptop2(:,2) = 4.5e2
331      htop2(:,2) = 2.0e3
332
333! Mode 2p
334      pbot(:,3)  = 1.2e5
335      hbot(:,3)  = 0.1e3
336      nmode(:,3) = 5.0e7 ! m-3
337      ptop(:,3)  = 2.3e4
338      htop(:,3)  = 1.0e3
339
340! Mode 3
341      pbot(:,4)  = 1.2e5
342      hbot(:,4)  = 0.5e3
343      nmode(:,4) = 1.4e7 ! m-3
344      ptop(:,4)  = 3.9e4
345      htop(:,4)  = 1.0e3
346
347! UV absorber
348      pbot(:,5)  = 3.3e4  ! 58 km
349      hbot(:,5)  = 1.0e3
350      nmode(:,5) = 1.00e7 !m-3
351      ptop(:,5)  = 3.7e3  ! 70 km
352      htop(:,5)  = 1.0e3
353
354!----------------------------
355!!   Latitudinal variations 
356!----------------------------
357      lat_H16(:) = (/0.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,90./)
358      MF12(:) = (/0.98,0.98,0.99,1.00,0.98,0.94,0.86,0.81,0.73,0.67,0.64,0.61,0.59,0.47,0.36,0.36/)
359      MF3(:)  = (/1.30,1.30,1.26,1.23,1.17,1.13,1.06,1.03,1.04,1.09,1.22,1.51,1.82,2.02,2.09,2.09/)
360      pbotM2(:) = (/1.0e4,1.0e4,1.0e4,1.0e4,1.0e4,1.0e4,1.0e4,1.0e4,9.4e3, &
361                    9.2e3,9.1e3,9.6e3,1.14e4,1.21e4,1.25e4,1.25e4/)
362      ptopM2(:) = (/8.0e3,8.0e3,8.0e3,8.0e3,8.0e3,8.0e3,8.0e3,8.0e3,7.7e3, &
363                    7.5e3,7.5e3,8.0e3,9.2e3,1.0e4,1.06e4,1.06e4/)
364      htopM2(:) = (/3.5e3,3.5e3,3.5e3,3.5e3,3.5e3,3.5e3,3.5e3,3.5e3,3.4e3, &
365                    3.2e3,2.6e3,2.0e3,2.0e3,0.6e3,0.5e3,0.5e3/)
366     
367      do i=1,ngrid
368        latit = abs(latitude_deg(i))
369        ilat=2
370        do j=2,nlat_H16
371           if (latit.gt.lat_H16(j-1)) ilat=j
372        enddo
373        factlat    = (lat_H16(ilat)-latit)/(lat_H16(ilat)-lat_H16(ilat-1))
374        pbot(i,2)  = pbotM2(ilat)*(1.-factlat) + pbotM2(ilat-1)*factlat
375        ptop(i,2)  = ptopM2(ilat)*(1.-factlat) + ptopM2(ilat-1)*factlat
376        htop(i,2)  = htopM2(ilat)*(1.-factlat) + htopM2(ilat-1)*factlat
377        htop2(i,2) = min(htop2(i,2),htop(i,2))
378        nmode(i,1) = nmode(i,1)*(MF12(ilat)*(1.-factlat)+MF12(ilat-1)*factlat)
379        nmode(i,2) = nmode(i,2)*(MF12(ilat)*(1.-factlat)+MF12(ilat-1)*factlat)
380        nmode(i,3) = nmode(i,3)*(MF3(ilat) *(1.-factlat)+MF3(ilat-1) *factlat)
381      enddo
382     
383      return
384end subroutine cloud_haus16
Note: See TracBrowser for help on using the repository browser.