Index: LMDZ6/trunk/libf/phylmd/blowing_snow_ini_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/blowing_snow_ini_mod.F90 (revision 4671) +++ LMDZ6/trunk/libf/phylmd/blowing_snow_ini_mod.F90 (revision 4672) @@ -8,5 +8,5 @@ real coef_eva_bs, expo_eva_bs, fallv_bs, zeta_bs real prt_bs, pbst_bs, qbst_bs - integer :: niter_bs, iflag_saltation_bs + integer :: iflag_saltation_bs !$OMP THREADPRIVATE(prt_level, lunout) @@ -14,5 +14,5 @@ !$OMP THREADPRIVATE(coef_eva_bs, expo_eva_bs, fallv_bs, zeta_bs) !$OMP THREADPRIVATE(pbst_bs, prt_bs, qbst_bs) - !$OMP THREADPRIVATE(niter_bs, iflag_saltation_bs) + !$OMP THREADPRIVATE(iflag_saltation_bs) contains @@ -45,5 +45,5 @@ CALL getin_p('qbst_bs',qbst_bs) - pbst_bs= 0.01 + pbst_bs= 0.0003 CALL getin_p('pbst_bs',pbst_bs) @@ -64,7 +64,4 @@ CALL getin_p('expo_eva_bs',expo_eva_bs) - niter_bs = 5 - CALL getin_p('niter_bs',niter_bs) - iflag_saltation_bs=0 CALL getin_p('iflag_saltation_bs',iflag_saltation_bs) Index: LMDZ6/trunk/libf/phylmd/blowing_snow_sublim_sedim.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/blowing_snow_sublim_sedim.F90 (revision 4671) +++ LMDZ6/trunk/libf/phylmd/blowing_snow_sublim_sedim.F90 (revision 4672) @@ -9,5 +9,5 @@ use blowing_snow_ini_mod, only : coef_eva_bs,RTT,RD,RG,expo_eva_bs, fallv_bs -use blowing_snow_ini_mod, only : RCPD, RLSTT, RLMLT, RLVTT, RVTMP2, niter_bs +use blowing_snow_ini_mod, only : RCPD, RLSTT, RLMLT, RLVTT, RVTMP2 USE lmdz_lscp_tools, only : calc_qsat_ecmwf @@ -125,5 +125,7 @@ ! vapor, temperature, precip fluxes update zq(i) = zq(i) - (sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime - zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime + zq(i) = max(0., zq(i)) + !zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime + !zqbs(i) = max(0., zqbs(i)) ! melting zt(i) = zt(i) & @@ -151,7 +153,9 @@ - ! vapor, temperature, precip fluxes update + ! vapor, temperature, precip fluxes update following sublimation zq(i) = zq(i) - (sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime - zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime + zq(i) = max(0., zq(i)) + !zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime + !zqbs(i) = max(0., zqbs(i)) zt(i) = zt(i) & + (sedimn(i)-sedim(i)) & @@ -173,20 +177,40 @@ ENDDO ! loop on ngrid - ! Now sedimention scheme (alike that in lscp) - DO n = 1, niter_bs - DO i = 1, ngrid - zrho(i,k) = pplay(i,k) / zt(i) / RD - zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG) - velo(i,k) = fallv_bs - dqsedim = dtime/REAL(niter_bs)/zdz(i)*zqbs(i)*velo(i,k) ! dqice/dt=1/rho*d(rho*wice*qice)/dz - precbs = MIN(MAX(dqsedim,0.0),zqbs(i)) - zqbs(i) = MAX(zqbs(i)-1.*precbs , 0.0) - ENDDO !loop on ngrid - ENDDO ! loop on niter_bs - - ! add to non sublimated precip - DO i=1,ngrid - sedim(i) = sedim(i)+max(0.,zqbsi(i)-zqbs(i))*(paprs(i,k)-paprs(i,k+1))/(RG*dtime) + ! Now sedimention scheme + + ! exact resolution of the conservation equation for qbs with the updated flux from the top (after evap) + ! valid only if the fall velocity is constant + + DO i = 1, ngrid + + zrho(i,k) = pplay(i,k) / zt(i) / RD + zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG) + velo(i,k) = fallv_bs + zqbs(i) = zqbsi(i)*exp(-velo(i,k)/zdz(i)*dtime)+sedim(i)/zrho(i,k)/velo(i,k) + zqbs(i) = max(zqbs(i),0.) + ! flux update + sedim(i) = sedim(i) + zrho(i,k)*zdz(i)/dtime*(zqbsi(i)-zqbs(i)) + sedim(i) = max(0.,sedim(i)) + ENDDO + +! old version with bugs +! DO n = 1, niter_bs +! DO i = 1, ngrid +! zrho(i,k) = pplay(i,k) / zt(i) / RD +! zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG) +! velo(i,k) = fallv_bs +! dqsedim = dtime/REAL(niter_bs)/zdz(i)*zqbs(i)*velo(i,k) ! dqice/dt=1/rho*d(rho*wice*qice)/dz +! precbs = MIN(MAX(dqsedim,0.0),zqbs(i)) +! zqbs(i) = MAX(zqbs(i)-1.*precbs , 0.0) +! ENDDO !loop on ngrid +! ENDDO ! loop on niter_bs +! ! add to non sublimated precip +! DO i=1,ngrid +! sedim(i) = sedim(i)+max(0.,zqbsi(i)-zqbs(i))*(paprs(i,k)-paprs(i,k+1))/(RG*dtime) +! ENDDO +! + + ! Outputs: Index: LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90 (revision 4671) +++ LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90 (revision 4672) @@ -140,8 +140,10 @@ REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow REAL, DIMENSION (klon,6) :: alb6 - REAL :: rho0, rhoice, ustart0, hsalt, esalt, rhod + REAL :: rho0, rhoice, ustart0,esalt, rhod REAL :: lambdasalt,fluxsalt, csalt, nunu, aa, bb, cc REAL :: tau_dens, tau_dens0, tau_densmin, rhomax, rhohard - REAL, DIMENSION(klon) :: ws1, rhos, ustart, qsalt + REAL, DIMENSION(klon) :: ws1, rhos, ustart, qsalt, hsalt, snowini + REAL :: maxerosion ! in kg/s + ! End definition !**************************************************************************************** @@ -364,22 +366,24 @@ ustart0 = 0.211 rhoice = 920.0 - rho0 = 200.0 + rho0 = 300.0 rhomax=450.0 - rhohard=400.0 + rhohard=450.0 tau_dens0=86400.0*10. ! 10 days by default, in s - tau_densmin=86400.0 ! 1 days according to in situ obs by C. Amory + tau_densmin=21600.0 ! 1/4 days according to in situ obs by C. Amory during blowing snow + + ! Marshall et al. 1999 (snow densification during rain) ! computation of threshold friction velocity ! which depends on surface snow density - do i = 1, knon + do j = 1, knon ! estimation of snow density ! snow density increases with snow age and - ! increases even faster in case of sedimentation of blowing snow - tau_dens=max(tau_densmin, tau_dens0*exp(-abs(precip_bs(i))/pbst_bs-abs(precip_rain(i))/prt_bs)) - rhos(i)=rho0+(rhohard-rho0)*(1.-exp(-agesno(i)*86400.0/tau_dens)) + ! increases even faster in case of sedimentation of blowing snow or rain + tau_dens=max(tau_densmin, & tau_dens0*exp(-abs(precip_bs(j))/pbst_bs-abs(precip_rain(j))/prt_bs) & + *exp(-max(tsurf(j)-272.0,0.))) + rhos(j)=rho0+(rhohard-rho0)*(1.-exp(-agesno(j)*86400.0/tau_dens)) ! blowing snow flux formula used in MAR - ws1(i)=(u1(i)**2+v1(i)**2)**0.5 - ustar(i)=(cdragm(i)*(u1(i)**2+v1(i)**2))**0.5 - ustart(i)=ustart0*exp(max(rhoice/rho0-rhoice/rhos(i),0.))*exp(max(0.,rhos(i)-rhomax)) + ws1(j)=(u1(j)**2+v1(j)**2)**0.5 + ustar(j)=(cdragm(j)*(u1(j)**2+v1(j)**2))**0.5 + ustart(j)=ustart0*exp(max(rhoice/rho0-rhoice/rhos(j),0.))*exp(max(0.,rhos(j)-rhomax)) ! we have multiplied by exp to prevent erosion when rhos>rhomax (usefull till ! rhohard<450) @@ -397,14 +401,14 @@ cc=2.0 lambdasalt=0.45 - do i =1, knon - rhod=p1lay(i)/RD/temp_air(i) - nunu=max(ustar(i)/ustart(i),1.e-3) - fluxsalt=rhod/RG*(ustar(i)**3)*(1.-nunu**(-2)) * & + do j =1, knon + rhod=p1lay(j)/RD/temp_air(j) + nunu=max(ustar(j)/ustart(j),1.e-3) + fluxsalt=rhod/RG*(ustar(j)**3)*(1.-nunu**(-2)) * & (aa+bb*nunu**(-2)+cc*nunu**(-1)) - csalt=fluxsalt/(2.8*ustart(i)) - hsalt=0.08436*ustar(i)**1.27 - qsalt(i)=1./rhod*csalt*lambdasalt*RG/(max(ustar(i)**2,1E-6)) & - * exp(-lambdasalt*RG*hsalt/max(ustar(i)**2,1E-6)) - qsalt(i)=max(qsalt(i),0.) + csalt=fluxsalt/(2.8*ustart(j)) + hsalt(j)=0.08436*ustar(j)**1.27 + qsalt(j)=1./rhod*csalt*lambdasalt*RG/(max(ustar(j)**2,1E-6)) & + * exp(-lambdasalt*RG*hsalt(j)/max(ustar(j)**2,1E-6)) + qsalt(j)=max(qsalt(j),0.) enddo @@ -412,28 +416,33 @@ else ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001) - do i=1, knon - esalt=1./(3.25*max(ustar(i),0.001)) - hsalt=0.08436*ustar(i)**1.27 - qsalt(i)=(max(ustar(i)**2-ustart(i)**2,0.))/(RG*hsalt)*esalt - !ep=qsalt*cdragm(i)*sqrt(u1(i)**2+v1(i)**2) + do j=1, knon + esalt=1./(3.25*max(ustar(j),0.001)) + hsalt(j)=0.08436*ustar(j)**1.27 + qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt + !ep=qsalt*cdragm(j)*sqrt(u1(j)**2+v1(j)**2) enddo endif - ! calculation of erosion (emission flux towards the first atmospheric level) + ! calculation of erosion (flux positive towards the surface here) ! consistent with implicit resolution of turbulent mixing equation - do i=1, knon - rhod=p1lay(i)/RD/temp_air(i) - fluxbs(i)=rhod*ws1(i)*cdragm(i)*zeta_bs*(AcoefQBS(i)-qsalt(i)) & - / (1.-rhod*ws1(i)*zeta_bs*cdragm(i)*BcoefQBS(i)*dtime) - !fluxbs(i)= zeta_bs*rhod*ws1(i)*cdragm(i)*(qbs1(i)-qsalt(i)) + do j=1, knon + i = knindex(j) + rhod=p1lay(j)/RD/temp_air(j) + ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within an interval of about 10s + ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level) + ! to exceed (in intensity) integrated over 10s to exceed the total mass of snow particle in the saltation layer + ! (rho*qsalt*hsalt) + maxerosion=hsalt(j)*qsalt(j)*rhod/10. + fluxbs(j)=rhod*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) & + / (1.-rhod*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime) + !fluxbs(j)= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j)) + fluxbs(j)=max(-maxerosion,fluxbs(j)) + + ! for outputs + + zxustartlic(i) = ustart(j) + zxrhoslic(i) = rhos(j) enddo - ! for outputs - do j = 1, knon - i = knindex(j) - zxustartlic(i) = ustart(j) - zxrhoslic(i) = rhos(j) - enddo - endif @@ -441,5 +450,5 @@ !**************************************************************************************** -! Calculate surface snow amount +! Calculate snow amount ! !**************************************************************************************** @@ -451,9 +460,11 @@ evap_totsnow(:)=evap(:) ENDIF - + + CALL fonte_neige(knon, is_lic, knindex, dtime, & tsurf, precip_rain, precip_totsnow, & snow, qsol, tsurf_new, evap_totsnow) - + + WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0. zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0))) Index: LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90 (revision 4671) +++ LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90 (revision 4672) @@ -185,8 +185,10 @@ REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow REAL, DIMENSION (klon,6) :: alb6 - REAL :: rho0, rhoice, ustart0, hsalt, esalt, rhod + REAL :: rho0, rhoice, ustart0,esalt, rhod REAL :: lambdasalt,fluxsalt, csalt, nunu, aa, bb, cc REAL :: tau_dens, tau_dens0, tau_densmin, rhomax, rhohard - REAL, DIMENSION(klon) :: ws1, rhos, ustart, qsalt + REAL, DIMENSION(klon) :: ws1, rhos, ustart, qsalt, hsalt, snowini + REAL :: maxerosion ! in kg/s + ! End definition !**************************************************************************************** @@ -453,22 +455,24 @@ ustart0 = 0.211 rhoice = 920.0 - rho0 = 200.0 + rho0 = 300.0 rhomax=450.0 - rhohard=400.0 + rhohard=450.0 tau_dens0=86400.0*10. ! 10 days by default, in s - tau_densmin=86400.0 ! 1 days according to in situ obs by C. Amory + tau_densmin=21600.0 ! 1/4 days according to in situ obs by C. Amory during blowing snow + + ! Marshall et al. 1999 (snow densification during rain) ! computation of threshold friction velocity ! which depends on surface snow density - do i = 1, knon + do j = 1, knon ! estimation of snow density ! snow density increases with snow age and - ! increases even faster in case of sedimentation of blowing snow - tau_dens=max(tau_densmin, tau_dens0*exp(-abs(precip_bs(i))/pbst_bs-abs(precip_rain(i))/prt_bs)) - rhos(i)=rho0+(rhohard-rho0)*(1.-exp(-agesno(i)*86400.0/tau_dens)) + ! increases even faster in case of sedimentation of blowing snow or rain + tau_dens=max(tau_densmin, & tau_dens0*exp(-abs(precip_bs(j))/pbst_bs-abs(precip_rain(j))/prt_bs) & + *exp(-max(tsurf(j)-272.0,0.))) + rhos(j)=rho0+(rhohard-rho0)*(1.-exp(-agesno(j)*86400.0/tau_dens)) ! blowing snow flux formula used in MAR - ws1(i)=(u1(i)**2+v1(i)**2)**0.5 - ustar(i)=(cdragm(i)*(u1(i)**2+v1(i)**2))**0.5 - ustart(i)=ustart0*exp(max(rhoice/rho0-rhoice/rhos(i),0.))*exp(max(0.,rhos(i)-rhomax)) + ws1(j)=(u1(j)**2+v1(j)**2)**0.5 + ustar(j)=(cdragm(j)*(u1(j)**2+v1(j)**2))**0.5 + ustart(j)=ustart0*exp(max(rhoice/rho0-rhoice/rhos(j),0.))*exp(max(0.,rhos(j)-rhomax)) ! we have multiplied by exp to prevent erosion when rhos>rhomax (usefull till ! rhohard<450) @@ -486,14 +490,14 @@ cc=2.0 lambdasalt=0.45 - do i =1, knon - rhod=p1lay(i)/RD/temp_air(i) - nunu=max(ustar(i)/ustart(i),1.e-3) - fluxsalt=rhod/RG*(ustar(i)**3)*(1.-nunu**(-2)) * & + do j =1, knon + rhod=p1lay(j)/RD/temp_air(j) + nunu=max(ustar(j)/ustart(j),1.e-3) + fluxsalt=rhod/RG*(ustar(j)**3)*(1.-nunu**(-2)) * & (aa+bb*nunu**(-2)+cc*nunu**(-1)) - csalt=fluxsalt/(2.8*ustart(i)) - hsalt=0.08436*ustar(i)**1.27 - qsalt(i)=1./rhod*csalt*lambdasalt*RG/(max(ustar(i)**2,1E-6)) & - * exp(-lambdasalt*RG*hsalt/max(ustar(i)**2,1E-6)) - qsalt(i)=max(qsalt(i),0.) + csalt=fluxsalt/(2.8*ustart(j)) + hsalt(j)=0.08436*ustar(j)**1.27 + qsalt(j)=1./rhod*csalt*lambdasalt*RG/(max(ustar(j)**2,1E-6)) & + * exp(-lambdasalt*RG*hsalt(j)/max(ustar(j)**2,1E-6)) + qsalt(j)=max(qsalt(j),0.) enddo @@ -501,26 +505,31 @@ else ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001) - do i=1, knon - esalt=1./(3.25*max(ustar(i),0.001)) - hsalt=0.08436*ustar(i)**1.27 - qsalt(i)=(max(ustar(i)**2-ustart(i)**2,0.))/(RG*hsalt)*esalt - !ep=qsalt*cdragm(i)*sqrt(u1(i)**2+v1(i)**2) + do j=1, knon + esalt=1./(3.25*max(ustar(j),0.001)) + hsalt(j)=0.08436*ustar(j)**1.27 + qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt + !ep=qsalt*cdragm(j)*sqrt(u1(j)**2+v1(j)**2) enddo endif - ! calculation of erosion (emission flux towards the first atmospheric level) + ! calculation of erosion (flux positive towards the surface here) ! consistent with implicit resolution of turbulent mixing equation - do i=1, knon - rhod=p1lay(i)/RD/temp_air(i) - fluxbs(i)=rhod*ws1(i)*cdragm(i)*zeta_bs*(AcoefQBS(i)-qsalt(i)) & - / (1.-rhod*ws1(i)*zeta_bs*cdragm(i)*BcoefQBS(i)*dtime) - !fluxbs(i)= zeta_bs*rhod*ws1(i)*cdragm(i)*(qbs1(i)-qsalt(i)) - enddo - - ! for outputs - do j = 1, knon - i = knindex(j) - zxustartlic(i) = ustart(j) - zxrhoslic(i) = rhos(j) + do j=1, knon + i = knindex(j) + rhod=p1lay(j)/RD/temp_air(j) + ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within an interval of about 10s + ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level) + ! to exceed (in intensity) integrated over 10s to exceed the total mass of snow particle in the saltation layer + ! (rho*qsalt*hsalt) + maxerosion=hsalt(j)*qsalt(j)*rhod/10. + fluxbs(j)=rhod*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) & + / (1.-rhod*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime) + !fluxbs(j)= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j)) + fluxbs(j)=max(-maxerosion,fluxbs(j)) + + ! for outputs + + zxustartlic(i) = ustart(j) + zxrhoslic(i) = rhos(j) enddo