source: TOOLS/CMIP6_FORCING/AER_OPTICS/refr_ind_seasalt_MACC_v4.f @ 3385

Last change on this file since 3385 was 3385, checked in by oboucher, 4 years ago

Uploading package to prepare tropospheric aerosol optics LUT in LMDz fo rCMIP6

File size: 3.8 KB
Line 
1      PROGRAM sea_salt
2!--based on Tang and Munkelwitz JGR 1994 and Tang JGR 1997
3!--then correct for observed growth factor
4      IMPLICIT NONE
5      REAL rhop, rhow, nrw, niw, nrp, nip, xx
6      PARAMETER (rhop=2.16,rhow=1.0) !--g cm-3
7      INTEGER nbre_RH, Nmode
8      PARAMETER(nbre_rh=12, Nmode=3)
9      REAL RH_tab(nbre_RH)
10      DATA RH_tab/0.,10.,20.,30.,40.,50.,60.,70.,80.,85.,90.,95./
11c--observed growth factor SS - CS - AS
12      REAL growth(nbre_RH,Nmode)
13      DATA growth/
14     .  1.00,1.11,1.15,1.20,1.25,1.31,1.39,1.51,1.70,1.85,2.10,2.61, !--SS
15     .  1.00,1.10,1.14,1.18,1.23,1.29,1.37,1.48,1.66,1.81,2.05,2.55, !--CS
16     .  1.00,1.09,1.13,1.17,1.21,1.27,1.34,1.44,1.61,1.75,1.98,2.45/ !--AS
17c
18c--wavelength in m
19      INTEGER Nwvmax, Nwvhale, Nwvdry, wve, wv
20      PARAMETER (Nwvmax=200,Nwvhale=60,Nwvdry=61)
21      REAL lambda(Nwvmax)
22      REAL lambda_hale(Nwvhale), n_r_hale(Nwvhale), n_i_hale(Nwvhale)
23      REAL lambda_dry(Nwvdry),   n_r_dry(Nwvdry),  n_i_dry(Nwvdry)
24c
25      INTEGER IRH, mode
26      REAL lambda_min, lambda_max
27      PARAMETER (lambda_min=0.2E-6, lambda_max=5.0E-6)
28c
29      REAL n_r(nbre_rh,Nmode,Nwvmax)
30      REAL n_i(nbre_rh,Nmode,Nwvmax)
31c
32      REAL ratio(nbre_rh,Nmode)      !--diameter growth factor
33      REAL rhod(nbre_rh,Nmode)       !--density
34c
35      DO wve=1,Nwvmax
36        lambda(wve)=
37     .  lambda_min+FLOAT(wve-1)*(lambda_max-lambda_min)/FLOAT(Nwvmax-1)
38      ENDDO 
39c
40c--reading water refractive index
41      OPEN(unit=11,file='r_water_hale.dat')
42      DO wve=1,Nwvhale
43        READ(11,*) lambda_hale(wve), n_r_hale(wve)
44      ENDDO
45      CLOSE(11)
46c      print *,'n_r_hale=', n_r_hale
47c
48      OPEN(unit=11,file='i_water_hale.dat')
49      DO wve=1,Nwvhale
50        READ(11,*) lambda_hale(wve), n_i_hale(wve)
51      ENDDO
52      CLOSE(11)
53c      print *,'n_i_hale=', n_i_hale
54c
55c--reading dry seasalt refractive index
56      OPEN(unit=11,file='r_seasal_she.dat')
57      DO wve=1,Nwvdry
58        READ(11,*) lambda_dry(wve), n_r_dry(wve)
59      ENDDO
60      CLOSE(11)
61c      print *,'n_r_dry=', n_r_dry
62c
63      OPEN(unit=11,file='i_seasal_she.dat')
64      DO wve=1,Nwvdry
65        READ(11,*) lambda_dry(wve), n_i_dry(wve)
66      ENDDO
67      CLOSE(11)
68c      print *,'n_i_dry=', n_i_dry
69c
70      OPEN (unit=11,file='growth_ss_SS_v2')
71      OPEN (unit=12,file='growth_ss_CS_v2')
72      OPEN (unit=13,file='growth_ss_AS_v2')
73c
74      OPEN (unit=21,file='ri_ss_SS_v2')
75      OPEN (unit=22,file='ri_ss_AS_v2')
76      OPEN (unit=23,file='ri_ss_CS_v2')
77c
78      DO mode=1, Nmode
79c
80      DO IRH=1, nbre_RH
81c
82      xx=1./growth(IRH,mode)**3.
83      rhod(IRH,mode)=xx*rhop+(1.-xx)*rhow 
84c
85      DO wve=1,Nwvmax
86c
87c-plus proche voisin
88          nrw=n_r_hale(1)
89          niw=n_i_hale(1)
90          DO wv=1, Nwvhale
91            IF (lambda_hale(wv)/1.e9.LT.lambda(wve)) THEN
92              nrw=n_r_hale(wv)
93              niw=n_i_hale(wv)
94            ENDIF 
95          ENDDO
96
97          nrp=n_r_dry(1)
98          nip=n_i_dry(1)
99          DO wv=1, Nwvdry
100            IF (lambda_dry(wv)/1.e9.LT.lambda(wve)) THEN
101              nrp=n_r_dry(wv)
102              nip=n_i_dry(wv)
103            ENDIF 
104          ENDDO
105c
106          n_r(IRH,mode,wve) = xx*nrp +(1.-xx)*nrw 
107          n_i(IRH,mode,wve) = xx*nip +(1.-xx)*niw 
108c
109          WRITE(20+mode,*) RH_tab(IRH), lambda(wve),
110     .                     n_r(IRH,mode,wve), n_i(IRH,mode,wve)
111c
112      ENDDO !-wve
113      ENDDO !-IRH
114c
115      WRITE(10+mode,20)'RH', (RH_tab(IRH),IRH=1,nbre_rh)
116      WRITE(10+mode,10)"growth factor",(growth(IRH,mode),IRH=1,nbre_rh)
117      WRITE(10+mode,10)"density factor",
118     .            (rhod(IRH,mode)/rhod(1,mode),IRH=1,nbre_rh)
119      WRITE(10+mode,10)'rho @ RH', (rhod(IRH,mode),IRH=1,nbre_rh)
120c
121      ENDDO
122
12310    FORMAT(A14, F6.3,11(',',F6.3))
12420    FORMAT(A14, F6.1,11(',',F6.1))
12563    FORMAT(F10.6,E13.6)
126 
127      CLOSE(11)
128      CLOSE(12)
129      CLOSE(13)
130
131      CLOSE(21)
132      CLOSE(22)
133      CLOSE(23)
134c
135      END
Note: See TracBrowser for help on using the repository browser.