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

      PROGRAM sea_salt
!--based on Tang and Munkelwitz JGR 1994 and Tang JGR 1997 
!--then correct for observed growth factor
      IMPLICIT NONE 
      REAL rhop, rhow, nrw, niw, nrp, nip, xx
      PARAMETER (rhop=2.16,rhow=1.0) !--g cm-3
      INTEGER nbre_RH, Nmode
      PARAMETER(nbre_rh=12, Nmode=3)
      REAL RH_tab(nbre_RH)
      DATA RH_tab/0.,10.,20.,30.,40.,50.,60.,70.,80.,85.,90.,95./
c--observed growth factor SS - CS - AS
      REAL growth(nbre_RH,Nmode)
      DATA growth/
     .  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
     .  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
     .  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
c
c--wavelength in m
      INTEGER Nwvmax, Nwvhale, Nwvdry, wve, wv
      PARAMETER (Nwvmax=200,Nwvhale=60,Nwvdry=61)
      REAL lambda(Nwvmax)
      REAL lambda_hale(Nwvhale), n_r_hale(Nwvhale), n_i_hale(Nwvhale)
      REAL lambda_dry(Nwvdry),   n_r_dry(Nwvdry),  n_i_dry(Nwvdry)
c
      INTEGER IRH, mode
      REAL lambda_min, lambda_max
      PARAMETER (lambda_min=0.2E-6, lambda_max=5.0E-6)
c
      REAL n_r(nbre_rh,Nmode,Nwvmax)
      REAL n_i(nbre_rh,Nmode,Nwvmax)
c
      REAL ratio(nbre_rh,Nmode)      !--diameter growth factor
      REAL rhod(nbre_rh,Nmode)       !--density
c
      DO wve=1,Nwvmax
        lambda(wve)=
     .  lambda_min+FLOAT(wve-1)*(lambda_max-lambda_min)/FLOAT(Nwvmax-1)
      ENDDO 
c
c--reading water refractive index
      OPEN(unit=11,file='r_water_hale.dat')
      DO wve=1,Nwvhale
        READ(11,*) lambda_hale(wve), n_r_hale(wve)
      ENDDO
      CLOSE(11)
c      print *,'n_r_hale=', n_r_hale
c
      OPEN(unit=11,file='i_water_hale.dat')
      DO wve=1,Nwvhale
        READ(11,*) lambda_hale(wve), n_i_hale(wve)
      ENDDO
      CLOSE(11)
c      print *,'n_i_hale=', n_i_hale
c
c--reading dry seasalt refractive index
      OPEN(unit=11,file='r_seasal_she.dat')
      DO wve=1,Nwvdry
        READ(11,*) lambda_dry(wve), n_r_dry(wve)
      ENDDO
      CLOSE(11)
c      print *,'n_r_dry=', n_r_dry
c
      OPEN(unit=11,file='i_seasal_she.dat')
      DO wve=1,Nwvdry
        READ(11,*) lambda_dry(wve), n_i_dry(wve)
      ENDDO
      CLOSE(11)
c      print *,'n_i_dry=', n_i_dry
c
      OPEN (unit=11,file='growth_ss_SS_v2')
      OPEN (unit=12,file='growth_ss_CS_v2')
      OPEN (unit=13,file='growth_ss_AS_v2')
c
      OPEN (unit=21,file='ri_ss_SS_v2')
      OPEN (unit=22,file='ri_ss_AS_v2')
      OPEN (unit=23,file='ri_ss_CS_v2')
c
      DO mode=1, Nmode
c
      DO IRH=1, nbre_RH
c
      xx=1./growth(IRH,mode)**3.
      rhod(IRH,mode)=xx*rhop+(1.-xx)*rhow 
c
      DO wve=1,Nwvmax
c
c-plus proche voisin
          nrw=n_r_hale(1)
          niw=n_i_hale(1)
          DO wv=1, Nwvhale
            IF (lambda_hale(wv)/1.e9.LT.lambda(wve)) THEN
              nrw=n_r_hale(wv)
              niw=n_i_hale(wv)
            ENDIF  
          ENDDO

          nrp=n_r_dry(1)
          nip=n_i_dry(1)
          DO wv=1, Nwvdry
            IF (lambda_dry(wv)/1.e9.LT.lambda(wve)) THEN
              nrp=n_r_dry(wv)
              nip=n_i_dry(wv)
            ENDIF  
          ENDDO
c
          n_r(IRH,mode,wve) = xx*nrp +(1.-xx)*nrw 
          n_i(IRH,mode,wve) = xx*nip +(1.-xx)*niw 
c
          WRITE(20+mode,*) RH_tab(IRH), lambda(wve),
     .                     n_r(IRH,mode,wve), n_i(IRH,mode,wve)
c
      ENDDO !-wve
      ENDDO !-IRH
c
      WRITE(10+mode,20)'RH', (RH_tab(IRH),IRH=1,nbre_rh)
      WRITE(10+mode,10)"growth factor",(growth(IRH,mode),IRH=1,nbre_rh)
      WRITE(10+mode,10)"density factor",
     .            (rhod(IRH,mode)/rhod(1,mode),IRH=1,nbre_rh)
      WRITE(10+mode,10)'rho @ RH', (rhod(IRH,mode),IRH=1,nbre_rh)
c
      ENDDO

10    FORMAT(A14, F6.3,11(',',F6.3))
20    FORMAT(A14, F6.1,11(',',F6.1))
63    FORMAT(F10.6,E13.6)
 
      CLOSE(11)
      CLOSE(12)
      CLOSE(13)

      CLOSE(21)
      CLOSE(22)
      CLOSE(23)
c
      END