subroutine SnOptP(jjtime) ! +------------------------------------------------------------------------+ ! | MAR/SISVAT SnOptP 12-08-2019 MAR | ! | SubRoutine SnOptP computes the Snow Pack optical Properties | ! +------------------------------------------------------------------------+ ! | | ! | PARAMETERS: knonv: Total Number of columns = | ! | ^^^^^^^^^^ = Total Number of continental Grid Boxes | ! | X Number of Mosaic Cell per Grid Box | ! | | ! | INPUT: isnoSV = total Nb of Ice/Snow Layers | ! | ^^^^^ ispiSV = 0,...,nsno: Uppermost Superimposed Ice Layer | ! | | ! | | ! | INPUT: G1snSV : Dendricity (<0) or Sphericity (>0) of Snow Layer | ! | ^^^^^ G2snSV : Sphericity (>0) or Size of Snow Layer | ! | agsnSV : Snow Age [day] | ! | ro__SV : Snow/Soil Volumic Mass [kg/m3] | ! | eta_SV : Water Content [m3/m3] | ! | rusnSV : Surficial Water Thickness [kg/m2] .OR. [mm] | ! | SWS_SV : Surficial Water Status | ! | dzsnSV : Snow Layer Thickness [m] | ! | | ! | albssv : Soil Albedo [-] | ! | zzsnsv : Snow Pack Thickness [m] | ! | | ! | OUTPUT: albisv : Snow/Ice/Water/Soil Integrated Albedo [-] | ! | ^^^^^^ sEX_sv : Verticaly Integrated Extinction Coefficient | ! | DOPsnSV : Snow Optical diameter [m] | ! | | ! | Internal Variables: | ! | ^^^^^^^^^^^^^^^^^^ | ! | SnOpSV : Snow Grain optical Size [m] | ! | EX1_sv : Integrated Snow Extinction (0.3--0.8micr.m) | ! | EX2_sv : Integrated Snow Extinction (0.8--1.5micr.m) | ! | EX3_sv : Integrated Snow Extinction (1.5--2.8micr.m) | ! | | ! | METHODE: Calcul de la taille optique des grains ? partir de | ! | ^^^^^^^ -leur type decrit par les deux variables descriptives | ! | continues sur la plage -99/+99 passees en appel. | ! | -la taille optique (1/10mm) des etoiles, | ! | des grains fins et | ! | des jeunes faces planes | ! | | ! | METHOD: Computation of the optical diameter of the grains | ! | ^^^^^^ described with the CROCUS formalism G1snSV / G2snSV | ! | | ! | REFERENCE: Brun et al. 1989, J. Glaciol 35 pp. 333--342 | ! | ^^^^^^^^^ Brun et al. 1992, J. Glaciol 38 pp. 13-- 22 | ! | Eric Martin Sept.1996 | ! | | ! | | ! +------------------------------------------------------------------------+ ! +--Global Variables ! + ================ use VARphy use VAR_SV use VARdSV use VARxSV use VARySV use VARtSV USE surface_data, only: iflag_albcalc,correc_alb IMPLICIT NONE ! + -- INPUT integer :: jjtime ! +--Internal Variables ! + ================== real :: coalb1(knonv) ! weighted Coalbedo, Vis. real :: coalb2(knonv) ! weighted Coalbedo, nIR 1 real :: coalb3(knonv) ! weighted Coalbedo, nIR 2 real :: coalbm ! weighted Coalbedo, mean real :: sExt_1(knonv) ! Extinction Coeff., Vis. real :: sExt_2(knonv) ! Extinction Coeff., nIR 1 real :: sExt_3(knonv) ! Extinction Coeff., nIR 2 real :: SnOpSV(knonv, nsno) ! Snow Grain optical Size ! #AG real agesno integer :: isn ,ikl ,isn1, i real :: sbeta1,sbeta2,sbeta3,sbeta4,sbeta5 real :: AgeMax,AlbMin,HSnoSV,HIceSV,doptmx,SignG1,Sph_OK real :: dalbed,dalbeS,dalbeW real :: bsegal,czemax,csegal,csza real :: RoFrez,DiffRo,SignRo,SnowOK,OpSqrt real :: albSn1,albIc1,a_SnI1,a_SII1 real :: albSn2,albIc2,a_SnI2,a_SII2 real :: albSn3,albIc3,a_SnI3,a_SII3 real :: albSno,albIce,albSnI,albSII,albWIc,albmax real :: doptic,Snow_H,SIce_H,SnownH,SIcenH real :: exarg1,exarg2,exarg3,sign_0,sExt_0 real :: albedo_old,albCor real :: ro_ave,dz_ave,minalb real :: l1min,l1max,l2min,l2max,l3min,l3max real :: l6min(6), l6max(6), albSn6(6), a_SII6(6) real :: lmintmp,lmaxtmp,albtmp ! +--Local DATA ! + ============ ! +--For the computation of the solar irradiance extinction in snow ! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ data sbeta1/0.0192/,sbeta2/0.4000/,sbeta3/0.1098/ data sbeta4/1.0000/ data sbeta5/2.00e1/ ! +--Snow Age Maximum (Taiga, e.g. Col de Porte) ! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ data AgeMax /60.0/ ! +... AgeMax: Snow Age Maximum [day] data AlbMin /0.94/ ! +... AlbMin: Albedo Minimum / visible (0.3--0.8 micrometers) data HSnoSV /0.01/ ! +... HSnoSV: Snow Thickness through witch ! + Albedo is interpolated to Ice Albedo data HIceSV /0.10/ ! +... HIceSV: Snow/Ice Thickness through witch ! + Albedo is interpolated to Soil Albedo data doptmx /2.3e-3/ ! +... doptmx: Maximum optical Diameter (pi * R**2) [m] ! + data czeMAX /0.173648178/ ! 80.deg (Segal et al., 1991 JAS) data bsegal /4.00 / ! data albmax /0.99 / ! Albedo max ! +-- wavelength limits [m] for each broad band data l1min/400.0e-9/,l1max/800.0e-9/ data l2min/800.0e-9/,l2max/1500.0e-9/ data l3min/1500.0e-9/,l3max/2800.0e-9/ data l6min/185.0e-9,250.0e-9,400.0e-9, & 690.0e-9,1190.0e-9,2380.0e-9/ data l6max/250.0e-9,400.0e-9, & 690.0e-9,1190.0e-9,2380.0e-9,4000.0e-9/ ! +--Snow Grain optical Size ! + ======================= DO ikl=1,knonv DO isn=1,max(1,isnoSV(ikl)) G2snSV(ikl,isn) = max(epsi,G2snSV(ikl,isn)) ! +... Avoid non physical Values SignG1 = sign(unun,G1snSV(ikl,isn)) Sph_OK = max(zero,SignG1) SnOpSV(ikl,isn) = 1.e-4 * & ! +... SI: (from 1/10 mm to m) ! +--Contribution of Non Dendritic Snow ! + ---------------------------------- ( Sph_OK *( G2snSV(ikl,isn)*G1snSV(ikl,isn)/G1_dSV & +max(demi*G2snSV(ikl,isn),DFcdSV) & *(unun-G1snSV(ikl,isn) /G1_dSV)) & ! +--Contribution of Dendritic Snow ! + ---------------------------------- +(1.-Sph_OK)*( -G1snSV(ikl,isn)*DDcdSV /G1_dSV & +(unun+G1snSV(ikl,isn) /G1_dSV) & * (G2snSV(ikl,isn)*DScdSV /G1_dSV & +(unun-G2snSV(ikl,isn) /G1_dSV) & *DFcdSV ))) SnOpSV(ikl,isn) = max(zero,SnOpSV(ikl,isn)) ! + --For outputs (Etienne) ! + ------------------------ DOPsnSV(ikl,isn)=SnOpSV(ikl,isn) END DO END DO ! +--Snow/Ice Albedo ! + =============== ! +--Uppermost effective Snow Layer ! + ------------------------------ DO ikl=1,knonv isn = max(iun,isnoSV(ikl)) SignRo = sign(unun, rocdSV - ro__SV(ikl,isn)) SnowOK = max(zero,SignRo) ! Ice Density Threshold OpSqrt = sqrt(SnOpSV(ikl,isn)) !CA +--Correction of snow albedo for Antarctica/Greenland !CA -------------------------------------------------- albCor = correc_alb ! #GL albCor = 1.01 ! #AC albCor = 1.01 IF (iflag_albcalc .GE. 1) THEN ! Albedo calculation according to Kokhanovsky and Zege 2004 dalbed = 0.0 doptic=SnOpSV(ikl,isn) csza=coszSV(ikl) CALL albedo_kokhanovsky(l1min,l1max,csza,doptic,albSn1) CALL albedo_kokhanovsky(l2min,l2max,csza,doptic,albSn2) CALL albedo_kokhanovsky(l3min,l3max,csza,doptic,albSn3) DO i=1,6 lmintmp=l6min(i) lmaxtmp=l6max(i) CALL albedo_kokhanovsky(lmintmp,lmaxtmp,csza,doptic,albtmp) albSn6(i)=albtmp ENDDO ELSE ! Default calculation in SISVAT ! Zenith Angle Correction (Segal et al., 1991, JAS 48, p.1025) !--------------------------- (Wiscombe & Warren, dec1980, JAS , p.2723) ! (Warren, 1982, RG , p. 81) ! ------------------------------------------------- dalbed = 0.0 csegal = max(czemax ,coszSV(ikl)) ! #cz dalbeS = ((bsegal+unun)/(unun+2.0*bsegal*csegal) ! #cz. - unun )*0.32 ! #cz. / bsegal ! #cz dalbeS = max(dalbeS,zero) ! #cz dalbed = dalbeS * min(1,isnoSV(ikl)) dalbeW =(0.64 - csegal )*0.0625 ! Warren 1982, RevGeo, fig.12b ! ! 0.0625 = 5% * 1/0.8, p.81 ! ! 0.64 = cos(50) dalbed = dalbeW * min(1,isnoSV(ikl)) !------------------------------------------------------------------------- albSn1 = 0.96-1.580*OpSqrt albSn1 = max(albSn1,AlbMin) albSn1 = max(albSn1,zero) albSn1 = min(albSn1*albCor,unun) albSn2 = 0.95-15.40*OpSqrt albSn2 = max(albSn2,zero) albSn2 = min(albSn2*albCor,unun) doptic = min(SnOpSV(ikl,isn),doptmx) albSn3 = 346.3*doptic -32.31*OpSqrt +0.88 albSn3 = max(albSn3,zero) albSn3 = min(albSn3*albCor,unun) albSn6(1:3)=albSn1 albSn6(4:6)=albSn2 ! !snow albedo corection if wetsnow ! #GL albSn1 = albSn1*max(0.9,(1.-1.5*eta_SV(ikl,isn))) ! #GL albSn2 = albSn2*max(0.9,(1.-1.5*eta_SV(ikl,isn))) ! #GL albSn3 = albSn3*max(0.9,(1.-1.5*eta_SV(ikl,isn))) ENDIF albSno = So1dSV*albSn1 & + So2dSV*albSn2 & + So3dSV*albSn3 !XF minalb = (aI2dSV + (aI3dSV -aI2dSV) & * (ro__SV(ikl,isn)-ro_Ice)/(roSdSV-ro_Ice)) minalb = min(aI3dSV,max(aI2dSV,minalb)) ! pure/firn albedo SnowOK = SnowOK*max(zero,sign(unun, albSno-minalb)) albSn1 = SnowOK*albSn1+(1.0-SnowOK)*max(albSno,minalb) albSn2 = SnowOK*albSn2+(1.0-SnowOK)*max(albSno,minalb) albSn3 = SnowOK*albSn3+(1.0-SnowOK)*max(albSno,minalb) albSn6(:) = SnowOK*albSn6(:)+(1.0-SnowOK)*max(albSno,minalb) ! + ro < roSdSV | min al > aI3dSV ! + roSdSV < ro < rocdSV | aI2dSV < min al < aI3dSV (fct of density) ! +--Snow/Ice Pack Thickness ! + ----------------------- isn = max(min(isnoSV(ikl) ,ispiSV(ikl)),0) Snow_H = zzsnsv(ikl,isnoSV(ikl))-zzsnsv(ikl,isn) SIce_H = zzsnsv(ikl,isnoSV(ikl)) SnownH = Snow_H / HSnoSV SnownH = min(unun, SnownH) SIcenH = SIce_H / (HIceSV) SIcenH = min(unun, SIcenH) ! + The value of SnownH is set to 1 in case of ice lenses above ! + 1m of dry snow (ro<600kg/m3) for using CROCUS albedo ! ro_ave = 0. ! dz_ave = 0. ! SnowOK = 1. ! do isn = isnoSV(ikl),1,-1 ! ro_ave = ro_ave + ro__SV(ikl,isn) * dzsnSV(ikl,isn) * SnowOK ! dz_ave = dz_ave + dzsnSV(ikl,isn) * SnowOK ! SnowOK = max(zero,sign(unun,1.-dz_ave)) ! enddo ! ro_ave = ro_ave / max(dz_ave,epsi) ! SnowOK = max(zero,sign(unun,600.-ro_ave)) ! SnownH = SnowOK + SnownH * (1. - SnowOK) ! +--Integrated Snow/Ice Albedo: Case of Water on Bare Ice ! + ----------------------------------------------------- isn = max(min(isnoSV(ikl) ,ispiSV(ikl)),0) albWIc = aI1dSV-(aI1dSV-aI2dSV) & * exp(-(rusnSV(ikl) & ! * (1. -SWS_SV(ikl) & ! 0 <=> freezing * (1 -min(1,iabs(isn-isnoSV(ikl))))) & ! 1 <=> isn=isnoSV / ru_dSV)**0.50) ! ! albWIc = max(aI1dSV,min(aI2dSV,albWIc+slopSV(ikl)* ! . min(5.,max(1.,dx/10000.)))) SignRo = sign(unun,ro_Ice-5.-ro__SV(ikl,isn)) ! RoSN<920kg/m3 SnowOK = max(zero,SignRo) albWIc = (1. - SnowOK) * albWIc + SnowOK & * (aI2dSV + (aI3dSV -aI2dSV) & * (ro__SV(ikl,isn)-ro_Ice)/(roSdSV-ro_Ice)) ! + rocdSV < ro < ro_ice | aI2dSV< al <aI3dSV (fct of density) ! + ro > ro_ice | aI1dSV< al <aI2dSV (fct of superficial water content ! +--Integrated Snow/Ice Albedo ! + ------------------------------- a_SII1 = albWIc +(albSn1-albWIc) *SnownH a_SII1 = min(a_SII1 ,albSn1) a_SII2 = albWIc +(albSn2-albWIc) *SnownH a_SII2 = min(a_SII2 ,albSn2) a_SII3 = albWIc +(albSn3-albWIc) *SnownH a_SII3 = min(a_SII3 ,albSn3) DO i=1,6 a_SII6(i) = albWIc +(albSn6(i)-albWIc) *SnownH a_SII6(i) = min(a_SII6(i) ,albSn6(i)) ENDDO !c #AG agesno = min(agsnSV(ikl,isn) ,AgeMax) !c #AG a_SII1 = a_SII1 -0.175*agesno/AgeMax ! +... Impurities: Col de Porte Parameter. ! +--Elsewhere Integrated Snow/Ice Albedo ! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! #cp ELSE albSII = So1dSV*a_SII1 & + So2dSV*a_SII2 & + So3dSV*a_SII3 ! #cp END IF ! +--Integrated Snow/Ice/Soil Albedo ! + ------------------------------- alb1sv(ikl) = albssv(ikl) +(a_SII1-albssv(ikl))*SIcenH alb1sv(ikl) = min(alb1sv(ikl) ,a_SII1) alb2sv(ikl) = albssv(ikl) +(a_SII2-albssv(ikl))*SIcenH alb2sv(ikl) = min(alb2sv(ikl) ,a_SII2) alb3sv(ikl) = albssv(ikl) +(a_SII3-albssv(ikl))*SIcenH alb3sv(ikl) = min(alb3sv(ikl) ,a_SII3) albisv(ikl) = albssv(ikl) +(albSII-albssv(ikl))*SIcenH albisv(ikl) = min(albisv(ikl) ,albSII) DO i=1,6 alb6sv(ikl,i) = albssv(ikl) +(a_SII6(i)-albssv(ikl))*SIcenH alb6sv(ikl,i) = min(alb6sv(ikl,i) ,a_SII6(i)) ENDDO ! +--Integrated Snow/Ice/Soil Albedo: Clouds Correction! Greuell & all., 1994 ! + --------------------------------------------------! Glob.&t Planet.Change ! ! (9):91-114 alb1sv(ikl) = alb1sv(ikl) + 0.05 *(cld_SV(ikl)-0.5)*SIcenH & + dalbed * (1.-cld_SV(ikl)) alb2sv(ikl) = alb2sv(ikl) + 0.05 *(cld_SV(ikl)-0.5)*SIcenH & + dalbed * (1.-cld_SV(ikl)) alb3sv(ikl) = alb3sv(ikl) + 0.05 *(cld_SV(ikl)-0.5)*SIcenH & + dalbed * (1.-cld_SV(ikl)) alb6sv(ikl,:) =alb6sv(ikl,:)+0.05 *(cld_SV(ikl)-0.5)*SIcenH & + dalbed * (1.-cld_SV(ikl)) albisv(ikl) = albisv(ikl) + 0.05 *(cld_SV(ikl)-0.5)*SIcenH & + dalbed * (1.-cld_SV(ikl)) ! +--Integrated Snow/Ice/Soil Albedo: Minimum snow albedo = aI1dSV ! + ------------------------------------------------------------- albedo_old = albisv(ikl) albisv(ikl) = max(albisv(ikl),aI1dSV * SIcenH & + albssv(ikl) *(1.0 - SIcenH)) alb1sv(ikl) = alb1sv(ikl) - 1.0/3.0 & ! 33 % * (albedo_old-albisv(ikl)) / So1dSV alb2sv(ikl) = alb2sv(ikl) - 1.0/3.0 & ! 33 % * (albedo_old-albisv(ikl)) / So2dSV alb3sv(ikl) = alb3sv(ikl) - 1.0/3.0 & ! 33 % * (albedo_old-albisv(ikl)) / So3dSV alb6sv(ikl,1:3) = alb6sv(ikl,1:3) - 1.0/6.0 & ! 16 % * (albedo_old-albisv(ikl)) / (So1dSV/3) alb6sv(ikl,4:6) = alb6sv(ikl,4:6) - 1.0/6.0 & ! 16 % * (albedo_old-albisv(ikl)) / (So2dSV/3) ! +--Integrated Snow/Ice/Soil Albedo: Maximum albedo = 95% ! + ----------------------------------------------------- albedo_old = albisv(ikl) albisv(ikl) = min(albisv(ikl),0.95) alb1sv(ikl) = alb1sv(ikl) - 1.0/3.0 & ! 33 % * (albedo_old-albisv(ikl)) / So1dSV alb2sv(ikl) = alb2sv(ikl) - 1.0/3.0 & ! 33 % * (albedo_old-albisv(ikl)) / So2dSV alb3sv(ikl) = alb3sv(ikl) - 1.0/3.0 & ! 33 % * (albedo_old-albisv(ikl)) / So3dSV alb6sv(ikl,1:3) = alb6sv(ikl,1:3) - 1.0/6.0 & ! 16 % * (albedo_old-albisv(ikl)) / (So1dSV/3) alb6sv(ikl,4:6) = alb6sv(ikl,4:6) - 1.0/6.0 & ! 16 % * (albedo_old-albisv(ikl)) / (So2dSV/3) !Sea Ice/snow permanent-interractive prescription from Nemo !AO_CK 20/02/2020 ! ! No check if coupling update since MAR and NEMO albedo are too different ! and since MAR albedo is computed on properties that are not in NEMO ! ! prescription for each time step with NEMO values ! #AO if (LSmask(ikl) .eq. 0 .and. coupling_ao .eq. .true.) then ! #AO if (AOmask(ikl) .eq. 0) then ! #AO albisv(ikl) = (1.-AOmask(ikl))* albAOsisv(ikl) ! #AO. +(AOmask(ikl)*albisv(ikl)) ! #AO alb1sv(ikl) = (1.-AOmask(ikl))* albAOsisv(ikl) ! #AO. +(AOmask(ikl)*alb1sv(ikl)) ! #AO alb2sv(ikl) = (1.-AOmask(ikl))* albAOsisv(ikl) ! #AO. +(AOmask(ikl)*alb2sv(ikl)) ! #AO alb3sv(ikl) = (1.-AOmask(ikl))* albAOsisv(ikl) ! #AO. +(AOmask(ikl)*alb3sv(ikl)) ! #AO endif ! #AO endif alb1sv(ikl) = min(max(zero,alb1sv(ikl)),albmax) alb2sv(ikl) = min(max(zero,alb2sv(ikl)),albmax) alb3sv(ikl) = min(max(zero,alb3sv(ikl)),albmax) DO i=1,6 alb6sv(ikl,i) = min(max(zero,alb6sv(ikl,i)),albmax) ENDDO END DO ! +--Extinction Coefficient: Exponential Factor ! + ========================================== DO ikl=1,knonv sExt_1(ikl) = 1. sExt_2(ikl) = 1. sExt_3(ikl) = 1. sEX_sv(ikl,nsno+1) = 1. coalb1(ikl) = (1. -alb1sv(ikl))*So1dSV coalb2(ikl) = (1. -alb2sv(ikl))*So2dSV coalb3(ikl) = (1. -alb3sv(ikl))*So3dSV coalbm = coalb1(ikl) +coalb2(ikl) +coalb3(ikl) coalb1(ikl) = coalb1(ikl) /coalbm coalb2(ikl) = coalb2(ikl) /coalbm coalb3(ikl) = coalb3(ikl) /coalbm END DO !XF DO isn=nsno,1,-1 DO ikl=1,knonv sEX_sv(ikl,isn) = 1.0 ! !sEX_sv(ikl,isn) = 0.95 ! if MAR is too warm in summer END DO END DO DO ikl=1,knonv DO isn=max(1,isnoSV(ikl)),1,-1 SignRo = sign(unun, rocdSV - ro__SV(ikl,isn)) SnowOK = max(zero,SignRo) ! Ice Density Threshold RoFrez = 1.e-3 * ro__SV(ikl,isn) * (1.0-eta_SV(ikl,isn)) OpSqrt = sqrt(max(epsi,SnOpSV(ikl,isn))) exarg1 = SnowOK *1.e2 *max(sbeta1*RoFrez/OpSqrt,sbeta2) & +(1.0-SnowOK) *sbeta5 exarg2 = SnowOK *1.e2 *max(sbeta3*RoFrez/OpSqrt,sbeta4) & +(1.0-SnowOK) *sbeta5 exarg3 = SnowOK *1.e2 *sbeta5 & +(1.0-SnowOK) *sbeta5 ! +--Integrated Extinction of Solar Irradiance (Normalized Value) ! + ============================================================ sExt_1(ikl) = sExt_1(ikl) & * exp(min(0.0,-exarg1 *dzsnSV(ikl,isn))) sign_0 = sign(unun,eps9 -sExt_1(ikl)) sExt_0 = max(zero,sign_0)*sExt_1(ikl) sExt_1(ikl) = sExt_1(ikl) -sExt_0 sExt_2(ikl) = sExt_2(ikl) & * exp(min(0.0,-exarg2 *dzsnSV(ikl,isn))) sign_0 = sign(unun,eps9 -sExt_2(ikl)) sExt_0 = max(zero,sign_0)*sExt_2(ikl) sExt_2(ikl) = sExt_2(ikl) -sExt_0 sExt_3(ikl) = sExt_3(ikl) & * exp(min(0.0,-exarg3 *dzsnSV(ikl,isn))) sign_0 = sign(unun,eps9 -sExt_3(ikl)) sExt_0 = max(zero,sign_0)*sExt_3(ikl) sExt_3(ikl) = sExt_3(ikl) -sExt_0 sEX_sv(ikl,isn) = coalb1(ikl) * sExt_1(ikl) & + coalb2(ikl) * sExt_2(ikl) & + coalb3(ikl) * sExt_3(ikl) END DO END DO DO isn=0,-nsol,-1 DO ikl=1,knonv sEX_sv(ikl,isn) = 0.0 END DO END DO return end subroutine snoptp !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE albedo_kokhanovsky(lambdamin,lambdamax, & cossza,dopt,albint) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! Authors: Hajar El Habchi El Fenniri, Etienne Vignon, Cecile Agosta ! Ghislain Picard ! Routine that calculates the integrated snow spectral albedo between ! lambdamin and lambdamax following Kokhanisky and Zege 2004, ! Scattering optics of snow, Applied Optics, Vol 43, No7 ! Code inspired from the snowoptics package of Ghislain Picard: ! https://github.com/ghislainp/snowoptics !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ USE VARphy IMPLICIT NONE ! Inputs !-------- REAL :: lambdamin ! minimum wavelength for integration [m] REAL :: lambdamax ! maximum wavelength for integration [m] REAL :: cossza ! solar zenith angle cosinus REAL :: dopt ! optical diameter [m] !Outputs !------- REAL :: albint ! Local Variables !----------------- REAL :: ropt,cosalb,norm,Pas REAL :: SSA,alpha,gamm,R,cos30,alb30 INTEGER :: i REAL :: B_amp ! amplification factor PARAMETER (B_amp=1.6) REAL :: g_asy ! asymetry factor PARAMETER (g_asy=0.845) INTEGER :: nlambda ! length of wavelength vector PARAMETER(nlambda=200) REAL :: lmin PARAMETER(lmin=185.0E-9) REAL :: lmax PARAMETER(lmax=4000.0E-9) REAL :: albmax PARAMETER(albmax=0.99) REAL :: wavelengths(nlambda) REAL :: ni(nlambda) DATA wavelengths / 1.85000000e-07, 2.04170854e-07, & 2.23341709e-07, 2.42512563e-07, & 2.61683417e-07, 2.80854271e-07, 3.00025126e-07, 3.19195980e-07, & 3.38366834e-07, 3.57537688e-07, 3.76708543e-07, 3.95879397e-07, & 4.15050251e-07, 4.34221106e-07, 4.53391960e-07, 4.72562814e-07, & 4.91733668e-07, 5.10904523e-07, 5.30075377e-07, 5.49246231e-07, & 5.68417085e-07, 5.87587940e-07, 6.06758794e-07, 6.25929648e-07, & 6.45100503e-07, 6.64271357e-07, 6.83442211e-07, 7.02613065e-07, & 7.21783920e-07, 7.40954774e-07, 7.60125628e-07, 7.79296482e-07, & 7.98467337e-07, 8.17638191e-07, 8.36809045e-07, 8.55979899e-07, & 8.75150754e-07, 8.94321608e-07, 9.13492462e-07, 9.32663317e-07, & 9.51834171e-07, 9.71005025e-07, 9.90175879e-07, 1.00934673e-06, & 1.02851759e-06, 1.04768844e-06, 1.06685930e-06, 1.08603015e-06, & 1.10520101e-06, 1.12437186e-06, 1.14354271e-06, 1.16271357e-06, & 1.18188442e-06, 1.20105528e-06, 1.22022613e-06, 1.23939698e-06, & 1.25856784e-06, 1.27773869e-06, 1.29690955e-06, 1.31608040e-06, & 1.33525126e-06, 1.35442211e-06, 1.37359296e-06, 1.39276382e-06, & 1.41193467e-06, 1.43110553e-06, 1.45027638e-06, 1.46944724e-06, & 1.48861809e-06, 1.50778894e-06, 1.52695980e-06, 1.54613065e-06, & 1.56530151e-06, 1.58447236e-06, 1.60364322e-06, 1.62281407e-06, & 1.64198492e-06, 1.66115578e-06, 1.68032663e-06, 1.69949749e-06, & 1.71866834e-06, 1.73783920e-06, 1.75701005e-06, 1.77618090e-06, & 1.79535176e-06, 1.81452261e-06, 1.83369347e-06, 1.85286432e-06, & 1.87203518e-06, 1.89120603e-06, 1.91037688e-06, 1.92954774e-06, & 1.94871859e-06, 1.96788945e-06, 1.98706030e-06, 2.00623116e-06, & 2.02540201e-06, 2.04457286e-06, 2.06374372e-06, 2.08291457e-06, & 2.10208543e-06, 2.12125628e-06, 2.14042714e-06, 2.15959799e-06, & 2.17876884e-06, 2.19793970e-06, 2.21711055e-06, 2.23628141e-06, & 2.25545226e-06, 2.27462312e-06, 2.29379397e-06, 2.31296482e-06, & 2.33213568e-06, 2.35130653e-06, 2.37047739e-06, 2.38964824e-06, & 2.40881910e-06, 2.42798995e-06, 2.44716080e-06, 2.46633166e-06, & 2.48550251e-06, 2.50467337e-06, 2.52384422e-06, 2.54301508e-06, & 2.56218593e-06, 2.58135678e-06, 2.60052764e-06, 2.61969849e-06, & 2.63886935e-06, 2.65804020e-06, 2.67721106e-06, 2.69638191e-06, & 2.71555276e-06, 2.73472362e-06, 2.75389447e-06, 2.77306533e-06, & 2.79223618e-06, 2.81140704e-06, 2.83057789e-06, 2.84974874e-06, & 2.86891960e-06, 2.88809045e-06, 2.90726131e-06, 2.92643216e-06, & 2.94560302e-06, 2.96477387e-06, 2.98394472e-06, 3.00311558e-06, & 3.02228643e-06, 3.04145729e-06, 3.06062814e-06, 3.07979899e-06, & 3.09896985e-06, 3.11814070e-06, 3.13731156e-06, 3.15648241e-06, & 3.17565327e-06, 3.19482412e-06, 3.21399497e-06, 3.23316583e-06, & 3.25233668e-06, 3.27150754e-06, 3.29067839e-06, 3.30984925e-06, & 3.32902010e-06, 3.34819095e-06, 3.36736181e-06, 3.38653266e-06, & 3.40570352e-06, 3.42487437e-06, 3.44404523e-06, 3.46321608e-06, & 3.48238693e-06, 3.50155779e-06, 3.52072864e-06, 3.53989950e-06, & 3.55907035e-06, 3.57824121e-06, 3.59741206e-06, 3.61658291e-06, & 3.63575377e-06, 3.65492462e-06, 3.67409548e-06, 3.69326633e-06, & 3.71243719e-06, 3.73160804e-06, 3.75077889e-06, 3.76994975e-06, & 3.78912060e-06, 3.80829146e-06, 3.82746231e-06, 3.84663317e-06, & 3.86580402e-06, 3.88497487e-06, 3.90414573e-06, 3.92331658e-06, & 3.94248744e-06, 3.96165829e-06, 3.98082915e-06, 4.00000000e-06/ DATA ni /7.74508407e-10, 7.74508407e-10, & 7.74508407e-10, 7.74508407e-10, & 7.74508407e-10, 7.74508407e-10, 7.74508407e-10, 7.74508407e-10, & 6.98381122e-10, 6.23170274e-10, 5.97655992e-10, 5.84106004e-10, & 5.44327597e-10, 5.71923510e-10, 6.59723827e-10, 8.05183870e-10, & 1.03110161e-09, 1.36680386e-09, 1.85161253e-09, 2.56487751e-09, & 3.56462855e-09, 4.89450926e-09, 6.49252022e-09, 9.62029335e-09, & 1.32335015e-08, 1.75502184e-08, 2.19240625e-08, 3.03304156e-08, & 4.07715972e-08, 5.00414911e-08, 7.09722331e-08, 1.00773751e-07, & 1.31427409e-07, 1.42289041e-07, 1.49066787e-07, 2.01558515e-07, & 2.99106105e-07, 4.03902086e-07, 4.54292169e-07, 5.21906983e-07, & 6.27643362e-07, 9.43955678e-07, 1.33464494e-06, 1.97278315e-06, & 2.31801329e-06, 2.20584676e-06, 1.85568138e-06, 1.73395193e-06, & 1.73101406e-06, 1.91333936e-06, 2.26413019e-06, 3.23959718e-06, & 4.94316963e-06, 6.89378896e-06, 1.02237444e-05, 1.21439656e-05, & 1.31567585e-05, 1.33448288e-05, 1.32000000e-05, 1.31608040e-05, & 1.33048369e-05, 1.40321464e-05, 1.51526244e-05, 1.80342858e-05, & 3.82875736e-05, 1.07325259e-04, 2.11961637e-04, 3.82008054e-04, & 5.30897470e-04, 5.29244735e-04, 4.90876605e-04, 4.33905427e-04, & 3.77795349e-04, 3.17633099e-04, 2.81078564e-04, 2.57579485e-04, & 2.42203100e-04, 2.23789060e-04, 2.04306870e-04, 1.87909255e-04, & 1.73117146e-04, 1.61533186e-04, 1.53420328e-04, 1.47578033e-04, & 1.42334776e-04, 1.35691466e-04, 1.30495414e-04, 1.36065123e-04, & 1.70928821e-04, 2.66389730e-04, 4.80957955e-04, 8.25041961e-04, & 1.21654792e-03, 1.50232875e-03, 1.62316078e-03, 1.61649750e-03, & 1.53736801e-03, 1.42343711e-03, 1.24459117e-03, 1.02388611e-03, & 7.89112523e-04, 5.97204264e-04, 4.57152413e-04, 3.62341259e-04, & 2.99128332e-04, 2.57035569e-04, 2.26992203e-04, 2.07110355e-04, & 2.05835688e-04, 2.25108810e-04, 2.64262893e-04, 3.23594011e-04, & 3.93061117e-04, 4.62789970e-04, 5.19664416e-04, 5.59739628e-04, & 5.93476084e-04, 6.22797885e-04, 6.57484833e-04, 6.92849600e-04, & 7.26584901e-04, 7.56604648e-04, 7.68009488e-04, 7.65579073e-04, & 7.50526164e-04, 7.39809972e-04, 7.55622847e-04, 8.05099514e-04, & 9.67279246e-04, 1.16281559e-03, 1.42570247e-03, 2.04986585e-03, & 2.93971170e-03, 4.49827711e-03, 7.32537532e-03, 1.18889734e-02, & 1.85851805e-02, 2.86242532e-02, 4.34131035e-02, 6.37828307e-02, & 9.24145850e-02, 1.35547945e-01, 1.94143245e-01, 2.54542814e-01, & 3.02282024e-01, 3.42214181e-01, 3.85475065e-01, 4.38000000e-01, & 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, & 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, & 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, & 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, & 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, & 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, & 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, & 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, & 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, & 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, & 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, & 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, & 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01/ Pas =(lmax-lmin)/nlambda ropt = dopt/2.0 SSA = 3.0/(rhoIce*ropt) cos30 = cos(30.0/180.0*pi) albint=0. norm=0. DO i = 1,nlambda gamm = 4.0 * pi * ni(i) / wavelengths(i) cosalb = 2.0 / (SSA * rhoice) * B_amp * gamm alpha = 16. / 3 * cosalb / (1.0 - g_asy) R = exp(-(alpha**0.5) * 3.0 / 7.0 * (1.0 + 2.0 * cossza)) alb30 = exp(-(alpha**0.5)* 3.0 / 7.0 * (1.0 + 20 * cos30)) IF ((wavelengths(i).GE.lambdamin).AND. & (wavelengths(i).LT.lambdamax)) THEN albint=albint+R*Pas ! rectangle integration -> can be improved with trapezintegration norm=norm+Pas ENDIF END DO albint=max(0.,min(albint/max(norm,1E-10),albmax)) END SUBROUTINE albedo_kokhanovsky