Ignore:
Timestamp:
May 17, 2021, 3:35:58 PM (3 years ago)
Author:
evignon
Message:

Commission de la nouvelle interface LMDZ-SISVAT
la prochaine commission consistera a supprimer l'ancien repertoire sisvat
et a faire un peu de nettoyage.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_sop.F

    r3792 r3900  
    1  
    2  
    31      subroutine SnOptP(jjtime)
    42 
    … …  
    3028C |   OUTPUT:  albisv   : Snow/Ice/Water/Soil Integrated Albedo        [-] |
    3129C |   ^^^^^^   sEX_sv   : Verticaly Integrated Extinction Coefficient      |
     30C |            DOPsnSV   : Snow Optical diameter [m]                       |
    3231C |                                                                        |
    3332C |   Internal Variables:                                                  |
    … …  
    6867      use VARySV
    6968      use VARtSV
    70       USE surface_data, only: iflag_albzenith
     69      USE surface_data, only: iflag_albcalc,correc_alb
    7170
    7271      IMPLICIT NONE
    … …  
    8988c #AG real      agesno
    9089 
    91       integer   isn   ,ikl   ,isn1 
     90      integer   isn   ,ikl   ,isn1, i 
    9291      real      sbeta1,sbeta2,sbeta3,sbeta4,sbeta5
    9392      real      AgeMax,AlbMin,HSnoSV,HIceSV,doptmx,SignG1,Sph_OK
    9493      real      dalbed,dalbeS,dalbeW
    95       real      bsegal,czemax,csegal
     94      real      bsegal,czemax,csegal,csza
    9695      real      RoFrez,DiffRo,SignRo,SnowOK,OpSqrt
    9796      real      albSn1,albIc1,a_SnI1,a_SII1
    … …  
    103102      real      albedo_old,albCor
    104103      real      ro_ave,dz_ave,minalb
    105       real      thetazs,thetazs0,aap,bbp,ccp,alb0      ! parameters for zenoth angle effect at Dome C
    106 
    107 
     104      real      l1min,l1max,l2min,l2max,l3min,l3max
     105      real      l6min(6), l6max(6), albSn6(6), a_SII6(6)
     106      real      lmintmp,lmaxtmp,albtmp
    108107 
    109108C +--Local   DATA
    … …  
    138137      data      albmax  /0.99       /  ! Albedo max
    139138
    140 C +-- For the computation of solar zenoth angle effect from Dome C data
    141 C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    142       data      aap/0.00016/,bbp/-0.017/,ccp/1.2/
    143 
     139C +-- wavelength limits [m] for each broad band
     140
     141      data      l1min/400.0e-9/,l1max/800.0e-9/
     142      data      l2min/800.0e-9/,l2max/1500.0e-9/
     143      data      l3min/1500.0e-9/,l3max/2800.0e-9/
     144
     145      data      l6min/185.0e-9,250.0e-9,400.0e-9,
     146     .               690.0e-9,1190.0e-9,2380.0e-9/
     147      data      l6max/250.0e-9,400.0e-9,
     148     .          690.0e-9,1190.0e-9,2380.0e-9,4000.0e-9/
    144149 
    145150 
    … …  
    177182     .                                       *DFcdSV                 )))
    178183          SnOpSV(ikl,isn) =  max(zero,SnOpSV(ikl,isn))
     184
     185C + --For outputs (Etienne)
     186C + ------------------------
     187          DOPsnSV(ikl,isn)=SnOpSV(ikl,isn)
    179188        END DO
    180189      END DO
    181190 
     191
     192
    182193 
    183194C +--Snow/Ice Albedo
    … …  
    191202        DO ikl=1,knonv
    192203 
    193            isn   =  max(iun,isnoSV(ikl))
     204          isn   =  max(iun,isnoSV(ikl))
    194205 
    195206          SignRo = sign(unun, rocdSV - ro__SV(ikl,isn))
    … …  
    200211cCA    +--Correction of snow albedo for Antarctica/Greenland
    201212cCA       --------------------------------------------------
    202           albCor = 1.
     213
     214         
     215          albCor = correc_alb
    203216c #GL     albCor = 1.01
    204 c #AC     albCor = 1.01
    205  
     217c #AC    albCor = 1.01
     218
     219
     220          IF (iflag_albcalc .GE. 1) THEN  ! Albedo calculation according to Kokhanovsky and Zege 2004
     221
     222          dalbed = 0.0
     223          doptic=SnOpSV(ikl,isn)
     224          csza=coszSV(ikl)
     225
     226          CALL albedo_kokhanovsky(l1min,l1max,csza,doptic,albSn1)
     227          CALL albedo_kokhanovsky(l2min,l2max,csza,doptic,albSn2)
     228          CALL albedo_kokhanovsky(l3min,l3max,csza,doptic,albSn3)
     229
     230          DO i=1,6
     231             lmintmp=l6min(i)
     232             lmaxtmp=l6max(i)
     233             CALL albedo_kokhanovsky(lmintmp,lmaxtmp,csza,doptic,albtmp)
     234             albSn6(i)=albtmp
     235          ENDDO
     236
     237
     238          ELSE ! Default calculation in SISVAT
     239
     240!    Zenith Angle Correction (Segal et al.,         1991, JAS 48, p.1025)
     241!--------------------------- (Wiscombe & Warren, dec1980, JAS   , p.2723)
     242!                            (Warren,               1982,  RG   , p.  81)
     243!                            -------------------------------------------------
     244         
     245          dalbed = 0.0
     246
     247          csegal = max(czemax                   ,coszSV(ikl))
     248c #cz     dalbeS =   ((bsegal+unun)/(unun+2.0*bsegal*csegal)
     249c #cz.                -       unun                          )*0.32
     250c #cz.              /  bsegal
     251c #cz     dalbeS = max(dalbeS,zero)
     252c #cz     dalbed =     dalbeS      *       min(1,isnoSV(ikl))
     253 
     254          dalbeW =(0.64 - csegal  )*0.0625  ! Warren 1982, RevGeo, fig.12b
     255                                            ! 0.0625 = 5% * 1/0.8,   p.81
     256                                            ! 0.64   = cos(50)
     257          dalbed =     dalbeW      *       min(1,isnoSV(ikl))
     258!-------------------------------------------------------------------------
     259
    206260          albSn1 =  0.96-1.580*OpSqrt
    207261          albSn1 =  max(albSn1,AlbMin)
    … …  
    219273          albSn3 =  min(albSn3*albCor,unun)
    220274 
     275          albSn6(1:3)=albSn1
     276          albSn6(4:6)=albSn2
     277
    221278          !snow albedo corection if wetsnow
    222279c #GL     albSn1 =  albSn1*max(0.9,(1.-1.5*eta_SV(ikl,isn)))
    223280c #GL     albSn2 =  albSn2*max(0.9,(1.-1.5*eta_SV(ikl,isn)))
    224281c #GL     albSn3 =  albSn3*max(0.9,(1.-1.5*eta_SV(ikl,isn)))
     282
     283          ENDIF
     284
    225285 
    226286          albSno =  So1dSV*albSn1
    … …  
    237297          albSn2 =  SnowOK*albSn2+(1.0-SnowOK)*max(albSno,minalb)
    238298          albSn3 =  SnowOK*albSn3+(1.0-SnowOK)*max(albSno,minalb)
     299          albSn6(:) =  SnowOK*albSn6(:)+(1.0-SnowOK)*max(albSno,minalb)
     300
    239301 
    240302c +           ro < roSdSV |          min al > aI3dSV
    … …  
    304366          a_SII3      =     albWIc      +(albSn3-albWIc)     *SnownH
    305367          a_SII3      = min(a_SII3       ,albSn3)
    306  
     368
     369          DO i=1,6
     370          a_SII6(i)   = albWIc      +(albSn6(i)-albWIc)     *SnownH
     371          a_SII6(i)   = min(a_SII6(i)       ,albSn6(i))
     372          ENDDO
     373
    307374cc #AG     agesno =      min(agsnSV(ikl,isn)          ,AgeMax)
    308375cc #AG     a_SII1      =     a_SII1      -0.175*agesno/AgeMax
    309376C +...                                   Impurities: Col de Porte Parameter.
    310377 
    311  
    312 !    Zenith Angle Correction (Segal et al.,         1991, JAS 48, p.1025)
    313 !    ----------------------- (Wiscombe & Warren, dec1980, JAS   , p.2723)
    314 !                            (Warren,               1982,  RG   , p.  81)
    315 !                            (Vignon, PhD Thesis, pp 150, conditions at Dome C)
    316 !                            -------------------------------------------------
    317  
    318          
    319           dalbed = 0.0
    320 
    321           if (iflag_albzenith .GE. 0) then
    322           csegal = max(czemax                   ,coszSV(ikl))
    323 c #cz     dalbeS =   ((bsegal+unun)/(unun+2.0*bsegal*csegal)
    324 c #cz.                -       unun                          )*0.32
    325 c #cz.              /  bsegal
    326 c #cz     dalbeS = max(dalbeS,zero)
    327 c #cz     dalbed =     dalbeS      *       min(1,isnoSV(ikl))
    328  
    329           dalbeW =(0.64 - csegal  )*0.0625  ! Warren 1982, RevGeo, fig.12b
    330                                             ! 0.0625 = 5% * 1/0.8,   p.81
    331                                             ! 0.64   = cos(50)
    332           dalbed =     dalbeW      *       min(1,isnoSV(ikl))
    333 
    334 
    335 ! Vignon PhD thesis, Dome C conditions
    336 
    337             if (iflag_albzenith .EQ. 1) then
    338             thetazs0=-bbp/(2.0*aap)
    339             alb0=(bbp**2)/4./aap-(bbp**2)/2./aap+ccp
    340             thetazs=max(acos(coszSV(ikl))/pi*180.0,thetazs0) ! in degrees
    341 
    342 
    343             dalbeW = aap*(thetazs**2)+bbp*thetazs+ccp-alb0       
    344             dalbed =     dalbeW      *       min(1,isnoSV(ikl))
    345             end if
    346 
    347 
    348             end if
     378
    349379 
    350380C +--Elsewhere    Integrated Snow/Ice Albedo
    … …  
    368398            alb3sv(ikl) =     albssv(ikl) +(a_SII3-albssv(ikl))*SIcenH
    369399            alb3sv(ikl) = min(alb3sv(ikl)  ,a_SII3)
    370  
     400
    371401            albisv(ikl) =     albssv(ikl) +(albSII-albssv(ikl))*SIcenH
    372402            albisv(ikl) = min(albisv(ikl)  ,albSII)
    373  
     403
     404            DO i=1,6
     405            alb6sv(ikl,i) = albssv(ikl) +(a_SII6(i)-albssv(ikl))*SIcenH
     406            alb6sv(ikl,i) = min(alb6sv(ikl,i)  ,a_SII6(i))
     407            ENDDO
     408
    374409 
    375410C +--Integrated Snow/Ice/Soil Albedo: Clouds Correction! Greuell & all., 1994
    … …  
    382417            alb3sv(ikl) = alb3sv(ikl) + 0.05 *(cld_SV(ikl)-0.5)*SIcenH
    383418     .                  + dalbed      *    (1.-cld_SV(ikl))
     419            alb6sv(ikl,:) =alb6sv(ikl,:)+0.05 *(cld_SV(ikl)-0.5)*SIcenH
     420     .                  + dalbed      *    (1.-cld_SV(ikl))
    384421            albisv(ikl) = albisv(ikl) + 0.05 *(cld_SV(ikl)-0.5)*SIcenH
    385422     .                  + dalbed      *    (1.-cld_SV(ikl))
    … …  
    397434            alb3sv(ikl) = alb3sv(ikl) - 1.0/3.0             ! 33 %
    398435     .                  * (albedo_old-albisv(ikl)) / So3dSV
    399  
    400  
    401 C +--Integrated Snow/Ice/Soil Albedo: Maximum albedo = 95%
     436            alb6sv(ikl,1:3) = alb6sv(ikl,1:3) - 1.0/6.0     ! 16 %
     437     .                  * (albedo_old-albisv(ikl)) / (So1dSV/3)
     438            alb6sv(ikl,4:6) = alb6sv(ikl,4:6) - 1.0/6.0     ! 16 %
     439     .                  * (albedo_old-albisv(ikl)) / (So2dSV/3)
     440
     441
     442C +--Integrated Snow/Ice/Soil Albedo: Maximum albedo = 95%
    402443C +  -----------------------------------------------------
    403444 
    … …  
    410451            alb3sv(ikl) = alb3sv(ikl) - 1.0/3.0             ! 33 %
    411452     .                  * (albedo_old-albisv(ikl)) / So3dSV
     453            alb6sv(ikl,1:3) = alb6sv(ikl,1:3) - 1.0/6.0     ! 16 %
     454     .                  * (albedo_old-albisv(ikl)) / (So1dSV/3)
     455            alb6sv(ikl,4:6) = alb6sv(ikl,4:6) - 1.0/6.0     ! 16 %
     456     .                  * (albedo_old-albisv(ikl)) / (So2dSV/3)
    412457
    413458
    … …  
    436481            alb2sv(ikl) = min(max(zero,alb2sv(ikl)),albmax)
    437482            alb3sv(ikl) = min(max(zero,alb3sv(ikl)),albmax)
    438  
     483           
     484            DO i=1,6
     485                alb6sv(ikl,i) = min(max(zero,alb6sv(ikl,i)),albmax)
     486            ENDDO
    439487        END DO
    440488 
    … …  
    519567 
    520568      return
     569
     570
    521571      end
     572
     573
     574!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     575      SUBROUTINE albedo_kokhanovsky(lambdamin,lambdamax,
     576     .                              cossza,dopt,albint)
     577!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     578! Authors: Hajar El Habchi El Fenniri, Etienne Vignon, Cecile Agosta
     579!          Ghislain Picard
     580! Routine that calculates the  integrated snow spectral albedo between
     581! lambdamin and lambdamax following Kokhanisky and Zege 2004,
     582! Scattering optics of snow, Applied Optics, Vol 43, No7
     583! Code inspired from the snowoptics package of Ghislain Picard:
     584! https://github.com/ghislainp/snowoptics
     585!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     586
     587       
     588      USE VARphy
     589
     590      IMPLICIT NONE
     591
     592! Inputs
     593!--------
     594      REAL lambdamin   ! minimum wavelength for integration [m]
     595      REAL lambdamax   ! maximum wavelength for integration [m]
     596      REAL cossza     ! solar zenith angle cosinus
     597      REAL dopt        ! optical diameter [m]
     598
     599!Outputs
     600!-------
     601      REAL albint
     602
     603! Local Variables
     604!-----------------
     605
     606      REAL ropt,cosalb,norm,Pas
     607      REAL SSA,alpha,gamm,R,cos30,alb30
     608      INTEGER i
     609
     610
     611      REAL B_amp                       ! amplification factor
     612      PARAMETER (B_amp=1.6) 
     613
     614      REAL g_asy                       ! asymetry factor
     615      PARAMETER (g_asy=0.845)
     616
     617      INTEGER nlambda                  ! length of wavelength vector
     618      PARAMETER(nlambda=200)
     619
     620      REAL lmin
     621      PARAMETER(lmin=185.0E-9)
     622
     623      REAL lmax
     624      PARAMETER(lmax=4000.0E-9)
     625
     626      REAL albmax
     627      PARAMETER(albmax=0.99)
     628
     629      REAL wavelengths(nlambda)
     630      REAL ni(nlambda)
     631
     632      DATA wavelengths / 1.85000000e-07, 2.04170854e-07,
     633     . 2.23341709e-07, 2.42512563e-07,
     634     . 2.61683417e-07, 2.80854271e-07, 3.00025126e-07, 3.19195980e-07,
     635     . 3.38366834e-07, 3.57537688e-07, 3.76708543e-07, 3.95879397e-07,
     636     . 4.15050251e-07, 4.34221106e-07, 4.53391960e-07, 4.72562814e-07,
     637     . 4.91733668e-07, 5.10904523e-07, 5.30075377e-07, 5.49246231e-07,
     638     . 5.68417085e-07, 5.87587940e-07, 6.06758794e-07, 6.25929648e-07,
     639     . 6.45100503e-07, 6.64271357e-07, 6.83442211e-07, 7.02613065e-07,
     640     . 7.21783920e-07, 7.40954774e-07, 7.60125628e-07, 7.79296482e-07,
     641     . 7.98467337e-07, 8.17638191e-07, 8.36809045e-07, 8.55979899e-07,
     642     . 8.75150754e-07, 8.94321608e-07, 9.13492462e-07, 9.32663317e-07,
     643     . 9.51834171e-07, 9.71005025e-07, 9.90175879e-07, 1.00934673e-06,
     644     . 1.02851759e-06, 1.04768844e-06, 1.06685930e-06, 1.08603015e-06,
     645     . 1.10520101e-06, 1.12437186e-06, 1.14354271e-06, 1.16271357e-06,
     646     . 1.18188442e-06, 1.20105528e-06, 1.22022613e-06, 1.23939698e-06,
     647     . 1.25856784e-06, 1.27773869e-06, 1.29690955e-06, 1.31608040e-06,
     648     . 1.33525126e-06, 1.35442211e-06, 1.37359296e-06, 1.39276382e-06,
     649     . 1.41193467e-06, 1.43110553e-06, 1.45027638e-06, 1.46944724e-06,
     650     . 1.48861809e-06, 1.50778894e-06, 1.52695980e-06, 1.54613065e-06,
     651     . 1.56530151e-06, 1.58447236e-06, 1.60364322e-06, 1.62281407e-06,
     652     . 1.64198492e-06, 1.66115578e-06, 1.68032663e-06, 1.69949749e-06,
     653     . 1.71866834e-06, 1.73783920e-06, 1.75701005e-06, 1.77618090e-06,
     654     . 1.79535176e-06, 1.81452261e-06, 1.83369347e-06, 1.85286432e-06,
     655     . 1.87203518e-06, 1.89120603e-06, 1.91037688e-06, 1.92954774e-06,
     656     . 1.94871859e-06, 1.96788945e-06, 1.98706030e-06, 2.00623116e-06,
     657     . 2.02540201e-06, 2.04457286e-06, 2.06374372e-06, 2.08291457e-06,
     658     . 2.10208543e-06, 2.12125628e-06, 2.14042714e-06, 2.15959799e-06,
     659     . 2.17876884e-06, 2.19793970e-06, 2.21711055e-06, 2.23628141e-06,
     660     . 2.25545226e-06, 2.27462312e-06, 2.29379397e-06, 2.31296482e-06,
     661     . 2.33213568e-06, 2.35130653e-06, 2.37047739e-06, 2.38964824e-06,
     662     . 2.40881910e-06, 2.42798995e-06, 2.44716080e-06, 2.46633166e-06,
     663     . 2.48550251e-06, 2.50467337e-06, 2.52384422e-06, 2.54301508e-06,
     664     . 2.56218593e-06, 2.58135678e-06, 2.60052764e-06, 2.61969849e-06,
     665     . 2.63886935e-06, 2.65804020e-06, 2.67721106e-06, 2.69638191e-06,
     666     . 2.71555276e-06, 2.73472362e-06, 2.75389447e-06, 2.77306533e-06,
     667     . 2.79223618e-06, 2.81140704e-06, 2.83057789e-06, 2.84974874e-06,
     668     . 2.86891960e-06, 2.88809045e-06, 2.90726131e-06, 2.92643216e-06,
     669     . 2.94560302e-06, 2.96477387e-06, 2.98394472e-06, 3.00311558e-06,
     670     . 3.02228643e-06, 3.04145729e-06, 3.06062814e-06, 3.07979899e-06,
     671     . 3.09896985e-06, 3.11814070e-06, 3.13731156e-06, 3.15648241e-06,
     672     . 3.17565327e-06, 3.19482412e-06, 3.21399497e-06, 3.23316583e-06,
     673     . 3.25233668e-06, 3.27150754e-06, 3.29067839e-06, 3.30984925e-06,
     674     . 3.32902010e-06, 3.34819095e-06, 3.36736181e-06, 3.38653266e-06,
     675     . 3.40570352e-06, 3.42487437e-06, 3.44404523e-06, 3.46321608e-06,
     676     . 3.48238693e-06, 3.50155779e-06, 3.52072864e-06, 3.53989950e-06,
     677     . 3.55907035e-06, 3.57824121e-06, 3.59741206e-06, 3.61658291e-06,
     678     . 3.63575377e-06, 3.65492462e-06, 3.67409548e-06, 3.69326633e-06,
     679     . 3.71243719e-06, 3.73160804e-06, 3.75077889e-06, 3.76994975e-06,
     680     . 3.78912060e-06, 3.80829146e-06, 3.82746231e-06, 3.84663317e-06,
     681     . 3.86580402e-06, 3.88497487e-06, 3.90414573e-06, 3.92331658e-06,
     682     . 3.94248744e-06, 3.96165829e-06, 3.98082915e-06, 4.00000000e-06/
     683
     684
     685      DATA ni /7.74508407e-10, 7.74508407e-10,
     686     .  7.74508407e-10, 7.74508407e-10,
     687     .  7.74508407e-10, 7.74508407e-10, 7.74508407e-10, 7.74508407e-10,
     688     .  6.98381122e-10, 6.23170274e-10, 5.97655992e-10, 5.84106004e-10,
     689     .  5.44327597e-10, 5.71923510e-10, 6.59723827e-10, 8.05183870e-10,
     690     .  1.03110161e-09, 1.36680386e-09, 1.85161253e-09, 2.56487751e-09,
     691     .  3.56462855e-09, 4.89450926e-09, 6.49252022e-09, 9.62029335e-09,
     692     .  1.32335015e-08, 1.75502184e-08, 2.19240625e-08, 3.03304156e-08,
     693     .  4.07715972e-08, 5.00414911e-08, 7.09722331e-08, 1.00773751e-07,
     694     .  1.31427409e-07, 1.42289041e-07, 1.49066787e-07, 2.01558515e-07,
     695     .  2.99106105e-07, 4.03902086e-07, 4.54292169e-07, 5.21906983e-07,
     696     .  6.27643362e-07, 9.43955678e-07, 1.33464494e-06, 1.97278315e-06,
     697     .  2.31801329e-06, 2.20584676e-06, 1.85568138e-06, 1.73395193e-06,
     698     .  1.73101406e-06, 1.91333936e-06, 2.26413019e-06, 3.23959718e-06,
     699     .  4.94316963e-06, 6.89378896e-06, 1.02237444e-05, 1.21439656e-05,
     700     .  1.31567585e-05, 1.33448288e-05, 1.32000000e-05, 1.31608040e-05,
     701     .  1.33048369e-05, 1.40321464e-05, 1.51526244e-05, 1.80342858e-05,
     702     .  3.82875736e-05, 1.07325259e-04, 2.11961637e-04, 3.82008054e-04,
     703     .  5.30897470e-04, 5.29244735e-04, 4.90876605e-04, 4.33905427e-04,
     704     .  3.77795349e-04, 3.17633099e-04, 2.81078564e-04, 2.57579485e-04,
     705     .  2.42203100e-04, 2.23789060e-04, 2.04306870e-04, 1.87909255e-04,
     706     .  1.73117146e-04, 1.61533186e-04, 1.53420328e-04, 1.47578033e-04,
     707     .  1.42334776e-04, 1.35691466e-04, 1.30495414e-04, 1.36065123e-04,
     708     .  1.70928821e-04, 2.66389730e-04, 4.80957955e-04, 8.25041961e-04,
     709     .  1.21654792e-03, 1.50232875e-03, 1.62316078e-03, 1.61649750e-03,
     710     .  1.53736801e-03, 1.42343711e-03, 1.24459117e-03, 1.02388611e-03,
     711     .  7.89112523e-04, 5.97204264e-04, 4.57152413e-04, 3.62341259e-04,
     712     .  2.99128332e-04, 2.57035569e-04, 2.26992203e-04, 2.07110355e-04,
     713     .  2.05835688e-04, 2.25108810e-04, 2.64262893e-04, 3.23594011e-04,
     714     .  3.93061117e-04, 4.62789970e-04, 5.19664416e-04, 5.59739628e-04,
     715     .  5.93476084e-04, 6.22797885e-04, 6.57484833e-04, 6.92849600e-04,
     716     .  7.26584901e-04, 7.56604648e-04, 7.68009488e-04, 7.65579073e-04,
     717     .  7.50526164e-04, 7.39809972e-04, 7.55622847e-04, 8.05099514e-04,
     718     .  9.67279246e-04, 1.16281559e-03, 1.42570247e-03, 2.04986585e-03,
     719     .  2.93971170e-03, 4.49827711e-03, 7.32537532e-03, 1.18889734e-02,
     720     .  1.85851805e-02, 2.86242532e-02, 4.34131035e-02, 6.37828307e-02,
     721     .  9.24145850e-02, 1.35547945e-01, 1.94143245e-01, 2.54542814e-01,
     722     .  3.02282024e-01, 3.42214181e-01, 3.85475065e-01, 4.38000000e-01,
     723     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
     724     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
     725     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
     726     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
     727     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
     728     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
     729     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
     730     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
     731     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
     732     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
     733     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
     734     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
     735     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01/
     736
     737
     738      Pas     =(lmax-lmin)/nlambda
     739      ropt    = dopt/2.0
     740      SSA     = 3.0/(rhoIce*ropt)
     741      cos30   = cos(30.0/180.0*pi)
     742
     743
     744      albint=0.
     745      norm=0.
     746       
     747      DO i = 1,nlambda
     748          gamm = 4.0 * pi * ni(i) / wavelengths(i)
     749          cosalb = 2.0 / (SSA * rhoice) * B_amp * gamm
     750          alpha = 16. / 3 * cosalb / (1.0 - g_asy)
     751          R = exp(-(alpha**0.5) * 3.0 / 7.0 * (1.0 + 2.0 * cossza))
     752          alb30 = exp(-(alpha**0.5)* 3.0 / 7.0 * (1.0 + 20 * cos30))
     753
     754          IF ((wavelengths(i).GE.lambdamin).AND.
     755     .       (wavelengths(i).LT.lambdamax)) THEN
     756          albint=albint+R*Pas  ! rectangle integration -> can be improved with trapezintegration
     757          norm=norm+Pas
     758          ENDIF
     759
     760      END DO
     761
     762      albint=max(0.,min(albint/max(norm,1E-10),albmax))
     763
     764      END
     765 
Note: See TracChangeset for help on using the changeset viewer.