subroutine SnOptP(jjtime) C +------------------------------------------------------------------------+ C | MAR/SISVAT SnOptP 12-08-2019 MAR | C | SubRoutine SnOptP computes the Snow Pack optical Properties | C +------------------------------------------------------------------------+ C | | C | PARAMETERS: knonv: Total Number of columns = | C | ^^^^^^^^^^ = Total Number of continental Grid Boxes | C | X Number of Mosaic Cell per Grid Box | C | | C | INPUT: isnoSV = total Nb of Ice/Snow Layers | C | ^^^^^ ispiSV = 0,...,nsno: Uppermost Superimposed Ice Layer | C | | C | | C | INPUT: G1snSV : Dendricity (<0) or Sphericity (>0) of Snow Layer | C | ^^^^^ G2snSV : Sphericity (>0) or Size of Snow Layer | C | agsnSV : Snow Age [day] | C | ro__SV : Snow/Soil Volumic Mass [kg/m3] | C | eta_SV : Water Content [m3/m3] | C | rusnSV : Surficial Water Thickness [kg/m2] .OR. [mm] | C | SWS_SV : Surficial Water Status | C | dzsnSV : Snow Layer Thickness [m] | C | | C | albssv : Soil Albedo [-] | C | zzsnsv : Snow Pack Thickness [m] | C | | C | OUTPUT: albisv : Snow/Ice/Water/Soil Integrated Albedo [-] | C | ^^^^^^ sEX_sv : Verticaly Integrated Extinction Coefficient | C | DOPsnSV : Snow Optical diameter [m] | C | | C | Internal Variables: | C | ^^^^^^^^^^^^^^^^^^ | C | SnOpSV : Snow Grain optical Size [m] | C | EX1_sv : Integrated Snow Extinction (0.3--0.8micr.m) | C | EX2_sv : Integrated Snow Extinction (0.8--1.5micr.m) | C | EX3_sv : Integrated Snow Extinction (1.5--2.8micr.m) | C | | C | METHODE: Calcul de la taille optique des grains ? partir de | C | ^^^^^^^ -leur type decrit par les deux variables descriptives | C | continues sur la plage -99/+99 passees en appel. | C | -la taille optique (1/10mm) des etoiles, | C | des grains fins et | C | des jeunes faces planes | C | | C | METHOD: Computation of the optical diameter of the grains | C | ^^^^^^ described with the CROCUS formalism G1snSV / G2snSV | C | | C | REFERENCE: Brun et al. 1989, J. Glaciol 35 pp. 333--342 | C | ^^^^^^^^^ Brun et al. 1992, J. Glaciol 38 pp. 13-- 22 | C | Eric Martin Sept.1996 | C | | C | | C +------------------------------------------------------------------------+ C +--Global Variables C + ================ use VARphy use VAR_SV use VARdSV use VARxSV use VARySV use VARtSV USE surface_data, only: iflag_albcalc,correc_alb IMPLICIT NONE C + -- INPUT integer jjtime C +--Internal Variables C + ================== 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 c #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 C +--Local DATA C + ============ C +--For the computation of the solar irradiance extinction in snow C + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ data sbeta1/0.0192/,sbeta2/0.4000/,sbeta3/0.1098/ data sbeta4/1.0000/ data sbeta5/2.00e1/ C +--Snow Age Maximum (Taiga, e.g. Col de Porte) C + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ data AgeMax /60.0/ C +... AgeMax: Snow Age Maximum [day] data AlbMin /0.94/ C +... AlbMin: Albedo Minimum / visible (0.3--0.8 micrometers) data HSnoSV /0.01/ C +... HSnoSV: Snow Thickness through witch C + Albedo is interpolated to Ice Albedo data HIceSV /0.10/ C +... HIceSV: Snow/Ice Thickness through witch C + Albedo is interpolated to Soil Albedo data doptmx /2.3e-3/ C +... doptmx: Maximum optical Diameter (pi * R**2) [m] C + data czeMAX /0.173648178/ ! 80.deg (Segal et al., 1991 JAS) data bsegal /4.00 / ! data albmax /0.99 / ! Albedo max C +-- 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/ C +--Snow Grain optical Size C + ======================= DO ikl=1,knonv DO isn=1,max(1,isnoSV(ikl)) G2snSV(ikl,isn) = max(epsi,G2snSV(ikl,isn)) C +... Avoid non physical Values SignG1 = sign(unun,G1snSV(ikl,isn)) Sph_OK = max(zero,SignG1) SnOpSV(ikl,isn) = 1.e-4 * C +... SI: (from 1/10 mm to m) C +--Contribution of Non Dendritic Snow C + ---------------------------------- . ( Sph_OK *( G2snSV(ikl,isn)*G1snSV(ikl,isn)/G1_dSV . +max(demi*G2snSV(ikl,isn),DFcdSV) . *(unun-G1snSV(ikl,isn) /G1_dSV)) C +--Contribution of Dendritic Snow C + ---------------------------------- . +(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)) C + --For outputs (Etienne) C + ------------------------ DOPsnSV(ikl,isn)=SnOpSV(ikl,isn) END DO END DO C +--Snow/Ice Albedo C + =============== C +--Uppermost effective Snow Layer C + ------------------------------ 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)) cCA +--Correction of snow albedo for Antarctica/Greenland cCA -------------------------------------------------- albCor = correc_alb c #GL albCor = 1.01 c #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)) c #cz dalbeS = ((bsegal+unun)/(unun+2.0*bsegal*csegal) c #cz. - unun )*0.32 c #cz. / bsegal c #cz dalbeS = max(dalbeS,zero) c #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 c #GL albSn1 = albSn1*max(0.9,(1.-1.5*eta_SV(ikl,isn))) c #GL albSn2 = albSn2*max(0.9,(1.-1.5*eta_SV(ikl,isn))) c #GL albSn3 = albSn3*max(0.9,(1.-1.5*eta_SV(ikl,isn))) ENDIF albSno = So1dSV*albSn1 . + So2dSV*albSn2 . + So3dSV*albSn3 cXF 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) c + ro < roSdSV | min al > aI3dSV c + roSdSV < ro < rocdSV | aI2dSV < min al < aI3dSV (fct of density) C +--Snow/Ice Pack Thickness C + ----------------------- 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) C + The value of SnownH is set to 1 in case of ice lenses above C + 1m of dry snow (ro<600kg/m3) for using CROCUS albedo c ro_ave = 0. c dz_ave = 0. c SnowOK = 1. c do isn = isnoSV(ikl),1,-1 c ro_ave = ro_ave + ro__SV(ikl,isn) * dzsnSV(ikl,isn) * SnowOK c dz_ave = dz_ave + dzsnSV(ikl,isn) * SnowOK c SnowOK = max(zero,sign(unun,1.-dz_ave)) c enddo c ro_ave = ro_ave / max(dz_ave,epsi) c SnowOK = max(zero,sign(unun,600.-ro_ave)) c SnownH = SnowOK + SnownH * (1. - SnowOK) C +--Integrated Snow/Ice Albedo: Case of Water on Bare Ice C + ----------------------------------------------------- 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) ! c albWIc = max(aI1dSV,min(aI2dSV,albWIc+slopSV(ikl)* c . 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)) c + rocdSV < ro < ro_ice | aI2dSV< al ro_ice | aI1dSV< al can be improved with trapezintegration norm=norm+Pas ENDIF END DO albint=max(0.,min(albint/max(norm,1E-10),albmax)) END