Ignore:
Timestamp:
Feb 29, 2024, 7:42:12 PM (3 months ago)
Author:
evignon
Message:

commission pour la suite du travail sur la mise a jour
de la param de neige soufflee

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90

    r4828 r4835  
    3131    USE cpl_mod,          ONLY : cpl_send_landice_fields
    3232    USE calcul_fluxs_mod
    33     USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic
     33    USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic, zxqsaltlic
    3434    USE phys_output_var_mod, ONLY : snow_o,zfra_o
    3535!FC
    3636    USE ioipsl_getin_p_mod, ONLY : getin_p
    37     USE lmdz_blowing_snow_ini, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs
    38 
     37    USE lmdz_blowing_snow_ini, ONLY : c_esalt_bs, zeta_bs, pbst_bs, prt_bs, rhoice_bs, rhohard_bs
     38    USE lmdz_blowing_snow_ini, ONLY : rhofresh_bs, tau_eqsalt_bs, tau_dens0_bs, tau_densmin_bs
    3939#ifdef CPP_INLANDSIS
    4040    USE surf_inlandsis_mod,  ONLY : surf_inlandsis
    … …  
    140140    REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow
    141141    REAL, DIMENSION (klon,6) :: alb6
    142     REAL                   :: rho0, rhoice, ustart0,esalt, rhod
     142    REAL                   :: esalt
    143143    REAL                   :: lambdasalt,fluxsalt, csalt, nunu, aa, bb, cc
    144     REAL                   :: tau_dens, tau_dens0, tau_densmin, rhomax, rhohard
    145     REAL, DIMENSION(klon)  :: ws1, rhos, ustart, qsalt, hsalt, snowini
    146     REAL                   :: ustarmin, maxerosion, tau_eq_salt
     144    REAL                   :: tau_dens, maxerosion, fluxbs_2
     145    REAL, DIMENSION(klon)  :: ws1, rhod, rhos, ustart0, ustart, qsalt, hsalt
    147146
    148147! End definition
    … …  
    331330!
    332331!****************************************************************************************
    333  
    334     CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)) 
    335 
    336 
    337 ! EV: following lines are obsolete since we set alb1 and alb2 to constant values
    338 ! I therefore comment them
    339 !    alb1(1:knon) = alb_neig(1:knon)*zfra(1:knon) + &
    340 !         0.6 * (1.0-zfra(1:knon))
     332
    341333!
    342334!IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
    … …  
    361353
    362354
    363 
    364   ! Simple blowing snow param
    365   if (ok_bs) then
    366        ustarmin=1.e-3
    367        ustart0 = 0.211
    368        rhoice = 920.0
    369        rho0 = 300.0
    370        rhomax=450.0
    371        rhohard=450.0
    372        tau_dens0=86400.0*10.  ! 10 days by default, in s
    373        tau_densmin=21600.0    ! 1/4 days according to in situ obs by C. Amory during blowing snow +
    374                               ! Marshall et al. 1999 (snow densification during rain)
    375        tau_eq_salt=10.
    376 
     355!****************************************************************************************
     356! Simple blowing snow param
     357!****************************************************************************************
     358! we proceed in 2 steps:
     359! first we erode - if possible -the accumulated snow during the time step
     360! then we update the density of the underlying layer and see if we can also erode
     361! this layer
     362
     363
     364   if (ok_bs) then
     365       fluxbs(:)=0.
     366       do j=1,klon
     367          ws1(j)=(u1(j)**2+v1(j)**2)**0.5
     368          ustar(j)=(cdragm(j)*(u1(j)**2+v1(j)**2))**0.5
     369          rhod(j)=p1lay(j)/RD/temp_air(j)
     370          ustart0(j) =(log(2.868)-log(1.625))/0.085*sqrt(cdragm(j))
     371       enddo
     372
     373       ! 1st step: erosion of fresh snow accumulated during the time step
     374       do j=1, knon
     375
     376           rhos(j)=rhofresh_bs
     377           ! blowing snow flux formula used in MAR
     378           ustart(j)=ustart0(j)*exp(max(rhoice_bs/rhofresh_bs-rhoice_bs/rhos(j),0.))*exp(max(0.,rhos(j)-rhohard_bs))
     379           ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs
     380           ! computation of qbs at the top of the saltation layer
     381           ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)
     382           esalt=1./(c_esalt_bs*max(1.e-6,ustar(j)))
     383           hsalt(j)=0.08436*(max(1.e-6,ustar(j))**1.27)
     384           qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt
     385           ! calculation of erosion (flux positive towards the surface here)
     386           ! consistent with implicit resolution of turbulent mixing equation
     387           ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s
     388           ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
     389           ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer
     390           ! (rho*qsalt*hsalt)
     391           ! during this first step we also lower bound the erosion to the amount of fresh snow accumulated during the time step
     392           maxerosion=min(precip_snow(j),hsalt(j)*qsalt(j)*rhod(j)/tau_eqsalt_bs)
     393
     394           fluxbs(j)=rhod(j)*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) &
     395                   / (1.-rhod(j)*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime)
     396           !fluxbs(j)= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j))
     397           fluxbs(j)=max(-maxerosion,fluxbs(j))
     398       enddo
     399
     400
     401       ! we now compute the snow age of the overlying layer (snow surface after erosion of the fresh snow accumulated during the time step)
     402       ! this is done through the routine albsno
     403       CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)+fluxbs(:))
     404
     405       ! 2nd step:
    377406       ! computation of threshold friction velocity
    378407       ! which depends on surface snow density
    … …  
    381410           ! snow density increases with snow age and
    382411           ! increases even faster in case of sedimentation of blowing snow or rain
    383            tau_dens=max(tau_densmin, tau_dens0*exp(-abs(precip_bs(j))/pbst_bs - &
    384                     abs(precip_rain(j))/prt_bs)*exp(-max(tsurf(j)-272.0,0.)))
    385            rhos(j)=rho0+(rhohard-rho0)*(1.-exp(-agesno(j)*86400.0/tau_dens))
     412           tau_dens=max(tau_densmin_bs, tau_dens0_bs*exp(-abs(precip_bs(j))/pbst_bs - &
     413                    abs(precip_rain(j))/prt_bs)*exp(-max(tsurf(j)-RTT,0.)))
     414           rhos(j)=rhofresh_bs+(rhohard_bs-rhofresh_bs)*(1.-exp(-agesno(j)*86400.0/tau_dens))
    386415           ! blowing snow flux formula used in MAR
    387            ws1(j)=(u1(j)**2+v1(j)**2)**0.5
    388            ustar(j)=(cdragm(j)*(u1(j)**2+v1(j)**2))**0.5
    389            ustart(j)=ustart0*exp(max(rhoice/rho0-rhoice/rhos(j),0.))*exp(max(0.,rhos(j)-rhomax))
    390            ! we have multiplied by exp to prevent erosion when rhos>rhomax (usefull till
    391            ! rhohard<450)
     416           ustart(j)=ustart0(j)*exp(max(rhoice_bs/rhofresh_bs-rhoice_bs/rhos(j),0.))*exp(max(0.,rhos(j)-rhohard_bs))
     417           ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs
     418           ! computation of qbs at the top of the saltation layer
     419           ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)
     420           esalt=c_esalt_bs*ustar(j)
     421           hsalt(j)=0.08436*(max(1.e-6,ustar(j))**1.27)
     422           qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt
     423           ! calculation of erosion (flux positive towards the surface here)
     424           ! consistent with implicit resolution of turbulent mixing equation
     425           ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s
     426           ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
     427           ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer
     428           ! (rho*qsalt*hsalt)
     429           maxerosion=hsalt(j)*qsalt(j)*rhod(j)/tau_eqsalt_bs
     430           fluxbs_2=rhod(j)*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) &
     431                   / (1.-rhod(j)*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime)
     432           !fluxbs_2= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j))
     433           fluxbs_2=max(-maxerosion,fluxbs_2)
     434           fluxbs(j)=fluxbs(j)+fluxbs_2
    392435       enddo
    393        
    394        ! computation of qbs at the top of the saltation layer
    395        ! two formulations possible
    396        ! pay attention that qbs is a mixing ratio and has to be converted
    397        ! to specific content
    398        
    399        if (iflag_saltation_bs .eq. 1) then
    400        ! expression from CRYOWRF (Sharma et al. 2022)
    401           aa=2.6
    402           bb=2.5
    403           cc=2.0
    404           lambdasalt=0.45
    405           do j =1, knon
    406                rhod=p1lay(j)/RD/temp_air(j)
    407                nunu=max(ustar(j)/ustart(j),1.e-3)
    408                fluxsalt=rhod/RG*(ustar(j)**3)*(1.-nunu**(-2)) * &
    409                         (aa+bb*nunu**(-2)+cc*nunu**(-1)) 
    410                csalt=fluxsalt/(2.8*ustart(j))
    411                hsalt(j)=0.08436*ustar(j)**1.27
    412                qsalt(j)=1./rhod*csalt*lambdasalt*RG/(max(ustar(j)**2,ustarmin)) &
    413                        * exp(-lambdasalt*RG*hsalt(j)/max(ustar(j)**2,ustarmin))
    414                qsalt(j)=max(qsalt(j),0.)
    415           enddo
    416 
    417 
    418        else
    419        ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)       
    420           do j=1, knon
    421               esalt=1./(3.25*max(ustarmin,ustar(j)))
    422               hsalt(j)=0.08436*(max(ustarmin,ustar(j))**1.27)
    423               qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt
    424               !ep=qsalt*cdragm(j)*sqrt(u1(j)**2+v1(j)**2)
    425           enddo
    426        endif
    427 
    428         ! calculation of erosion (flux positive towards the surface here)
    429         ! consistent with implicit resolution of turbulent mixing equation
    430        do j=1, knon
     436
     437       ! for outputs       
     438        do j=1, knon
    431439              i = knindex(j)
    432               rhod=p1lay(j)/RD/temp_air(j)
    433               ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eq_salt of about 10s
    434               ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
    435               ! integrated over tau_eq_salt to exceed the total mass of snow particle in the saltation layer
    436               ! (rho*qsalt*hsalt)
    437               maxerosion=hsalt(j)*qsalt(j)*rhod/tau_eq_salt
    438               fluxbs(j)=rhod*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) &
    439                        / (1.-rhod*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime)
    440               !fluxbs(j)= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j))
    441               fluxbs(j)=min(0.,max(-maxerosion,fluxbs(j)))
    442               ! for outputs
    443 
    444440              zxustartlic(i) = ustart(j)
    445441              zxrhoslic(i) = rhos(j)
    446        enddo
    447 
    448   endif
     442              zxqsaltlic(i)=qsalt(j)
     443        enddo
     444
     445
     446  else
     447  ! those lines are useful to calculate the snow age
     448       CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))
     449
     450  endif ! if ok_bs
    449451
    450452
Note: See TracChangeset for help on using the changeset viewer.