source: LMDZ4/trunk/libf/phylmd/aeropt.F @ 674

Last change on this file since 674 was 524, checked in by lmdzadmin, 21 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.8 KB
Line 
1!
2! $Header$
3!
4      SUBROUTINE aeropt(pplay, paprs, t_seri, msulfate, RHcl,
5     .            tau_ae, piz_ae, cg_ae, ai        )
6c
7      IMPLICIT none
8c
9c
10c     
11#include "dimensions.h"
12#include "dimphy.h"
13#include "YOMCST.h"
14c
15c Arguments:
16c
17      REAL paprs(klon,klev+1)
18      REAL pplay(klon,klev), t_seri(klon,klev)
19      REAL msulfate(klon,klev) ! masse sulfate ug SO4/m3  [ug/m^3]
20      REAL RHcl(klon,klev)     ! humidite relative ciel clair
21      REAL tau_ae(klon,klev,2) ! epaisseur optique aerosol
22      REAL piz_ae(klon,klev,2) ! single scattering albedo aerosol
23      REAL cg_ae(klon,klev,2)  ! asymmetry parameter aerosol
24      REAL ai(klon)            ! POLDER aerosol index
25c
26c Local
27c
28      INTEGER i, k, inu
29      INTEGER RH_num, nbre_RH
30      PARAMETER (nbre_RH=12)
31      REAL RH_tab(nbre_RH)
32      REAL RH_MAX, DELTA, rh
33      PARAMETER (RH_MAX=95.)
34      DATA RH_tab/0.,10.,20.,30.,40.,50.,60.,70.,80.,85.,90.,95./
35      REAL zrho, zdz
36      REAL taue670(klon)       ! epaisseur optique aerosol absorption 550 nm
37      REAL taue865(klon)       ! epaisseur optique aerosol extinction 865 nm
38      REAL alpha_aer_sulfate(nbre_RH,5)   !--unit m2/g SO4
39      REAL alphasulfate     
40c
41c Proprietes optiques
42c
43      REAL alpha_aer(nbre_RH,2)   !--unit m2/g SO4
44      REAL cg_aer(nbre_RH,2)
45      DATA alpha_aer/.500130E+01,  .500130E+01,  .500130E+01, 
46     .               .500130E+01,  .500130E+01,  .616710E+01, 
47     .               .826850E+01,  .107687E+02,  .136976E+02, 
48     .               .162972E+02,  .211690E+02,  .354833E+02, 
49     .               .139460E+01,  .139460E+01,  .139460E+01, 
50     .               .139460E+01,  .139460E+01,  .173910E+01, 
51     .               .244380E+01,  .332320E+01,  .440120E+01, 
52     .               .539570E+01,  .734580E+01,  .136038E+02 /
53      DATA cg_aer/.619800E+00,  .619800E+00,  .619800E+00, 
54     .            .619800E+00,  .619800E+00,  .662700E+00, 
55     .            .682100E+00,  .698500E+00,  .712500E+00, 
56     .            .721800E+00,  .734600E+00,  .755800E+00, 
57     .            .545600E+00,  .545600E+00,  .545600E+00, 
58     .            .545600E+00,  .545600E+00,  .583700E+00, 
59     .            .607100E+00,  .627700E+00,  .645800E+00, 
60     .            .658400E+00,  .676500E+00,  .708500E+00 /
61      DATA alpha_aer_sulfate/
62     . 4.910,4.910,4.910,4.910,6.547,7.373,
63     . 8.373,9.788,12.167,14.256,17.924,28.433,
64     . 1.453,1.453,1.453,1.453,2.003,2.321,
65     . 2.711,3.282,4.287,5.210,6.914,12.305,
66     . 4.308,4.308,4.308,4.308,5.753,6.521,
67     . 7.449,8.772,11.014,12.999,16.518,26.772,
68     . 3.265,3.265,3.265,3.265,4.388,5.016,
69     . 5.775,6.868,8.745,10.429,13.457,22.538,
70     . 2.116,2.116,2.116,2.116,2.882,3.330,
71     . 3.876,4.670,6.059,7.327,9.650,16.883/
72c
73      DO i=1, klon
74         taue670(i)=0.0
75         taue865(i)=0.0
76      ENDDO
77c     
78      DO k=1, klev
79      DO i=1, klon
80         if (t_seri(i,k).eq.0) write (*,*) 'aeropt T ',i,k,t_seri(i,k)
81         if (pplay(i,k).eq.0) write (*,*) 'aeropt p ',i,k,pplay(i,k)         
82        zrho=pplay(i,k)/t_seri(i,k)/RD                  ! kg/m3
83        zdz=(paprs(i,k)-paprs(i,k+1))/zrho/RG           ! m
84        rh=MIN(RHcl(i,k)*100.,RH_MAX)
85        RH_num = INT( rh/10. + 1.)
86        IF (rh.LT.0.) STOP 'aeropt: RH < 0 not possible'
87        IF (rh.gt.85.) RH_num=10
88        IF (rh.gt.90.) RH_num=11
89        DELTA=(rh-RH_tab(RH_num))/(RH_tab(RH_num+1)-RH_tab(RH_num))
90c               
91        inu=1
92        tau_ae(i,k,inu)=alpha_aer(RH_num,inu) +
93     .             DELTA*(alpha_aer(RH_num+1,inu)-alpha_aer(RH_num,inu))
94        tau_ae(i,k,inu)=tau_ae(i,k,inu)*msulfate(i,k)*zdz*1.e-6
95        piz_ae(i,k,inu)=1.0
96        cg_ae(i,k,inu)=cg_aer(RH_num,inu) +
97     .                 DELTA*(cg_aer(RH_num+1,inu)-cg_aer(RH_num,inu))
98c
99        inu=2
100        tau_ae(i,k,inu)=alpha_aer(RH_num,inu) +
101     .             DELTA*(alpha_aer(RH_num+1,inu)-alpha_aer(RH_num,inu))
102        tau_ae(i,k,inu)=tau_ae(i,k,inu)*msulfate(i,k)*zdz*1.e-6
103        piz_ae(i,k,inu)=1.0
104        cg_ae(i,k,inu)=cg_aer(RH_num,inu) +
105     .                 DELTA*(cg_aer(RH_num+1,inu)-cg_aer(RH_num,inu))
106cjq
107cjq for aerosol index
108c
109        alphasulfate=alpha_aer_sulfate(RH_num,4) +
110     .       DELTA*(alpha_aer_sulfate(RH_num+1,4)-
111     .       alpha_aer_sulfate(RH_num,4)) !--m2/g
112c     
113        taue670(i)=taue670(i)+
114     .       alphasulfate*msulfate(i,k)*zdz*1.e-6
115c
116        alphasulfate=alpha_aer_sulfate(RH_num,5) +
117     .       DELTA*(alpha_aer_sulfate(RH_num+1,5)-
118     .       alpha_aer_sulfate(RH_num,5)) !--m2/g
119c
120        taue865(i)=taue865(i)+
121     .         alphasulfate*msulfate(i,k)*zdz*1.e-6
122       
123      ENDDO
124      ENDDO
125c     
126      DO i=1, klon
127        ai(i)=(-log(MAX(taue670(i),0.0001)/
128     .                MAX(taue865(i),0.0001))/log(670./865.)) *
129     .        taue865(i)
130      ENDDO     
131c
132      RETURN
133      END
Note: See TracBrowser for help on using the repository browser.