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 | ! -------------------------------------------------------------------------- |
---|
244 | subroutine 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 |
---|
384 | end subroutine cloud_haus16 |
---|