Ignore:
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:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/LMDZ_ECRad

  • LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/surf_landice_mod.F90

    r4482 r4727  
    1212       rmu0, lwdownm, albedo, pphi1, &
    1313       swnet, lwnet, tsurf, p1lay, &
    14        cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
     14       cdragh, cdragm, precip_rain, precip_snow, precip_bs, temp_air, spechum, &
    1515       AcoefH, AcoefQ, BcoefH, BcoefQ, &
    1616       AcoefU, AcoefV, BcoefU, BcoefV, &
     17       AcoefQBS, BcoefQBS, &
    1718       ps, u1, v1, gustiness, rugoro, pctsrf, &
    18        snow, qsurf, qsol, agesno, &
    19        tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, fluxsens, fluxlat, &
     19       snow, qsurf, qsol, qbs1, agesno, &
     20       tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, fluxsens, fluxlat, fluxbs, &
    2021       tsurf_new, dflux_s, dflux_l, &
    2122       alt, slope, cloudf, &
    … …  
    3435    USE cpl_mod,          ONLY : cpl_send_landice_fields
    3536    USE calcul_fluxs_mod
    36     USE phys_output_var_mod
     37    USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic
     38    USE phys_output_var_mod, ONLY : snow_o,zfra_o
    3739#ifdef ISO   
    3840    USE fonte_neige_mod,  ONLY : xtrun_off_lic
    … …  
    4850!FC
    4951    USE ioipsl_getin_p_mod, ONLY : getin_p
    50 
     52    USE lmdz_blowing_snow_ini, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs
    5153
    5254#ifdef CPP_INLANDSIS
    … …  
    7274    REAL, DIMENSION(klon), INTENT(IN)             :: p1lay
    7375    REAL, DIMENSION(klon), INTENT(IN)             :: cdragh, cdragm
    74     REAL, DIMENSION(klon), INTENT(IN)             :: precip_rain, precip_snow
     76    REAL, DIMENSION(klon), INTENT(IN)             :: precip_rain, precip_snow, precip_bs
    7577    REAL, DIMENSION(klon), INTENT(IN)             :: temp_air, spechum
    7678    REAL, DIMENSION(klon), INTENT(IN)             :: AcoefH, AcoefQ
    7779    REAL, DIMENSION(klon), INTENT(IN)             :: BcoefH, BcoefQ
    7880    REAL, DIMENSION(klon), INTENT(IN)             :: AcoefU, AcoefV, BcoefU, BcoefV
     81    REAL, DIMENSION(klon), INTENT(IN)             :: AcoefQBS, BcoefQBS
    7982    REAL, DIMENSION(klon), INTENT(IN)             :: ps
    80     REAL, DIMENSION(klon), INTENT(IN)             :: u1, v1, gustiness
     83    REAL, DIMENSION(klon), INTENT(IN)             :: u1, v1, gustiness, qbs1
    8184    REAL, DIMENSION(klon), INTENT(IN)             :: rugoro
    8285    REAL, DIMENSION(klon,nbsrf), INTENT(IN)       :: pctsrf
    … …  
    118121!albedo SB <<<
    119122    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
     123    REAL, DIMENSION(klon), INTENT(OUT)            :: fluxbs
    120124    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
    121125    REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l     
    … …  
    179183
    180184    REAL,DIMENSION(klon) :: alb1,alb2
     185    REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow
    181186    REAL, DIMENSION (klon,6) :: alb6
     187    REAL                   :: rho0, rhoice, ustart0,esalt, rhod
     188    REAL                   :: lambdasalt,fluxsalt, csalt, nunu, aa, bb, cc
     189    REAL                   :: tau_dens, tau_dens0, tau_densmin, rhomax, rhohard
     190    REAL, DIMENSION(klon)  :: ws1, rhos, ustart, qsalt, hsalt, snowini
     191    REAL                   :: maxerosion ! in kg/s
     192
    182193! End definition
    183194!****************************************************************************************
    … …  
    214225           PRINT*, 'alb_nir_sno_lic',alb_nir_sno_lic
    215226 
    216 !  z0m=1.e-3
    217 !  z0h = z0m
    218227  firstcall=.false.
    219228  ENDIF
    220229!******************************************************************************************
    221 !
     230
    222231! Initialize output variables
    223232    alb3(:) = 999999.
    224233    alb2(:) = 999999.
    225234    alb1(:) = 999999.
    226    
     235    fluxbs(:)=0. 
    227236    runoff(:) = 0.
    228237!****************************************************************************************
    … …  
    232241    radsol(:) = 0.0
    233242    radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
     243
     244!****************************************************************************************
    234245
    235246!****************************************************************************************
    … …  
    327338
    328339
    329 
    330340    ELSE
    331341
    … …  
    342352    IF (soil_model) THEN
    343353       CALL soil(dtime, is_lic, knon, snow, tsurf, qsol, &
    344     & longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux)
    345 
     354        & longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux)
    346355       cal(1:knon) = RCPD / soilcap(1:knon)
    347356       radsol(1:knon)  = radsol(1:knon) + soilflux(1:knon)
    … …  
    406415         flux_u1, flux_v1)
    407416
    408 !****************************************************************************************
    409 ! Calculate snow height, age, run-off,..
     417
     418!****************************************************************************************
     419! Calculate albedo
     420!
     421!****************************************************************************************
     422 
     423    CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)) 
     424
     425
     426! EV: following lines are obsolete since we set alb1 and alb2 to constant values
     427! I therefore comment them
     428!    alb1(1:knon) = alb_neig(1:knon)*zfra(1:knon) + &
     429!         0.6 * (1.0-zfra(1:knon))
     430!
     431!IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
     432!       alb1(1 : knon)  = 0.6 !IM cf FH/GK
     433!       alb1(1 : knon)  = 0.82
     434!       alb1(1 : knon)  = 0.77 !211003 Ksta0.77
     435!       alb1(1 : knon)  = 0.8 !KstaTER0.8 & LMD_ARMIP5
     436!IM: KstaTER0.77 & LMD_ARMIP6   
     437
     438! Attantion: alb1 and alb2 are not the same!
     439    alb1(1:knon)  = alb_vis_sno_lic
     440    alb2(1:knon)  = alb_nir_sno_lic
     441
     442
     443!****************************************************************************************
     444! Rugosity
     445!
     446!****************************************************************************************
     447    z0m = z0m_landice
     448    z0h = z0h_landice
     449    !z0m = SQRT(z0m**2+rugoro**2)
     450
     451
     452
     453  ! Simple blowing snow param
     454  if (ok_bs) then
     455       ustart0 = 0.211
     456       rhoice = 920.0
     457       rho0 = 300.0
     458       rhomax=450.0
     459       rhohard=450.0
     460       tau_dens0=86400.0*10.  ! 10 days by default, in s
     461       tau_densmin=21600.0    ! 1/4 days according to in situ obs by C. Amory during blowing snow +
     462                              ! Marshall et al. 1999 (snow densification during rain)
     463
     464       ! computation of threshold friction velocity
     465       ! which depends on surface snow density
     466       do j = 1, knon
     467           ! estimation of snow density
     468           ! snow density increases with snow age and
     469           ! increases even faster in case of sedimentation of blowing snow or rain
     470           tau_dens=max(tau_densmin, tau_dens0*exp(-abs(precip_bs(j))/pbst_bs &
     471                   -abs(precip_rain(j))/prt_bs) *exp(-max(tsurf(j)-272.0,0.)))
     472           rhos(j)=rho0+(rhohard-rho0)*(1.-exp(-agesno(j)*86400.0/tau_dens))
     473           ! blowing snow flux formula used in MAR
     474           ws1(j)=(u1(j)**2+v1(j)**2)**0.5
     475           ustar(j)=(cdragm(j)*(u1(j)**2+v1(j)**2))**0.5
     476           ustart(j)=ustart0*exp(max(rhoice/rho0-rhoice/rhos(j),0.))*exp(max(0.,rhos(j)-rhomax))
     477           ! we have multiplied by exp to prevent erosion when rhos>rhomax (usefull till
     478           ! rhohard<450)
     479       enddo
     480       
     481       ! computation of qbs at the top of the saltation layer
     482       ! two formulations possible
     483       ! pay attention that qbs is a mixing ratio and has to be converted
     484       ! to specific content
     485       
     486       if (iflag_saltation_bs .eq. 1) then
     487       ! expression from CRYOWRF (Sharma et al. 2022)
     488          aa=2.6
     489          bb=2.5
     490          cc=2.0
     491          lambdasalt=0.45
     492          do j =1, knon
     493               rhod=p1lay(j)/RD/temp_air(j)
     494               nunu=max(ustar(j)/ustart(j),1.e-3)
     495               fluxsalt=rhod/RG*(ustar(j)**3)*(1.-nunu**(-2)) * &
     496                        (aa+bb*nunu**(-2)+cc*nunu**(-1)) 
     497               csalt=fluxsalt/(2.8*ustart(j))
     498               hsalt(j)=0.08436*ustar(j)**1.27
     499               qsalt(j)=1./rhod*csalt*lambdasalt*RG/(max(ustar(j)**2,1E-6)) &
     500                       * exp(-lambdasalt*RG*hsalt(j)/max(ustar(j)**2,1E-6))
     501               qsalt(j)=max(qsalt(j),0.)
     502          enddo
     503
     504
     505       else
     506       ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)       
     507          do j=1, knon
     508              esalt=1./(3.25*max(ustar(j),0.001))
     509              hsalt(j)=0.08436*ustar(j)**1.27
     510              qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt
     511              !ep=qsalt*cdragm(j)*sqrt(u1(j)**2+v1(j)**2)
     512          enddo
     513       endif
     514
     515        ! calculation of erosion (flux positive towards the surface here)
     516        ! consistent with implicit resolution of turbulent mixing equation
     517       do j=1, knon
     518              i = knindex(j)
     519              rhod=p1lay(j)/RD/temp_air(j)
     520              ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within an interval of about 10s
     521              ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
     522              ! to exceed (in intensity) integrated over 10s to exceed the total mass of snow particle in the saltation layer
     523              ! (rho*qsalt*hsalt)
     524              maxerosion=hsalt(j)*qsalt(j)*rhod/10.
     525              fluxbs(j)=rhod*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) &
     526                       / (1.-rhod*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime)
     527              !fluxbs(j)= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j))
     528              fluxbs(j)=max(-maxerosion,fluxbs(j))
     529
     530              ! for outputs
     531
     532              zxustartlic(i) = ustart(j)
     533              zxrhoslic(i) = rhos(j)
     534       enddo
     535
     536  endif
     537
     538
     539
     540!****************************************************************************************
     541! Calculate surface snow amount
    410542!   
    411543!****************************************************************************************
     544    IF (ok_bs) THEN
     545      precip_totsnow(:)=precip_snow(:)+precip_bs(:)
     546      evap_totsnow(:)=evap(:)-fluxbs(:) ! flux bs is positive towards the surface (snow erosion)
     547    ELSE
     548      precip_totsnow(:)=precip_snow(:)
     549      evap_totsnow(:)=evap(:)
     550    ENDIF
     551
    412552    CALL fonte_neige(knon, is_lic, knindex, dtime, &
    413          tsurf, precip_rain, precip_snow, &
    414          snow, qsol, tsurf_new, evap &
     553         tsurf, precip_rain, precip_totsnow, &
     554         snow, qsol, tsurf_new, evap_totsnow &
    415555#ifdef ISO   
    416556     & ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag   &
    … …  
    446586
    447587
    448 !****************************************************************************************
    449 ! Calculate albedo
    450 !
    451 !****************************************************************************************
    452     CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)) 
    453 
    454     WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
    455     zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))
    456     alb1(1:knon) = alb_neig(1:knon)*zfra(1:knon) + &
    457          0.6 * (1.0-zfra(1:knon))
    458 !
    459 !IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
    460 !       alb1(1 : knon)  = 0.6 !IM cf FH/GK
    461 !       alb1(1 : knon)  = 0.82
    462 !       alb1(1 : knon)  = 0.77 !211003 Ksta0.77
    463 !       alb1(1 : knon)  = 0.8 !KstaTER0.8 & LMD_ARMIP5
    464 !IM: KstaTER0.77 & LMD_ARMIP6   
    465 
    466 ! Attantion: alb1 and alb2 are not the same!
    467     alb1(1:knon)  = alb_vis_sno_lic
    468     alb2(1:knon)  = alb_nir_sno_lic
    469 
    470 
    471 !****************************************************************************************
    472 ! Rugosity
    473 !
    474 !****************************************************************************************
    475     z0m=1.e-3
    476     z0h = z0m
    477     z0m = SQRT(z0m**2+rugoro**2)
    478 
     588    WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.                                         
     589    zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0))) 
    479590
    480591
    … …  
    494605          run_off_lic_frac(j) = pctsrf(i,is_lic)
    495606       ENDDO
    496        
     607
    497608       CALL cpl_send_landice_fields(itime, knon, knindex, run_off_lic, run_off_lic_frac)
    498609    ENDIF
    … …  
    501612    runoff(1:knon)=run_off_lic(1:knon)/dtime
    502613
    503  
    504 !****************************************************************************************
    505 ! Etienne: comment these lines because of duplication just below
    506 !       snow_o=0.
    507 !       zfra_o = 0.
    508 !       DO j = 1, knon
    509 !           i = knindex(j)
    510 !           snow_o(i) = snow(j)
    511 !           zfra_o(i) = zfra(j)
    512 !       ENDDO
    513 !
    514 !****************************************************************************************
    515614       snow_o=0.
    516615       zfra_o = 0.
    … …  
    553652!albedo SB <<<
    554653
    555  
    556  
    557654
    558655  END SUBROUTINE surf_landice
Note: See TracChangeset for help on using the changeset viewer.