Changeset 4727 for LMDZ6/branches/LMDZ_ECRad/libf/phylmd/surf_landice_mod.F90

Timestamp:

Oct 19, 2023, 4:02:57 PM (8 months ago)

Author:

idelkadi

Message:

Merged trunk changes -r4488:4726 LMDZ_ECRad branch

Location:

LMDZ6/branches/LMDZ_ECRad

Files:

• . (modified) (1 prop)
• libf/phylmd/surf_landice_mod.F90 (modified) (13 diffs)

LMDZ6/branches/LMDZ_ECRad/libf/phylmd/surf_landice_mod.F90

@@ -12 +12 @@
        rmu0, lwdownm, albedo, pphi1, &
        swnet, lwnet, tsurf, p1lay, &
-       cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
+       cdragh, cdragm, precip_rain, precip_snow, precip_bs, temp_air, spechum, &
        AcoefH, AcoefQ, BcoefH, BcoefQ, &
        AcoefU, AcoefV, BcoefU, BcoefV, &
+       AcoefQBS, BcoefQBS, &
        ps, u1, v1, gustiness, rugoro, pctsrf, &
-       snow, qsurf, qsol, agesno, &
-       tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, fluxsens, fluxlat, &
+       snow, qsurf, qsol, qbs1, agesno, &
+       tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, fluxsens, fluxlat, fluxbs, &
        tsurf_new, dflux_s, dflux_l, &
        alt, slope, cloudf, &
@@ -30 +31 @@
     USE cpl_mod,          ONLY : cpl_send_landice_fields
     USE calcul_fluxs_mod
+    USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic
     USE phys_output_var_mod, ONLY : snow_o,zfra_o
 !FC
     USE ioipsl_getin_p_mod, ONLY : getin_p
-
+    USE lmdz_blowing_snow_ini, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs
 
 #ifdef CPP_INLANDSIS
@@ -56 +58 @@
     REAL, DIMENSION(klon), INTENT(IN)             :: p1lay
     REAL, DIMENSION(klon), INTENT(IN)             :: cdragh, cdragm
-    REAL, DIMENSION(klon), INTENT(IN)             :: precip_rain, precip_snow
+    REAL, DIMENSION(klon), INTENT(IN)             :: precip_rain, precip_snow, precip_bs
     REAL, DIMENSION(klon), INTENT(IN)             :: temp_air, spechum
     REAL, DIMENSION(klon), INTENT(IN)             :: AcoefH, AcoefQ
     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
+    REAL, DIMENSION(klon), INTENT(IN)             :: u1, v1, gustiness, qbs1
     REAL, DIMENSION(klon), INTENT(IN)             :: rugoro
     REAL, DIMENSION(klon,nbsrf), INTENT(IN)       :: pctsrf
@@ -94 +97 @@
 !albedo SB <<<
     REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
+    REAL, DIMENSION(klon), INTENT(OUT)            :: fluxbs
     REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l      
@@ -134 +138 @@
 
     REAL,DIMENSION(klon) :: alb1,alb2
+    REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow
     REAL, DIMENSION (klon,6) :: alb6
+    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, hsalt, snowini
+    REAL                   :: maxerosion ! in kg/s   
+
 ! End definition
 !****************************************************************************************
@@ -165 +176 @@
     alb2(:) = 999999.
     alb1(:) = 999999.
-    
+    fluxbs(:)=0.  
     runoff(:) = 0.
 !****************************************************************************************
@@ -173 +184 @@
     radsol(:) = 0.0
     radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
+
+!****************************************************************************************
 
 !****************************************************************************************
@@ -264 +277 @@
 
 
-
     ELSE 
 
@@ -314 +326 @@
          flux_u1, flux_v1)
 
-!****************************************************************************************
-! Calculate snow height, age, run-off,..
-!    
-!****************************************************************************************
-    CALL fonte_neige(knon, is_lic, knindex, dtime, &
-         tsurf, precip_rain, precip_snow, &
-         snow, qsol, tsurf_new, evap)
-
 
 !****************************************************************************************
@@ -327 +331 @@
 !
 !****************************************************************************************
+ 
     CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))  
 
-    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)))
-    alb1(1:knon) = alb_neig(1:knon)*zfra(1:knon) + &
-         0.6 * (1.0-zfra(1:knon))
+
+! EV: following lines are obsolete since we set alb1 and alb2 to constant values
+! I therefore comment them
+!    alb1(1:knon) = alb_neig(1:knon)*zfra(1:knon) + &
+!         0.6 * (1.0-zfra(1:knon))
 !
 !IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
@@ -354 +360 @@
     !z0m = SQRT(z0m**2+rugoro**2)
 
+
+
+  ! Simple blowing snow param
+  if (ok_bs) then
+       ustart0 = 0.211
+       rhoice = 920.0
+       rho0 = 300.0
+       rhomax=450.0
+       rhohard=450.0
+       tau_dens0=86400.0*10.  ! 10 days by default, in s
+       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 j = 1, knon
+           ! estimation of snow density
+           ! snow density increases with snow age and
+           ! 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(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)
+       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 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(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
+
+
+       else
+       ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)        
+          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 (flux positive towards the surface here) 
+        ! consistent with implicit resolution of turbulent mixing equation
+       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
+
+  endif
+
+
+
+!****************************************************************************************
+! Calculate snow amount
+!    
+!****************************************************************************************
+    IF (ok_bs) THEN
+      precip_totsnow(:)=precip_snow(:)+precip_bs(:)
+      evap_totsnow(:)=evap(:)-fluxbs(:) ! flux bs is positive towards the surface (snow erosion)
+    ELSE
+      precip_totsnow(:)=precip_snow(:)
+      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)))  
 
 
@@ -378 +493 @@
     runoff(1:knon)=run_off_lic(1:knon)/dtime
 
-  
-!****************************************************************************************
        snow_o=0.
        zfra_o = 0.
@@ -420 +533 @@
 !albedo SB <<<
 
-  
-  
 
   END SUBROUTINE surf_landice