SUBROUTINE SISVAT_qSn & ( & ! #e1. EqSn_0,EqSn_1,EqSn_d ! #m1. ,SIsubl,SImelt,SIrnof ) ! +------------------------------------------------------------------------+ ! | MAR SISVAT_qSn Fri 29-Jul-2011 MAR | ! | SubRoutine SISVAT_qSn updates the Snow Water Content | ! +------------------------------------------------------------------------+ ! | | ! | 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 | ! | ^^^^^ | ! | | ! | INPUT: TaT_SV : SBL Top Temperature [K] | ! | ^^^^^ dt__SV : Time Step [s] | ! | | ! | INPUT / drr_SV : Rain Intensity [kg/m2/s] | ! | OUTPUT: dzsnSV : Snow Layer Thickness [m] | ! | ^^^^^^ eta_SV : Snow Water Content [m3/m3] | ! | ro__SV : Snow/Soil Volumic Mass [kg/m3] | ! | TsisSV : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)| ! | & Snow Temperatures (layers 1,2,...,nsno) [K] | ! | | ! | OUTPUT: SWS_SV : Surficial Water Status | ! | ^^^^^^ | ! | EExcsv : Snow Energy in Excess, initial Forcing [J/m2] | ! | EqSn_d : Snow Energy in Excess, remaining [J/m2] | ! | EqSn_0 : Snow Energy, before Phase Change [J/m2] | ! | EqSn_1 : Snow Energy, after Phase Change [J/m2] | ! | SIsubl : Snow sublimed/deposed Mass [mm w.e.] | ! | SImelt : Snow Melted Mass [mm w.e.] | ! | SIrnof : Surficial Water + Run OFF Change [mm w.e.] | ! | | ! | Internal Variables: | ! | ^^^^^^^^^^^^^^^^^^ | ! | | ! | # OPTIONS: #E0: IO for Verification: Energy Budget | ! | # ^^^^^^^ | ! | # #su: IO for Verification: Slush Diagnostic | ! | | ! | | ! +------------------------------------------------------------------------+ ! +--Global Variables ! + ================ use VARphy use VAR_SV use VARdSV use VAR0SV use VARxSV use VARySV use surface_data, only: is_ok_slush,opt_runoff_ac IMPLICIT NONE ! Energy Budget ! ~~~~~~~~~~~~~~~~~~~~~~ ! #e1 real EqSn_d(knonv) ! Energy in Excess, initial ! #e1 real EqSn_0(knonv) ! Snow Energy, befor Phase Change ! #vm real EqSn01(knonv) ! Snow Energy, after Phase Change ! #vm real EqSn02(knonv) ! Snow Energy, after Phase Change ! ! .AND. Last Melting ! #e1 real EqSn_1(knonv) ! Snow Energy, after Phase Change ! ! .AND. Mass Redistr. ! Snow/Ice (Mass) Budget ! ~~~~~~~~~~~~~~~~~~~~~~ ! #m1 real SIsubl(knonv) ! Snow Deposed Mass ! #m1 real SImelt(knonv) ! Snow Melted Mass ! #m1 real SIrnof(knonv) ! Local Surficial Water + Run OFF ! +--Internal Variables ! + ================== integer :: ikl ,isn ! integer :: nh ! Non erodible Snow: up.lay.Index integer :: LayrOK ! 1 (0) if In(Above) Snow Pack integer :: k_face ! 1 (0) if Crystal(no) faceted integer :: LastOK ! 1 ==> 1! Snow Layer integer :: NOLayr ! 1 Layer Update integer :: noSnow(knonv) ! Nb of Layers Updater integer :: kSlush ! Slush Switch real :: dTSnow ! Temperature [C] real :: EExdum(knonv) ! Energy in Excess when no Snow real :: OKmelt ! 1 (0) if (no) Melting real :: EnMelt ! Energy in excess, for Melting real :: SnHLat ! Energy consumed in Melting real :: AdEnrg,B_Enrg ! Additional Energy from Vapor real :: dzVap0,dzVap1 ! Vaporized Thickness [m] real :: dzMelt(knonv) ! Melted Thickness [m] real :: rosDry ! Snow volumic Mass if no Water in real :: PorVol ! Pore volume real :: PClose ! Pore Hole Close OFF Switch real :: SGDiam ! Snow Grain Diameter real :: SGDmax ! Max. Snow Grain Diameter real :: rWater ! Retained Water [kg/m2] real :: drrNEW ! New available Water [kg/m2] real :: rdzNEW ! Snow Mass [kg/m2] real :: rdzsno ! Snow Mass [kg/m2] real :: EnFrez ! Energy Release in Freezing real :: WaFrez ! Water consumed in Melting real :: RapdOK ! 1. ==> Snow melts rapidly real :: ThinOK ! 1. ==> Snow Layer is thin real :: dzepsi ! Minim. Snow Layer Thickness (!) real :: dz_Min ! Minim. Snow Layer Thickness real :: z_Melt ! Last (thin) Layer Melting real :: rusnew ! Surficial Water Thickness [mm] real :: zWater ! Max Slush Water Thickness [mm] real :: zSlush ! Slush Water Thickness [mm] real :: ro_new ! New Snow/ice Density [kg/m3] real :: zc,zt ! Non erod.Snow Thickness[mm w.e.] real :: rusnSV0(knonv) real :: Tsave ! +--OUTPUT of SISVAT Trace Statistics (see assignation in PHY_SISVAT) ! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ integer :: isnnew,isinew,isnUpD,isnitr ! OUTPUT in SISVAT at specified i,j,k,n (see assignation in PHY_SISVAT) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! #wx integer iSV_v1,jSV_v1,nSV_v1,kSV_v1,lSV_v1 ! #wx common/SISVAT_EV/ iSV_v1,jSV_v1,nSV_v1,kSV_v1,lSV_v1 ! +--Energy and Mass Budget ! + ~~~~~~~~~~~~~~~~~~~~~~ ! #vm real WqSn_0(knonv) ! Snow Water+Forcing Initial ! #vm real WqSn_1(knonv) ! Snow Water+Forcing, Final ! #vm logical emopen ! IO Switch ! #vm common/Se_qSn_L/emopen ! ! #vm integer no_err ! ! #vm common/Se_qSn_I/no_err ! ! #vm real hourer,timeer ! ! #vm common/Se_qSn_R/timeer ! ! +--Slush Diagnostic: IO ! + ~~~~~~~~~~~~~~~~~~~~ ! #vu logical su_opn ! IO Switch ! #vu common/SI_qSn_L/su_opn ! ! +--DATA ! + ==== data dzepsi/0.0001/ ! Minim. Snow Layer Thickness (!) ! #?? data dz_Min/0.005/ ! Minim. Snow Layer Thickness ! ... Warning: Too high for Col de Porte: precludes 1st snow (layer) apparition data dz_Min/2.5e-3/ ! Minim. Snow Layer Thickness data SGDmax/0.003/ ! Maxim. Snow Grain Diameter [m] ! ! (Rowe et al. 1995, JGR p.16268) ! +--Energy Budget (IN) ! + ================== ! #e1 DO ikl=1,knonv ! #e1 EqSn_0(ikl) = 0. ! #e1 END DO ! #e1 DO isn=nsno,1,-1 ! #e1 DO ikl=1,knonv ! #e1 EqSn_0(ikl) = EqSn_0(ikl) + ro__SV(ikl,isn) *dzsnSV(ikl,isn) ! #e1. *(Cn_dSV *(TsisSV(ikl,isn) -TfSnow ) ! #e1. -Lf_H2O *(1. -eta_SV(ikl,isn))) ! #e1 END DO ! #e1 END DO ! +--Water Budget (IN) ! + ================== ! #vm DO ikl=1,knonv ! #vm WqSn_0(ikl) = drr_SV(ikl) * dt__SV ! #vm. +rusnSV(ikl) ! #vm END DO ! #vm DO isn=nsno,1,-1 ! #vm DO ikl=1,knonv ! #vm WqSn_0(ikl) = WqSn_0(ikl) + ro__SV(ikl,isn) *dzsnSV(ikl,isn) ! #vm END DO ! #vm END DO ! +--Snow Melt Budget ! + ================ ! #m1 DO ikl=1,knonv ! #m1 SImelt(ikl) = 0. ! #m1 SIrnof(ikl) = rusnSV(ikl) + RnofSV(ikl) * dt__SV ! #m1 END DO ! +--Initialization ! + ============== DO ikl=1,knonv noSnow(ikl) = 0 ! Nb of Layers Updater ispiSV(ikl) = 0 ! Pore Hole Close OFF Index ! ! (assumed to be the Top of ! ! the surimposed Ice Layer) zn5_SV(ikl) = 0. rusnSV0(ikl) = 0. END DO ! +--Melting/Freezing Energy ! + ======================= ! +...REMARK: Snow liquid Water Temperature assumed = TfSnow ! + ^^^^^^ DO ikl=1,knonv EExdum(ikl) = drr_SV(ikl) * C__Wat *(TaT_SV(ikl)-TfSnow) & * dt__SV EExcsv(ikl) = EExdum(ikl) * min(1,isnoSV(ikl)) ! Snow exists EExdum(ikl) = EExdum(ikl) - EExcsv(ikl) ! ! #e1 EqSn_d(ikl) = EExcsv(ikl) ! END DO ! +--Surficial Water Status ! + ---------------------- DO ikl=1,knonv SWS_SV(ikl) = max(zero,sign(unun,TfSnow & -TsisSV(ikl,isnoSV(ikl)))) END DO DO ikl=1,knonv DO isn=min(nsno,isnoSV(ikl)+1),1,-1 ! EV DO isn=nsno,1,-1 ! +--Energy, store Previous Content ! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ dTSnow = TsisSV(ikl,isn) - TfSnow EExcsv(ikl) = EExcsv(ikl) & + ro__SV(ikl,isn) * Cn_dSV * dTSnow & * dzsnSV(ikl,isn) TsisSV(ikl,isn) = TfSnow ! +--Water, store Previous Content ! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ drr_SV(ikl) = drr_SV(ikl) & + ro__SV(ikl,isn) * eta_SV(ikl,isn) & * dzsnSV(ikl,isn) & / dt__SV ro__SV(ikl,isn) = & ro__SV(ikl,isn) *(1. - eta_SV(ikl,isn)) eta_SV(ikl,isn) = 0. ! +--Melting if EExcsv > 0 ! + ====================== EnMelt = max(zero, EExcsv(ikl) ) ! +--Energy Consumption ! + ^^^^^^^^^^^^^^^^^^ SnHLat = ro__SV(ikl,isn) * Lf_H2O dzMelt(ikl) = EnMelt / max(SnHLat, epsi ) noSnow(ikl) = noSnow(ikl) & + max(zero ,sign(unun,dzMelt(ikl) & ! -dzsnSV(ikl ,isn))) & ! 1 if full Melt *min(1 , max(0 ,1+isnoSV(ikl)-isn)) ! 1 in the Pack dzMelt(ikl) = & min(dzsnSV(ikl, isn),dzMelt(ikl)) dzsnSV(ikl,isn) = & dzsnSV(ikl,isn) -dzMelt(ikl) zn5_SV(ikl) = zn5_SV(ikl) +dzMelt(ikl) EExcsv(ikl) = EExcsv(ikl) -dzMelt(ikl)*SnHLat wem_SV(ikl) = wem_SV(ikl) -dzMelt(ikl)*ro__SV(ikl,isn) ! +--Water Production ! + ^^^^^^^^^^^^^^^^^ drr_SV(ikl) = drr_SV(ikl) & + ro__SV(ikl,isn) * dzMelt(ikl)/dt__SV ! #m1 SImelt(ikl) = SImelt(ikl) ! #m1. + ro__SV(ikl,isn) * dzMelt(ikl) OKmelt =max(zero,sign(unun,drr_SV(ikl)-epsi)) ! +--Snow History ! + ^^^^^^^^^^^^ k_face = min( istoSV(ikl,isn),istdSV(1)) & ! = 1 if *max(0,2-istoSV(ikl,isn) ) ! faceted istoSV(ikl,isn) = & ! (1.-OKmelt) * istoSV(ikl,isn) & ! + OKmelt *((1-k_face) * istdSV(2) & ! + k_face * istdSV(3) ) ! ! +--Freezing if EExcsv < 0 ! + ====================== rdzsno = ro__SV(ikl,isn) * dzsnSV(ikl ,isn) LayrOK = min( 1, max(0 , isnoSV(ikl)-isn+1)) EnFrez = min(zero, EExcsv(ikl)) WaFrez = -( EnFrez * LayrOK / Lf_H2O) drrNEW = max(zero,drr_SV(ikl) - WaFrez / dt__SV) WaFrez = ( drr_SV(ikl) - drrNEW)* dt__SV drr_SV(ikl) = drrNEW EExcsv(ikl) = EExcsv(ikl) + WaFrez * Lf_H2O EnFrez = min(zero,EExcsv(ikl)) * LayrOK rdzNEW = WaFrez + rdzsno ro__SV(ikl,isn) = rdzNEW /max(epsi, dzsnSV(ikl,isn)) TsisSV(ikl,isn) = TfSnow & + EnFrez /(Cn_dSV *max(epsi, rdzNEW) ) EExcsv(ikl) = EExcsv(ikl) - EnFrez wer_SV(ikl) = WaFrez & + wer_SV(ikl) ! +--Snow Water Content ! + ================== ! +--Percolation Velocity ! + ^^^^^^^^^^^^^^^^^^^^ ! #PW SGDiam = 1.6d-4 ! #PW. + 1.1d-13 *(ro__SV(ikl,isn)*ro__SV(ikl,isn) ! #PW. *ro__SV(ikl,isn)*ro__SV(ikl,isn)) ! +--Pore Volume [-] ! + ^^^^^^^^^^^^^^^^^ rosDry =(1. - eta_SV(ikl,isn))* ro__SV(ikl,isn) ! PorVol = 1. - rosDry / ro_Ice ! PorVol = max(PorVol , zero ) ! ! +--Water Retention ! + ^^^^^^^^^^^^^^^^ rWater = ws0dSV * PorVol * ro_Wat * dzsnSV(ikl,isn) drrNEW = max(zero,drr_SV(ikl) - rWater /dt__SV) rWater = ( drr_SV(ikl) - drrNEW)*dt__SV drr_SV(ikl) = drrNEW rdzNEW = rWater & + rosDry * dzsnSV(ikl,isn) eta_SV(ikl,isn) = rWater / max(epsi,rdzNEW) ro__SV(ikl,isn) = rdzNEW / max(epsi,dzsnSV(ikl,isn)) ! +--Pore Hole Close OFF ! + ^^^^^^^^^^^^^^^^^^^ PClose = max(zero, & sign(unun,ro__SV(ikl,isn) & -roCdSV )) ispiSV(ikl) = ispiSV(ikl) *(1.-PClose) & + max(ispiSV(ikl),isn) * Pclose PClose = max(0 , & ! Water under SuPer.Ice min (1 ,ispiSV(ikl) & ! contributes to -isn )) ! Surficial Water !XF if(ro__SV(ikl,isn) >= roCdSV.and.ro__SV(ikl,1)<900) & PClose = min(0.50,PClose * & (1.-(ro_ice-ro__SV(ikl,isn))/(ro_ice-roCdSV))) PClose = max(0.,min(1.,PClose)) if(isn==1) then PClose = 1 ispiSV(ikl)= max(ispiSV(ikl),1) endif if(drr_SV(ikl) >0 .and.TsisSV(ikl,isn)>273.14) then if((ro__SV(ikl,isn)>900.and.ro__SV(ikl,isn)<920).or. & ro__SV(ikl,isn)>950) then dzsnSV(ikl,isn) = dzsnSV(ikl,isn)*ro__SV(ikl,isn)/ro_ice ro__SV(ikl,isn) = ro_ice PClose = 1 endif endif ! if (isn>1.and.isn900 .and. ! . ro__SV(ikl,isn) >roCdSV .and. ! . ro__SV(ikl,isn) <900 .and. ! . TsisSV(ikl,isn) >273.14 .and. ! . TsisSV(ikl,isn+1)<273.15 .and. ! . drr_SV(ikl) >0) then ! TsisSV(ikl,isn)=273.14 ! PClose = 1 ! endif !XF rusnSV(ikl) = rusnSV(ikl) & + drr_SV(ikl) *dt__SV * PClose rusnSV0(ikl)= rusnSV0(ikl) & + drr_SV(ikl) *dt__SV * PClose drr_SV(ikl) = drr_SV(ikl) *(1.-PClose) END DO END DO ! +--Remove Zero-Thickness Layers ! + ============================ 1000 CONTINUE isnitr = 0 DO ikl=1,knonv isnUpD = 0 isinew = 0 !XF DO isn=1,min(nsno-1,isnoSV(ikl)) isnnew =(unun-max(zero ,sign(unun,dzsnSV(ikl,isn)-dzepsi))) & * max(0 , min(1 ,isnoSV(ikl) +1 -isn )) isnUpD = max(isnUpD, isnnew) isnitr = max(isnitr, isnnew) isinew = isn*isnUpD *max(0, 1-isinew) & ! LowerMost 0-Layer +isinew ! Index dzsnSV(ikl,isn) = dzsnSV(ikl,isn+isnnew) ro__SV(ikl,isn) = ro__SV(ikl,isn+isnnew) TsisSV(ikl,isn) = TsisSV(ikl,isn+isnnew) eta_SV(ikl,isn) = eta_SV(ikl,isn+isnnew) G1snSV(ikl,isn) = G1snSV(ikl,isn+isnnew) G2snSV(ikl,isn) = G2snSV(ikl,isn+isnnew) dzsnSV(ikl,isn+isnnew) =(1-isnnew)*dzsnSV(ikl,isn+isnnew) ro__SV(ikl,isn+isnnew) =(1-isnnew)*ro__SV(ikl,isn+isnnew) eta_SV(ikl,isn+isnnew) =(1-isnnew)*eta_SV(ikl,isn+isnnew) G1snSV(ikl,isn+isnnew) =(1-isnnew)*G1snSV(ikl,isn+isnnew) G2snSV(ikl,isn+isnnew) =(1-isnnew)*G2snSV(ikl,isn+isnnew) END DO isnoSV(ikl) = isnoSV(ikl)-isnUpD ! Nb of Snow Layer ispiSV(ikl) = ispiSV(ikl) & ! Nb of SuperI Layer -isnUpD *max(0,min(ispiSV(ikl)-isinew,1)) ! Update if I=0 END DO IF (isnitr>0) GO TO 1000 ! +--New upper Limit of the non erodible Snow (istoSV .GT. 1) ! + ======================================== DO ikl=1,knonv nh = 0 !XF DO isn= isnoSV(ikl),1,-1 nh = nh + isn* min(istoSV(ikl,isn)-1,1)*max(0,1-nh) ENDDO zc = 0. zt = 0. !XF DO isn=1,isnoSV(ikl) zc = zc + dzsnSV(ikl,isn) *ro__SV(ikl,isn) & * max(0,min(1,nh+1-isn)) zt = zt + dzsnSV(ikl,isn) *ro__SV(ikl,isn) END DO zWE_SV(ikl) = zt zWEcSV(ikl) = min(zWEcSV(ikl),zt) zWEcSV(ikl) = max(zWEcSV(ikl),zc) END DO ! +--Energy Budget (OUT) ! + =================== ! #vm DO ikl=1,knonv ! #vm EqSn01(ikl) =-EqSn_0(ikl) ! #vm. -EExcsv(ikl) ! #vm END DO ! #vm DO isn=nsno,1,-1 ! #vm DO ikl=1,knonv ! #vm EqSn01(ikl) = EqSn01(ikl) + ro__SV(ikl,isn) *dzsnSV(ikl,isn) ! #vm. *(Cn_dSV *(TsisSV(ikl,isn) -TfSnow ) ! #vm. -Lf_H2O *(1. -eta_SV(ikl,isn))) ! #vm END DO ! #vm END DO ! +--"Negative Heat" from supercooled rain ! + ------------------------------------ DO ikl=1,knonv EExcsv(ikl) = EExcsv(ikl) + EExdum(ikl) ! +--Surficial Water Run OFF ! + ----------------------- rusnew = rusnSV(ikl) * SWf_SV(ikl) if(isnoSV(ikl)<=1 .OR. opt_runoff_ac) rusnew = 0. ! !if(ivgtSV(ikl)>=1) rusnew = 0. ! #EU rusnew = 0. ! #AC rusnew = 0. RnofSV(ikl) = RnofSV(ikl) & +(rusnSV(ikl) - rusnew ) / dt__SV RuofSV(ikl,1) = RuofSV(ikl,1) & +(rusnSV(ikl) - rusnew ) / dt__SV RuofSV(ikl,4) = RuofSV(ikl,4) & +(rusnSV0(ikl) ) / dt__SV rusnSV(ikl) = rusnew END DO ! +--Percolation down the Continental Ice Pack ! + ----------------------------------------- DO ikl=1,knonv drr_SV(ikl) = drr_SV(ikl) + rusnSV(ikl) & * (1-min(1,ispiSV(ikl)))/ dt__SV rusnSV(ikl) = rusnSV(ikl) & * min(1,ispiSV(ikl)) END DO !XF removal of too thin snowlayers if TT> 275.15 + bug if TT>> 273.15 DO ikl=1,knonv zt=0. DO isn=1,isnoSV(ikl) zt=zt+dzsnSV(ikl,isn) ENDDO if(zt<0.005+(TaT_SV(ikl)-TfSnow)/1000..and. & isnoSV(ikl) >0 .and. & TaT_SV(ikl) >=TfSnow .and. & istoSV(ikl,isnoSV(ikl)) >1 ) then DO isn=1,isnoSV(ikl) drr_SV(ikl) = drr_SV(ikl) & + dzsnSV(ikl,isn)*ro__SV(ikl,isn) /dt__SV dzsnSV(ikl,isn)= 0. ENDDO isnoSV(ikl) = 0 endif ENDDO ! +--Slush Formation (Activated. CAUTION: ADD RunOff Possibility before Activation) ! + --------------- ^^^^^^^ ^^^ IF (is_ok_slush) THEN DO ikl=1,knonv DO isn=1,isnoSV(ikl) kSlush = min(1,max(0,isn+1-ispiSV(ikl))) ! Slush Switch ! +--Available Additional Pore Volume [-] ! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ PorVol = 1. - ro__SV(ikl,isn) & ! [--] *(1. - eta_SV(ikl,isn))/ ro_Ice & ! - eta_SV(ikl,isn) & ! *ro__SV(ikl,isn) / ro_Wat ! PorVol = max(PorVol , zero ) ! zWater = dzsnSV(ikl,isn) * PorVol * 1000. & ! [mm] OR [kg/m2] * (1. -SWS_SV(ikl) & ! 0 <=> freezing *(1 -min(1,iabs(isn-isnoSV(ikl))))) ! 1 <=> isn=isnoSV zSlush = min(rusnSV(ikl) , zWater) ! [mm] OR [kg/m2] ro_new =(dzsnSV(ikl,isn) * ro__SV(ikl,isn) & ! +zSlush ) & ! / max(dzsnSV(ikl,isn) , epsi ) ! if(ro_new