subroutine SISVAT_zSn ! +------------------------------------------------------------------------+ ! | MAR SISVAT_zSn 12-07-2019 MAR | ! | SubRoutine SISVAT_zSn manages the Snow Pack vertical Discretization | ! | | ! +------------------------------------------------------------------------+ ! | | ! | PARAMETERS: knonv: Total Number of columns = | ! | ^^^^^^^^^^ = Total Number of continental grid boxes | ! | X Number of Mosaic Cell per grid box | ! | | ! | INPUT / NLaysv = New Snow Layer Switch | ! | OUTPUT: isnoSV = total Nb of Ice/Snow Layers | ! | ^^^^^^ ispiSV = 0,...,nsno: Uppermost Superimposed Ice Layer | ! | iiceSV = total Nb of Ice Layers | ! | istoSV = 0,...,5 : Snow History (see istdSV data) | ! | | ! | INPUT / TsisSV : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)| ! | OUTPUT: & Snow Temperatures (layers 1,2,...,nsno) [K] | ! | ^^^^^^ ro__SV : Soil/Snow Volumic Mass [kg/m3] | ! | eta_SV : Soil/Snow Water Content [m3/m3] | ! | dzsnSV : Snow Layer Thickness [m] | ! | G1snSV : Dendricity (<0) or Sphericity (>0) of Snow Layer | ! | G2snSV : Sphericity (>0) or Size of Snow Layer | ! | agsnSV : Snow Age [day] | ! | | ! | METHOD: 1) Agregate the thinest Snow Layer | ! | ^^^^^^ if a new Snow Layer has been precipitated (NLaysv = 1) | ! | 2) Divide a too thick Snow Layer except | ! | if the maximum Number of Layer is reached | ! | in this case forces NLay_s = 1 | ! | 3) Agregate the thinest Snow Layer | ! | in order to divide a too thick Snow Layer | ! | at next Time Step when NLay_s = 1 | ! | | ! +------------------------------------------------------------------------+ ! +--Global Variables ! + ================ use VARphy use VAR_SV use VARdSV use VAR0SV use VARxSV use VARySV use surface_data, only: ok_zsn_ii IMPLICIT NONE ! +--Internal Variables ! + ================== integer :: ikl ,isn ,i ! integer :: NLay_s(knonv) ! Split Snow Layer Switch integer :: isagr1(knonv) ! 1st Layer History integer :: isagr2(knonv) ! 2nd Layer History integer :: LstLay ! 0 ====> isnoSV = 1 integer :: isno_n ! Snow Normal.Profile integer :: iice_n ! Ice Normal.Profile integer :: iiceOK ! Ice Switch integer :: icemix ! 0 ====> Agregated Snow+Ice=Snow ! + ! 1 Ice integer :: isn1 (knonv) ! 1st layer to stagger real :: staggr ! stagger Switch real :: WEagre(knonv) ! Snow Water Equivalent Thickness real :: dzthin(knonv) ! Thickness of the thinest layer real :: OKthin ! Swich ON a new thinest layer real :: dz_dif ! difference from ideal discret. real :: thickL ! Thick Layer Indicator real :: OK_ICE ! Swich ON uppermost Ice Layer real :: Agrege(knonv) ! 1. when Agregation constrained real :: dzepsi ! Min Single Snw Layer Thickness real :: dzxmin ! Min Acceptable Layer Thickness real :: dz_min ! Min Layer Thickness real :: dz_max ! Max Layer Thickness real :: dzagr1(knonv) ! 1st Layer Thickness real :: dzagr2(knonv) ! 2nd Layer Thickness real :: T_agr1(knonv) ! 1st Layer Temperature real :: T_agr2(knonv) ! 2nd Layer Temperature real :: roagr1(knonv) ! 1st Layer Density real :: roagr2(knonv) ! 2nd Layer Density real :: etagr1(knonv) ! 1st Layer Water Content real :: etagr2(knonv) ! 2nd Layer Water Content real :: G1agr1(knonv) ! 1st Layer Dendricity/Spher. real :: G1agr2(knonv) ! 2nd Layer Dendricity/Spher. real :: G2agr1(knonv) ! 1st Layer Sphericity/Size real :: G2agr2(knonv) ! 2nd Layer Sphericity/Size real :: agagr1(knonv) ! 1st Layer Age real :: agagr2(knonv) ! 2nd Layer Age ! +--DATA ! + ==== data icemix / 0 / ! 0 ====> Agregated Snow+Ice=Snow data dzepsi / 0.0020/ ! Min single Layer Thickness data dzxmin / 0.0025/ ! Min accept.Layer Thickness ! #EU data dz_min / 0.0050/ ! Min Local Layer Thickness < SMn data dz_min / 0.0040/ ! Min Local Layer Thickness < SMn data dz_max / 0.0300/ ! Min Gener. Layer Thickness ! + CAUTION: dz_max > dz_min*2 is required ! Otherwise re-agregation is ! + ! activated after splitting ! +--Constrains Agregation of too thin Layers ! + ================================================= ! +--Search the thinest non-zero Layer ! + ---------------------------------- DO ikl=1,knonv if(isnoSV(ikl)<=2) dz_min=max(0.0050,dz_min) dzepsi=0.0015 if(ro__SV(ikl,isnoSV(ikl))>920) dzepsi=0.0020 dzthin(ikl) = 0. ! Arbitrary unrealistic END DO ! Layer Thickness !XF DO ikl=1,knonv DO isn=1,isnoSV(ikl)-3 ! no agregation of 3 first snowlayers ! ! XF 04/07/2019 isno_n = isnoSV(ikl)-isn+1 ! Snow Normal.Profile iice_n = iiceSV(ikl)-isn ! Ice Normal.Profile iiceOK = min(1,max(0,iice_n +1)) ! Ice Switch ! #vz dz_ref(isn) = ! ! #vz. dz_min *((1-iiceOK)*isno_n*isno_n ! Theoretical Profile ! #vz. + iiceOK * 2**iice_n) ! ! #vz. /max(1,isnoSV(ikl)) ! dz_dif = max(zero, & ! Actual Profile dz_min & ! *((1-iiceOK)*isno_n*isno_n & ! Theoretical Profile + iiceOK *2. **iice_n) & ! - dzsnSV(ikl, isn) ) ! Actual Profile ! #vz dzwdif(isn) = dz_dif ! OKthin = max(zero, & ! sign(unun, & ! dz_dif-dzthin(ikl))) & ! 1.=> New thinest Lay. * max(0, & ! 1 => .le. isnoSV min(1, & ! 1 => isn is in the isnoSV(ikl)-isn +1 )) & ! Snow Pack * min(unun, & ! ! ! 1st additional Condition to accept OKthin max(zero, & ! combination sign(unun,G1snSV(ikl, isn ) & ! G1 with same *G1snSV(ikl,max(1,isn-1)))) & ! sign => OK ! ! 2nd additional Condition to accept OKthin + max(zero, & ! G1>0 sign(unun,G1snSV(ikl, isn ))) & ! =>OK ! ! 3rd additional Condition to accept OKthin + max(zero, & ! dz too small sign(unun,dzxmin & ! =>OK -dzsnSV(ikl, isn ))))! i_thin(ikl) = (1. - OKthin) * i_thin(ikl) & ! Update thinest Lay. + OKthin * isn ! Index dzthin(ikl) = (1. - OKthin) * dzthin(ikl) & ! + OKthin * dz_dif ! END DO END DO DO ikl=1,knonv DO isn=1,isnoSV(ikl) OKthin = max(zero, & ! sign(unun, & ! dz_min & ! -dzsnSV(ikl,isn))) & ! * max(zero, & ! ON if dz > 0 sign(unun, & ! dzsnSV(ikl,isn)-epsi)) & ! *min(1,max(0, & ! Multiple Snow Lay. min (1, & ! Switch = 1 isnoSV(ikl) & ! if isno > iice + 1 -iiceSV(ikl)-1)) & ! ! + ! +int(max(zero, & ! sign(unun, & ! dzepsi & ! Minimum accepted for -dzsnSV(ikl,isn)))) & ! 1 Snow Layer over Ice *int(max(zero, & ! ON if dz > 0 sign(unun, & ! dzsnSV(ikl,isn)-epsi))) & ! *(1 -min (abs(isnoSV(ikl) & ! Switch = 1 -iiceSV(ikl)-1),1)) & ! if isno = iice + 1 ! + ! +max(0, & ! Ice min (1, & ! Switch iiceSV(ikl)+1-isn))) & ! *min(unun, & ! max(zero, & ! combination sign(unun,G1snSV(ikl, isn ) & ! G1>0 + G1<0 *G1snSV(ikl,max(1,isn-1)))) & ! NO + max(zero, & ! sign(unun,G1snSV(ikl, isn ))) & ! + max(zero, & ! sign(unun,dzxmin & ! -dzsnSV(ikl, isn ))))! i_thin(ikl) = (1. - OKthin) * i_thin(ikl) & ! Update thinest Lay. + OKthin * isn ! Index END DO END DO ! + *************** call SISVAT_zCr ! + *************** ! +--Assign the 2 Layers to agregate ! + ------------------------------- DO ikl=1,knonv isn = i_thin(ikl) if(LIndsv(ikl)>0) isn=min(nsno-1,isn) ! cXF isagr1(ikl) = istoSV(ikl,isn) isagr2(ikl) = istoSV(ikl,isn+LIndsv(ikl)) dzagr1(ikl) = dzsnSV(ikl,isn) dzagr2(ikl) = dzsnSV(ikl,isn+LIndsv(ikl)) T_agr1(ikl) = TsisSV(ikl,isn) T_agr2(ikl) = TsisSV(ikl,isn+LIndsv(ikl)) roagr1(ikl) = ro__SV(ikl,isn) roagr2(ikl) = ro__SV(ikl,isn+LIndsv(ikl)) etagr1(ikl) = eta_SV(ikl,isn) etagr2(ikl) = eta_SV(ikl,isn+LIndsv(ikl)) G1agr1(ikl) = G1snSV(ikl,isn) G1agr2(ikl) = G1snSV(ikl,isn+LIndsv(ikl)) G2agr1(ikl) = G2snSV(ikl,isn) G2agr2(ikl) = G2snSV(ikl,isn+LIndsv(ikl)) agagr1(ikl) = agsnSV(ikl,isn) agagr2(ikl) = agsnSV(ikl,isn+LIndsv(ikl)) LstLay = min(1,max( 0,isnoSV(ikl) -1)) ! 0 if single Layer isnoSV(ikl) = isnoSV(ikl) & ! decrement isnoSV -(1-LstLay)* max(zero, & ! if downmost Layer sign(unun,eps_21 & ! < 1.e-21 m -dzsnSV(ikl,1))) ! isnoSV(ikl) = max( 0, isnoSV(ikl) ) ! Agrege(ikl) = max(zero, & ! sign(unun,dz_min & ! No Agregation -dzagr1(ikl) )) & ! if too thick Layer *LstLay & ! if a single Layer * min( max(0 ,isnoSV(ikl)+1 & ! if Agregation -i_thin(ikl) & ! with a Layer -LIndsv(ikl) ),1) ! above the Pack WEagre(ikl) = 0. END DO DO ikl=1,knonv DO isn=1,isnoSV(ikl) WEagre(ikl) = WEagre(ikl) + ro__SV(ikl,isn)*dzsnSV(ikl,isn) & *min(1,max(0,i_thin(ikl)+1-isn)) ENDDO ENDDO ! +--Agregates ! + --------- ! + *************** call SISVAT_zAg & (isagr1,isagr2,WEagre & ,dzagr1,dzagr2,T_agr1,T_agr2 & ,roagr1,roagr2,etagr1,etagr2 & ,G1agr1,G1agr2,G2agr1,G2agr2 & ,agagr1,agagr2,Agrege & ) ! + *************** ! +--Rearranges the Layers ! + --------------------- ! +--New (agregated) Snow layer ! + ^^^^^^^^^^^^^^^^^^^^^^^^^^ DO ikl=1,knonv isn = i_thin(ikl) isn = min(isn,isn+LIndsv(ikl)) isnoSV(ikl) = max(0.,isnoSV(ikl) -Agrege(ikl)) iiceSV(ikl) = iiceSV(ikl) & -max(0,sign(1,iiceSV(ikl) -isn +icemix)) & *Agrege(ikl) & *max(0,sign(1,iiceSV(ikl) -1 )) istoSV(ikl,isn) = (1.-Agrege(ikl))*istoSV(ikl,isn) & + Agrege(ikl) *isagr1(ikl) dzsnSV(ikl,isn) = (1.-Agrege(ikl))*dzsnSV(ikl,isn) & + Agrege(ikl) *dzagr1(ikl) TsisSV(ikl,isn) = (1.-Agrege(ikl))*TsisSV(ikl,isn) & + Agrege(ikl) *T_agr1(ikl) ro__SV(ikl,isn) = (1.-Agrege(ikl))*ro__SV(ikl,isn) & + Agrege(ikl) *roagr1(ikl) eta_SV(ikl,isn) = (1.-Agrege(ikl))*eta_SV(ikl,isn) & + Agrege(ikl) *etagr1(ikl) G1snSV(ikl,isn) = (1.-Agrege(ikl))*G1snSV(ikl,isn) & + Agrege(ikl) *G1agr1(ikl) G2snSV(ikl,isn) = (1.-Agrege(ikl))*G2snSV(ikl,isn) & + Agrege(ikl) *G2agr1(ikl) agsnSV(ikl,isn) = (1.-Agrege(ikl))*agsnSV(ikl,isn) & + Agrege(ikl) *agagr1(ikl) END DO ! +--Above ! + ^^^^^ DO ikl=1,knonv isn1(ikl)=max(i_thin(ikl),i_thin(ikl)+LIndsv(ikl)) END DO DO i= 1,nsno-1 DO ikl=1,knonv staggr = min(1,max(0,i +1 -isn1(ikl) )) istoSV(ikl,i) = (1.-staggr )*istoSV(ikl,i ) & + staggr*((1.-Agrege(ikl))*istoSV(ikl,i ) & + Agrege(ikl) *istoSV(ikl,i+1)) dzsnSV(ikl,i) = (1.-staggr )*dzsnSV(ikl,i ) & + staggr*((1.-Agrege(ikl))*dzsnSV(ikl,i ) & + Agrege(ikl) *dzsnSV(ikl,i+1)) TsisSV(ikl,i) = (1.-staggr )*TsisSV(ikl,i ) & + staggr*((1.-Agrege(ikl))*TsisSV(ikl,i ) & + Agrege(ikl) *TsisSV(ikl,i+1)) ro__SV(ikl,i) = (1.-staggr )*ro__SV(ikl,i ) & + staggr*((1.-Agrege(ikl))*ro__SV(ikl,i ) & + Agrege(ikl) *ro__SV(ikl,i+1)) eta_SV(ikl,i) = (1.-staggr )*eta_SV(ikl,i ) & + staggr*((1.-Agrege(ikl))*eta_SV(ikl,i ) & + Agrege(ikl) *eta_SV(ikl,i+1)) G1snSV(ikl,i) = (1.-staggr )*G1snSV(ikl,i ) & + staggr*((1.-Agrege(ikl))*G1snSV(ikl,i ) & + Agrege(ikl) *G1snSV(ikl,i+1)) G2snSV(ikl,i) = (1.-staggr )*G2snSV(ikl,i ) & + staggr*((1.-Agrege(ikl))*G2snSV(ikl,i ) & + Agrege(ikl) *G2snSV(ikl,i+1)) agsnSV(ikl,i) = (1.-staggr )*agsnSV(ikl,i ) & + staggr*((1.-Agrege(ikl))*agsnSV(ikl,i ) & + Agrege(ikl) *agsnSV(ikl,i+1)) END DO END DO DO ikl=1,knonv isn = min(isnoSV(ikl) +1,nsno) istoSV(ikl,isn) = (1.-Agrege(ikl))*istoSV(ikl,isn) dzsnSV(ikl,isn) = (1.-Agrege(ikl))*dzsnSV(ikl,isn) TsisSV(ikl,isn) = (1.-Agrege(ikl))*TsisSV(ikl,isn) ro__SV(ikl,isn) = (1.-Agrege(ikl))*ro__SV(ikl,isn) eta_SV(ikl,isn) = (1.-Agrege(ikl))*eta_SV(ikl,isn) G1snSV(ikl,isn) = (1.-Agrege(ikl))*G1snSV(ikl,isn) G2snSV(ikl,isn) = (1.-Agrege(ikl))*G2snSV(ikl,isn) agsnSV(ikl,isn) = (1.-Agrege(ikl))*agsnSV(ikl,isn) END DO ! +--Constrains Splitting of too thick Layers ! + ================================================= ! +--Search the thickest non-zero Layer ! + ---------------------------------- DO ikl=1,knonv dzthin(ikl) = 0. ! Arbitrary unrealistic END DO DO ikl=1,knonv DO isn=1,isnoSV(ikl) isno_n = isnoSV(ikl)-isn+1 ! Snow Normal.Profile iice_n = iiceSV(ikl)-isn ! Ice Normal.Profile iiceOK = min(1,max(0,iice_n +1)) ! Ice Switch dz_dif =( dzsnSV(ikl,isn) & ! Actual Profile - dz_max *((1-iiceOK)*isno_n*isno_n & ! Theoretical Profile + iiceOK *2. **iice_n) ) & ! /max(dzsnSV(ikl,isn),epsi) ! OKthin = max(zero, & ! sign(unun, & ! dz_dif-dzthin(ikl))) & ! 1.=>New thickest Lay. * max(0, & ! 1 =>.le. isnoSV min(1, & ! isnoSV(ikl)-isn +1 )) ! i_thin(ikl) = (1. - OKthin) * i_thin(ikl) & ! Update thickest Lay. + OKthin * isn ! Index dzthin(ikl) = (1. - OKthin) * dzthin(ikl) & ! + OKthin * dz_dif ! END DO isn=max(1,isnoSV(ikl)-3) if(dzsnSV(ikl,isn)>0.30) then ! surface layer > 30cm i_thin(ikl) = isn ! XF 04/07/2019 dzthin(ikl) = dzsnSV(ikl,isn) endif isn=max(1,isnoSV(ikl)-2) if(dzsnSV(ikl,isn)>0.10) then ! surface layer > 10cm i_thin(ikl) = isn ! XF 04/07/2019 dzthin(ikl) = dzsnSV(ikl,isn) endif isn=max(1,isnoSV(ikl)-1) if(dzsnSV(ikl,isn)>0.05) then ! surface layer > 5cm i_thin(ikl) = isn ! XF 04/07/2019 dzthin(ikl) = dzsnSV(ikl,isn) endif isn=max(1,isnoSV(ikl)) if(dzsnSV(ikl,isn)>0.02) then ! surface layer > 2cm i_thin(ikl) = isn ! XF 04/07/2019 dzthin(ikl) = dzsnSV(ikl,isn) endif END DO DO ikl=1,knonv ThickL = max(zero, & ! 1. => a too thick sign(unun,dzthin(ikl) & ! Layer exists -epsi )) & ! * max(0,1-max(0 , isnoSV(ikl) & ! No spliting allowed -nsno+1 )) ! if isno > nsno - 1 Agrege(ikl) = ThickL & ! 1. => effective split * max(0,1-max(0 , NLaysv(ikl) & ! +isnoSV(ikl) & ! -nsno+1 )) ! NLay_s(ikl) = ThickL & ! Agregation * max(0,1-max(0 , NLaysv(ikl) & ! to allow Splitting +isnoSV(ikl) & ! at next Time Step -nsno )) & ! -Agrege(ikl) ! NLay_s(ikl) = max(0 , NLay_s(ikl)) ! Agregation effective END DO ! +--Rearranges the Layers ! + --------------------- DO isn=nsno,2,-1 DO ikl=1,knonv IF (Agrege(ikl).gt.0..AND.i_thin(ikl).lt.isnoSV(ikl)) THEN staggr = min(1,max(0,isn-i_thin(ikl) -1)) & * min(1,max(0, isnoSV(ikl)-isn+2)) istoSV(ikl,isn) = staggr * istoSV(ikl ,isn-1) & + (1. - staggr) * istoSV(ikl ,isn ) dzsnSV(ikl,isn) = staggr * dzsnSV(ikl ,isn-1) & + (1. - staggr) * dzsnSV(ikl ,isn ) TsisSV(ikl,isn) = staggr * TsisSV(ikl ,isn-1) & + (1. - staggr) * TsisSV(ikl ,isn ) ro__SV(ikl,isn) = staggr * ro__SV(ikl ,isn-1) & + (1. - staggr) * ro__SV(ikl ,isn ) eta_SV(ikl,isn) = staggr * eta_SV(ikl ,isn-1) & + (1. - staggr) * eta_SV(ikl ,isn ) G1snSV(ikl,isn) = staggr * G1snSV(ikl ,isn-1) & + (1. - staggr) * G1snSV(ikl ,isn ) G2snSV(ikl,isn) = staggr * G2snSV(ikl ,isn-1) & + (1. - staggr) * G2snSV(ikl ,isn ) agsnSV(ikl,isn) = staggr * agsnSV(ikl ,isn-1) & + (1. - staggr) * agsnSV(ikl ,isn ) END IF END DO END DO DO ikl=1,knonv isn = i_thin(ikl) dzsnSV(ikl,isn) = 0.5*Agrege(ikl) *dzsnSV(ikl,isn) & + (1.-Agrege(ikl))*dzsnSV(ikl,isn) isn = min(i_thin(ikl) +1,nsno) istoSV(ikl,isn) = Agrege(ikl) *istoSV(ikl,isn-1) & + (1.-Agrege(ikl))*istoSV(ikl,isn) dzsnSV(ikl,isn) = Agrege(ikl) *dzsnSV(ikl,isn-1) & + (1.-Agrege(ikl))*dzsnSV(ikl,isn) TsisSV(ikl,isn) = Agrege(ikl) *TsisSV(ikl,isn-1) & + (1.-Agrege(ikl))*TsisSV(ikl,isn) ro__SV(ikl,isn) = Agrege(ikl) *ro__SV(ikl,isn-1) & + (1.-Agrege(ikl))*ro__SV(ikl,isn) eta_SV(ikl,isn) = Agrege(ikl) *eta_SV(ikl,isn-1) & + (1.-Agrege(ikl))*eta_SV(ikl,isn) G1snSV(ikl,isn) = Agrege(ikl) *G1snSV(ikl,isn-1) & + (1.-Agrege(ikl))*G1snSV(ikl,isn) G2snSV(ikl,isn) = Agrege(ikl) *G2snSV(ikl,isn-1) & + (1.-Agrege(ikl))*G2snSV(ikl,isn) agsnSV(ikl,isn) = Agrege(ikl) *agsnSV(ikl,isn-1) & + (1.-Agrege(ikl))*agsnSV(ikl,isn) isnoSV(ikl) = min(Agrege(ikl) +isnoSV(ikl),real(nsno)) iiceSV(ikl) = iiceSV(ikl) & + Agrege(ikl) *max(0,sign(1,iiceSV(ikl) & -isn +icemix)) & *max(0,sign(1,iiceSV(ikl) & -1 )) END DO ! +--Constrains Agregation in case of too much Layers ! + ================================================= ! +--Search the thinest non-zero Layer ! + ----------------------------------- DO ikl=1,knonv dzthin(ikl) = 0. ! Arbitrary unrealistic END DO ! Layer Thickness DO ikl=1,knonv DO isn=1,isnoSV(ikl)-3 ! no agregation of 3 first snowlayers ! ! XF 04/07/2019 isno_n = isnoSV(ikl)-isn+1 ! Snow Normal.Profile iice_n = iiceSV(ikl)-isn ! Ice Normal.Profile iiceOK = min(1,max(0,iice_n +1)) ! Ice Switch ! #vz dz_ref(isn) = ! ! #vz. dz_min *((1-iiceOK)*isno_n*isno_n ! Theoretical Profile ! #vz. + iiceOK * 2**iice_n) ! ! #vz. /max(1,isnoSV(ikl)) ! dz_dif = dz_min & ! Actual Profile - dzsnSV(ikl ,isn) & ! /max(epsi,((1-iiceOK)*isno_n*isno_n & ! Theoretical Profile + iiceOK *2. **iice_n)) ! ! #vz dzwdif(isn) = dz_dif ! OKthin = max(zero, & ! sign(unun, & ! dz_dif - dzthin(ikl))) & ! 1.=> New thinest Lay. * max(0, & ! 1 => .le. isnoSV min(1, & ! isnoSV(ikl)-isn +1 )) ! i_thin(ikl) = (1. - OKthin) * i_thin(ikl) & ! Update thinest Lay. + OKthin * isn ! Index dzthin(ikl) = (1. - OKthin) * dzthin(ikl) & ! + OKthin * dz_dif ! END DO END DO ! +--Index of the contiguous Layer to agregate ! + ----------------------------------------- ! + *************** call SISVAT_zCr ! + *************** ! +--Assign the 2 Layers to agregate ! + ------------------------------- DO ikl=1,knonv isn = i_thin(ikl) if(LIndsv(ikl)>0) isn=min(isn, nsno-1) !cXF isagr1(ikl) = istoSV(ikl,isn) isagr2(ikl) = istoSV(ikl,isn+LIndsv(ikl)) dzagr1(ikl) = dzsnSV(ikl,isn) dzagr2(ikl) = dzsnSV(ikl,isn+LIndsv(ikl)) T_agr1(ikl) = TsisSV(ikl,isn) T_agr2(ikl) = TsisSV(ikl,isn+LIndsv(ikl)) roagr1(ikl) = ro__SV(ikl,isn) roagr2(ikl) = ro__SV(ikl,isn+LIndsv(ikl)) etagr1(ikl) = eta_SV(ikl,isn) etagr2(ikl) = eta_SV(ikl,isn+LIndsv(ikl)) G1agr1(ikl) = G1snSV(ikl,isn) G1agr2(ikl) = G1snSV(ikl,isn+LIndsv(ikl)) G2agr1(ikl) = G2snSV(ikl,isn) G2agr2(ikl) = G2snSV(ikl,isn+LIndsv(ikl)) agagr1(ikl) = agsnSV(ikl,isn) agagr2(ikl) = agsnSV(ikl,isn+LIndsv(ikl)) LstLay = min(1,max( 0, isnoSV(ikl)-1 )) Agrege(ikl) = min(1, & max(0, & NLaysv(ikl) +isnoSV(ikl)-nsno & +NLay_s(ikl) ) & *LstLay ) ! + minimum uppermost layer thickness to guarantee a correct reproduction of the snow ! + atmosphere coupling if(dzsnSV(ikl,max(1,isnoSV(ikl)-0))>0.02 .or. & ! surface layers> 2-5-10 dzsnSV(ikl,max(1,isnoSV(ikl)-1))>0.05 .or. & ! XF 04/07/2019 dzsnSV(ikl,max(1,isnoSV(ikl)-2))>0.10 .or. & dzsnSV(ikl,max(1,isnoSV(ikl)-3))>0.30 )then Agrege(ikl) = min(1, & max(0, & NLaysv(ikl) +isnoSV(ikl)+1-nsno & ! nsno-1 layers ma +NLay_s(ikl) ) & *LstLay ) endif isnoSV(ikl) = isnoSV(ikl) & -(1-LstLay)*max(zero, & sign(unun, eps_21 & -dzsnSV(ikl,1) )) isnoSV(ikl) =max( 0, isnoSV(ikl) ) WEagre(ikl) = 0. END DO DO isn=1,nsno DO ikl=1,knonv WEagre(ikl) = WEagre(ikl) + ro__SV(ikl,isn)*dzsnSV(ikl,isn) & *min(1,max(0,i_thin(ikl)+1-isn)) ENDDO ENDDO ! +--Agregates ! + --------- ! + *************** call SISVAT_zAg & (isagr1,isagr2,WEagre & ,dzagr1,dzagr2,T_agr1,T_agr2 & ,roagr1,roagr2,etagr1,etagr2 & ,G1agr1,G1agr2,G2agr1,G2agr2 & ,agagr1,agagr2,Agrege & ) ! + *************** ! +--Rearranges the Layers ! + --------------------- ! +--New (agregated) Snow layer ! + ^^^^^^^^^^^^^^^^^^^^^^^^^^ DO ikl=1,knonv isn = i_thin(ikl) isn = min(isn,isn+LIndsv(ikl)) isnoSV(ikl) = max(0.,isnoSV(ikl) -Agrege(ikl)) iiceSV(ikl) = iiceSV(ikl) & -max(0,sign(1,iiceSV(ikl) -isn +icemix)) & *Agrege(ikl) & *max(0,sign(1,iiceSV(ikl) -1 )) istoSV(ikl,isn) = (1.-Agrege(ikl))*istoSV(ikl,isn) & + Agrege(ikl) *isagr1(ikl) dzsnSV(ikl,isn) = (1.-Agrege(ikl))*dzsnSV(ikl,isn) & + Agrege(ikl) *dzagr1(ikl) TsisSV(ikl,isn) = (1.-Agrege(ikl))*TsisSV(ikl,isn) & + Agrege(ikl) *T_agr1(ikl) ro__SV(ikl,isn) = (1.-Agrege(ikl))*ro__SV(ikl,isn) & + Agrege(ikl) *roagr1(ikl) eta_SV(ikl,isn) = (1.-Agrege(ikl))*eta_SV(ikl,isn) & + Agrege(ikl) *etagr1(ikl) G1snSV(ikl,isn) = (1.-Agrege(ikl))*G1snSV(ikl,isn) & + Agrege(ikl) *G1agr1(ikl) G2snSV(ikl,isn) = (1.-Agrege(ikl))*G2snSV(ikl,isn) & + Agrege(ikl) *G2agr1(ikl) agsnSV(ikl,isn) = (1.-Agrege(ikl))*agsnSV(ikl,isn) & + Agrege(ikl) *agagr1(ikl) END DO ! +--Above ! + ^^^^^ DO ikl=1,knonv isn1(ikl)=max(i_thin(ikl),i_thin(ikl)+LIndsv(ikl)) END DO DO i= 1,nsno-1 DO ikl=1,knonv staggr = min(1,max(0,i +1 -isn1(ikl) )) istoSV(ikl,i) = (1.-staggr )*istoSV(ikl,i ) & + staggr*((1.-Agrege(ikl))*istoSV(ikl,i ) & + Agrege(ikl) *istoSV(ikl,i+1)) dzsnSV(ikl,i) = (1.-staggr )*dzsnSV(ikl,i ) & + staggr*((1.-Agrege(ikl))*dzsnSV(ikl,i ) & + Agrege(ikl) *dzsnSV(ikl,i+1)) TsisSV(ikl,i) = (1.-staggr )*TsisSV(ikl,i ) & + staggr*((1.-Agrege(ikl))*TsisSV(ikl,i ) & + Agrege(ikl) *TsisSV(ikl,i+1)) ro__SV(ikl,i) = (1.-staggr )*ro__SV(ikl,i ) & + staggr*((1.-Agrege(ikl))*ro__SV(ikl,i ) & + Agrege(ikl) *ro__SV(ikl,i+1)) eta_SV(ikl,i) = (1.-staggr )*eta_SV(ikl,i ) & + staggr*((1.-Agrege(ikl))*eta_SV(ikl,i ) & + Agrege(ikl) *eta_SV(ikl,i+1)) G1snSV(ikl,i) = (1.-staggr )*G1snSV(ikl,i ) & + staggr*((1.-Agrege(ikl))*G1snSV(ikl,i ) & + Agrege(ikl) *G1snSV(ikl,i+1)) G2snSV(ikl,i) = (1.-staggr )*G2snSV(ikl,i ) & + staggr*((1.-Agrege(ikl))*G2snSV(ikl,i ) & + Agrege(ikl) *G2snSV(ikl,i+1)) agsnSV(ikl,i) = (1.-staggr )*agsnSV(ikl,i ) & + staggr*((1.-Agrege(ikl))*agsnSV(ikl,i ) & + Agrege(ikl) *agsnSV(ikl,i+1)) END DO END DO DO ikl=1,knonv isn = min(isnoSV(ikl) +1,nsno) istoSV(ikl,isn) = (1.-Agrege(ikl))*istoSV(ikl,isn) dzsnSV(ikl,isn) = (1.-Agrege(ikl))*dzsnSV(ikl,isn) TsisSV(ikl,isn) = (1.-Agrege(ikl))*TsisSV(ikl,isn) ro__SV(ikl,isn) = (1.-Agrege(ikl))*ro__SV(ikl,isn) eta_SV(ikl,isn) = (1.-Agrege(ikl))*eta_SV(ikl,isn) G1snSV(ikl,isn) = (1.-Agrege(ikl))*G1snSV(ikl,isn) G2snSV(ikl,isn) = (1.-Agrege(ikl))*G2snSV(ikl,isn) agsnSV(ikl,isn) = (1.-Agrege(ikl))*agsnSV(ikl,isn) END DO ! +--Search new Ice/Snow Interface (option II in MAR) ! + =============================================== IF (ok_zsn_ii) THEN DO ikl=1,knonv iiceSV(ikl) = 0 END DO DO ikl=1,knonv DO isn=1,isnoSV(ikl) OK_ICE = max(zero,sign(unun,ro__SV(ikl,isn)-ro_ice+20.)) & * max(zero,sign(unun,dzsnSV(ikl,isn)-epsi)) iiceSV(ikl) = (1.-OK_ICE) *iiceSV(ikl) & + OK_ICE *isn END DO END DO END IF return end subroutine sisvat_zsn