Index: LMDZ6/trunk/libf/phylmd/blowing_snow_ini_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/blowing_snow_ini_mod.F90 (revision 4528) +++ LMDZ6/trunk/libf/phylmd/blowing_snow_ini_mod.F90 (revision 4529) @@ -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 + integer :: niter_bs, 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) + !$OMP THREADPRIVATE(niter_bs, iflag_saltation_bs) contains @@ -67,4 +67,7 @@ CALL getin_p('niter_bs',niter_bs) + iflag_saltation_bs=0 + CALL getin_p('iflag_saltation_bs',iflag_saltation_bs) + end subroutine blowing_snow_ini Index: LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90 (revision 4528) +++ LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90 (revision 4529) @@ -2147,4 +2147,5 @@ AcoefH, AcoefQ, BcoefH, BcoefQ, & AcoefU, AcoefV, BcoefU, BcoefV, & + AcoefQBS, BcoefQBS, & ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, & ysnow, yqsurf, yqsol,yqbs1, yagesno, & Index: LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90 (revision 4528) +++ LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90 (revision 4529) @@ -15,4 +15,5 @@ AcoefH, AcoefQ, BcoefH, BcoefQ, & AcoefU, AcoefV, BcoefU, BcoefV, & + AcoefQBS, BcoefQBS, & ps, u1, v1, gustiness, rugoro, pctsrf, & snow, qsurf, qsol, qbs1, agesno, & @@ -34,5 +35,5 @@ !FC USE ioipsl_getin_p_mod, ONLY : getin_p - USE blowing_snow_ini_mod, ONLY : zeta_bs, pbst_bs, prt_bs + USE blowing_snow_ini_mod, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs #ifdef CPP_INLANDSIS @@ -62,4 +63,5 @@ REAL, DIMENSION(klon), INTENT(IN) :: BcoefH, BcoefQ REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV + REAL, DIMENSION(klon), INTENT(IN) :: AcoefQBS, BcoefQBS REAL, DIMENSION(klon), INTENT(IN) :: ps REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness, qbs1 @@ -138,7 +140,8 @@ REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow REAL, DIMENSION (klon,6) :: alb6 - REAL :: rho0, rhoice, ustart0, hsalt, esalt, qsalt + REAL :: rho0, rhoice, ustart0, hsalt, 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 + REAL, DIMENSION(klon) :: ws1, rhos, ustart, qsalt ! End definition !**************************************************************************************** @@ -366,4 +369,7 @@ 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 + + ! computation of threshold friction velocity + ! which depends on surface snow density do i = 1, knon ! estimation of snow density @@ -378,12 +384,47 @@ ! we have multiplied by exp to prevent erosion when rhos>rhomax (usefull till ! rhohard<450) - esalt=1./(3.25*max(ustar(i),0.001)) - hsalt=0.08436*ustar(i)**1.27 - qsalt=(max(ustar(i)**2-ustart(i)**2,0.))/(RG*hsalt)*esalt - !ep=qsalt*cdragm(i)*sqrt(u1(i)**2+v1(i)**2) - !Pay attention that this is an explicit formulation for the surface flux - !Formulation fully consistent with the implicit resolution of the turbulent - !mixing scheme remains to be coded - fluxbs(i)= zeta_bs*p1lay(i)/RD/temp_air(i)*ws1(i)*cdragm(i)*(qbs1(i)-qsalt) + enddo + + ! computation of qbs at the top of the saltation layer + ! two formulations possible + ! pay attention that qbs is a mixing ratio and has to be converted + ! to specific content + + if (iflag_saltation_bs .eq. 1) then + ! expression from CRYOWRF (Sharma et al. 2022) + aa=2.6 + bb=2.5 + 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)) * & + (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.) + enddo + + + 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) + enddo + endif + + ! calculation of erosion (emission flux towards the first atmospheric level) + ! 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*BcoefQBS(i)*dtime) + !fluxbs(i)= zeta_bs*rhod*ws1(i)*cdragm(i)*(qbs1(i)-qsalt(i)) enddo Index: LMDZ6/trunk/libf/phylmdiso/pbl_surface_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/pbl_surface_mod.F90 (revision 4528) +++ LMDZ6/trunk/libf/phylmdiso/pbl_surface_mod.F90 (revision 4529) @@ -2474,4 +2474,5 @@ AcoefH, AcoefQ, BcoefH, BcoefQ, & AcoefU, AcoefV, BcoefU, BcoefV, & + AcoefQBS, BcoefQBS, & ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, & ysnow, yqsurf, yqsol,yqbs1, yagesno, & Index: LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90 (revision 4528) +++ LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90 (revision 4529) @@ -15,4 +15,5 @@ AcoefH, AcoefQ, BcoefH, BcoefQ, & AcoefU, AcoefV, BcoefU, BcoefV, & + AcoefQBS, BcoefQBS, & ps, u1, v1, gustiness, rugoro, pctsrf, & snow, qsurf, qsol, qbs1, agesno, & @@ -49,5 +50,5 @@ !FC USE ioipsl_getin_p_mod, ONLY : getin_p - USE blowing_snow_ini_mod, ONLY : zeta_bs, pbst_bs, prt_bs + USE blowing_snow_ini_mod, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs #ifdef CPP_INLANDSIS @@ -78,4 +79,5 @@ REAL, DIMENSION(klon), INTENT(IN) :: BcoefH, BcoefQ REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV + REAL, DIMENSION(klon), INTENT(IN) :: AcoefQBS, BcoefQBS REAL, DIMENSION(klon), INTENT(IN) :: ps REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness, qbs1 @@ -183,7 +185,8 @@ REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow REAL, DIMENSION (klon,6) :: alb6 - REAL :: rho0, rhoice, ustart0, hsalt, esalt, qsalt + REAL :: rho0, rhoice, ustart0, hsalt, 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 + REAL, DIMENSION(klon) :: ws1, rhos, ustart, qsalt ! End definition !**************************************************************************************** @@ -455,4 +458,7 @@ 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 + + ! computation of threshold friction velocity + ! which depends on surface snow density do i = 1, knon ! estimation of snow density @@ -467,9 +473,47 @@ ! we have multiplied by exp to prevent erosion when rhos>rhomax (usefull till ! rhohard<450) - esalt=1./(3.25*max(ustar(i),0.001)) - hsalt=0.08436*ustar(i)**1.27 - qsalt=(max(ustar(i)**2-ustart(i)**2,0.))/(RG*hsalt)*esalt - !ep=qsalt*cdragm(i)*sqrt(u1(i)**2+v1(i)**2) - fluxbs(i)= zeta_bs*p1lay(i)/RD/temp_air(i)*ws1(i)*cdragm(i)*(qbs1(i)-qsalt) + enddo + + ! computation of qbs at the top of the saltation layer + ! two formulations possible + ! pay attention that qbs is a mixing ratio and has to be converted + ! to specific content + + if (iflag_saltation_bs .eq. 1) then + ! expression from CRYOWRF (Sharma et al. 2022) + aa=2.6 + bb=2.5 + 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)) * & + (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.) + enddo + + + 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) + enddo + endif + + ! calculation of erosion (emission flux towards the first atmospheric level) + ! 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*BcoefQBS(i)*dtime) + !fluxbs(i)= zeta_bs*rhod*ws1(i)*cdragm(i)*(qbs1(i)-qsalt(i)) enddo