Changeset 4506 for LMDZ6/branches


Ignore:
Timestamp:
Apr 11, 2023, 9:48:08 PM (13 months ago)
Author:
evignon
Message:

commit sur des modifications propres aux isotopes pour la neige soufflee

Location:
LMDZ6/branches/blowing_snow/libf/phylmdiso
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/blowing_snow/libf/phylmdiso/add_phys_tend_mod.F90

    r4143 r4506  
    1616CONTAINS
    1717
    18 SUBROUTINE add_pbl_tend(zdu, zdv, zdt, zdq, zdql, zdqi, paprs, text,abortphy,flag_inhib_tend, itap &
     18SUBROUTINE add_pbl_tend(zdu, zdv, zdt, zdq, zdql, zdqi, zdqbs, paprs, text,abortphy,flag_inhib_tend, itap &
    1919#ifdef ISO
    2020        ,zdxt,zdxtl,zdxti &
    … …  
    6262  ! ------------
    6363  REAL zdu(klon, klev), zdv(klon, klev)
    64   REAL zdt(klon, klev), zdq(klon, klev), zdql(klon, klev), zdqi(klon, klev)
     64  REAL zdt(klon, klev), zdq(klon, klev), zdql(klon, klev), zdqi(klon, klev), zdqbs(klon,klev)
    6565  CHARACTER *(*) text
    6666  REAL paprs(klon,klev+1)
    … …  
    104104    PRINT *, ' add_pbl_tend, zzdt ', zzdt
    105105    PRINT *, ' add_pbl_tend, zzdq ', zzdq
    106     CALL add_phys_tend(zdu, zdv, zzdt, zzdq, zdql, zdqi, paprs, text,abortphy,flag_inhib_tend, itap, 0 &
     106    CALL add_phys_tend(zdu, zdv, zzdt, zzdq, zdql, zdqi, zdqbs, paprs, text,abortphy,flag_inhib_tend, itap, 0 &
    107107#ifdef ISO
    108108        ,zzdxt,zdxtl,zdxti &
    … …  
    110110        )
    111111  ELSE
    112     CALL add_phys_tend(zdu, zdv, zdt, zdq, zdql, zdqi, paprs, text,abortphy,flag_inhib_tend, itap, 0 &
     112    CALL add_phys_tend(zdu, zdv, zdt, zdq, zdql, zdqi, zdqbs, paprs, text,abortphy,flag_inhib_tend, itap, 0 &
    113113#ifdef ISO
    114114        ,zdxt,zdxtl,zdxti &
    … …  
    123123! $Id$
    124124!
    125 SUBROUTINE add_phys_tend (zdu,zdv,zdt,zdq,zdql,zdqi,paprs,text, &
     125SUBROUTINE add_phys_tend (zdu,zdv,zdt,zdq,zdql,zdqi,zdqbs,paprs,text, &
    126126                          abortphy,flag_inhib_tend, itap, diag_mode &
    127127#ifdef ISO
    … …  
    142142USE dimphy, ONLY: klon, klev
    143143USE phys_state_var_mod, ONLY : phys_tstep
    144 USE phys_local_var_mod, ONLY: u_seri, v_seri, ql_seri, qs_seri, q_seri, t_seri
     144USE phys_local_var_mod, ONLY: u_seri, v_seri, ql_seri, qs_seri, qbs_seri, q_seri, t_seri
    145145#ifdef ISO
    146146USE phys_local_var_mod, ONLY: xtl_seri, xts_seri, xt_seri
    … …  
    150150USE print_control_mod, ONLY: prt_level
    151151USE cmp_seri_mod
    152 USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qt_col, d_ek_col, d_h_dair_col &
    153   &           , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_col
     152USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qbs_col, d_qt_col, d_ek_col, d_h_dair_col &
     153  &           , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_qbs_col, d_h_col
    154154
    155155#ifdef ISO
    … …  
    167167!------------
    168168REAL, DIMENSION(klon,klev),     INTENT(IN)    :: zdu, zdv
    169 REAL, DIMENSION(klon,klev),     INTENT(IN)    :: zdt, zdql, zdqi
     169REAL, DIMENSION(klon,klev),     INTENT(IN)    :: zdt, zdql, zdqi, zdqbs
    170170REAL, DIMENSION(klon,klev+1),   INTENT(IN)    :: paprs
    171171CHARACTER*(*),                  INTENT(IN)    :: text
    … …  
    197197! Save variables, used in diagnostic mode (diag_mode=1).
    198198REAL, DIMENSION(klon,klev)   :: sav_u_seri, sav_v_seri
    199 REAL, DIMENSION(klon,klev)   :: sav_ql_seri, sav_qs_seri, sav_q_seri
     199REAL, DIMENSION(klon,klev)   :: sav_ql_seri, sav_qs_seri, sav_qbs_seri, sav_q_seri
    200200REAL, DIMENSION(klon,klev)   :: sav_t_seri
    201201REAL, DIMENSION(klon,klev)   :: sav_zdq
    … …  
    228228! zh_ql_col----  total enthalpy of liquid watter (J/m2)
    229229! zh_qs_col----  total enthalpy of solid watter  (J/m2)
     230! zh_qbs_col----  total enthalpy of blowing snow  (J/m2)
    230231! zqw_col------  total mass of watter vapour (kg/m2)
    231232! zql_col------  total mass of liquid watter (kg/m2)
    232 ! zqs_col------  total mass of solid watter (kg/m2)
     233! zqs_col------  total mass of cloud ice (kg/m2)
     234! zqbs_col------  total mass of blowing snow (kg/m2)
    233235! zek_col------  total kinetic energy (kg/m2)
    234236!
    … …  
    237239REAL zql_col(klon,2)
    238240REAL zqs_col(klon,2)
     241REAL zqbs_col(klon,2)
    239242REAL zek_col(klon,2)
    240243REAL zh_dair_col(klon,2)
    241 REAL zh_qw_col(klon,2), zh_ql_col(klon,2), zh_qs_col(klon,2)
     244REAL zh_qw_col(klon,2), zh_ql_col(klon,2), zh_qs_col(klon,2), zh_qbs_col(klon,2)
    242245REAL zh_col(klon,2)
    243246
    … …  
    278281    sav_ql_seri(:,:) = ql_seri(:,:)
    279282    sav_qs_seri(:,:) = qs_seri(:,:)
     283    sav_qbs_seri(:,:) = qbs_seri(:,:)
    280284    sav_q_seri(:,:)  = q_seri(:,:)
    281285    sav_t_seri(:,:)  = t_seri(:,:)
    … …  
    304308
    305309    CALL integr_v(klon, klev, zcpvap, &
    306                   t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zairm, &
    307                     zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
    308                     zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n))
     310                  t_seri, q_seri, ql_seri, qs_seri, qbs_seri, u_seri, v_seri, zairm, &
     311                    zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zqbs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
     312                    zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_qbs_col(:,n), zh_col(:,n))
    309313
    310314  end if ! end if (fl_ebil .GT. 0)
    … …  
    318322     ql_seri(:,:)=ql_seri(:,:)+zdql(:,:)
    319323     qs_seri(:,:)=qs_seri(:,:)+zdqi(:,:)
     324     qbs_seri(:,:)=qbs_seri(:,:)+zdqbs(:,:)
    320325#ifdef ISO
    321326     xtl_seri(:,:,:)=xtl_seri(:,:,:)+zdxtl(:,:,:)
    … …  
    572577        ql_seri(:,:)=ql_seri(:,:)-zdql(:,:)
    573578        qs_seri(:,:)=qs_seri(:,:)-zdqi(:,:)
     579        qbs_seri(:,:)=qbs_seri(:,:)-zdqbs(:,:)
    574580      ENDIF
    575581
    … …  
    586592
    587593    CALL integr_v(klon, klev, zcpvap, &
    588                   t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zairm, &
    589                     zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
    590                     zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n))
     594                  t_seri, q_seri, ql_seri, qs_seri, qbs_seri, u_seri, v_seri, zairm, &
     595                    zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zqbs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
     596                    zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_qbs_col(:,n), zh_col(:,n))
    591597
    592598    ! ------------------------------------------------
    … …  
    597603    d_ql_col(:) = (zql_col(:,2)-zql_col(:,1))/phys_tstep
    598604    d_qs_col(:) = (zqs_col(:,2)-zqs_col(:,1))/phys_tstep
    599     d_qt_col(:) = d_qw_col(:) + d_ql_col(:) + d_qs_col(:)
     605    d_qbs_col(:) = (zqbs_col(:,2)-zqbs_col(:,1))/phys_tstep
     606    d_qt_col(:) = d_qw_col(:) + d_ql_col(:) + d_qs_col(:) + d_qbs_col(:)
    600607
    601608    d_ek_col(:) = (zek_col(:,2)-zek_col(:,1))/phys_tstep
    … …  
    605612    d_h_ql_col(:) = (zh_ql_col(:,2)-zh_ql_col(:,1))/phys_tstep
    606613    d_h_qs_col(:) = (zh_qs_col(:,2)-zh_qs_col(:,1))/phys_tstep
     614    d_h_qbs_col(:) = (zh_qbs_col(:,2)-zh_qbs_col(:,1))/phys_tstep
    607615
    608616    d_h_col = (zh_col(:,2)-zh_col(:,1))/phys_tstep
    … …  
    616624    ql_seri(:,:) = sav_ql_seri(:,:)
    617625    qs_seri(:,:) = sav_qs_seri(:,:)
     626    qbs_seri(:,:) = sav_qbs_seri(:,:)
    618627    q_seri(:,:)  = sav_q_seri(:,:)
    619628    t_seri(:,:)  = sav_t_seri(:,:)
    … …  
    630639END SUBROUTINE add_phys_tend
    631640
    632 SUBROUTINE diag_phys_tend (nlon, nlev, uu, vv, temp, qv, ql, qs, &
    633                           zdu,zdv,zdt,zdq,zdql,zdqs,paprs,text)
     641SUBROUTINE diag_phys_tend (nlon, nlev, uu, vv, temp, qv, ql, qs, qbs, &
     642                          zdu,zdv,zdt,zdq,zdql,zdqs,zdqbs,paprs,text)
    634643!======================================================================
    635644! Ajoute les tendances des variables physiques aux variables
    … …  
    647656USE print_control_mod, ONLY: prt_level
    648657USE cmp_seri_mod
    649 USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qt_col, d_ek_col, d_h_dair_col &
    650   &           , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_col
     658USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qbs_col, d_qt_col, d_ek_col, d_h_dair_col &
     659  &           , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_qbs_col, d_h_col
    651660IMPLICIT none
    652661  include "YOMCST.h"
    … …  
    657666INTEGER, INTENT(IN)                           :: nlon, nlev
    658667REAL, DIMENSION(nlon,nlev),     INTENT(IN)    :: uu, vv
    659 REAL, DIMENSION(nlon,nlev),     INTENT(IN)    :: temp, qv, ql, qs
     668REAL, DIMENSION(nlon,nlev),     INTENT(IN)    :: temp, qv, ql, qs, qbs
    660669REAL, DIMENSION(nlon,nlev),     INTENT(IN)    :: zdu, zdv
    661 REAL, DIMENSION(nlon,nlev),     INTENT(IN)    :: zdt, zdq, zdql, zdqs
     670REAL, DIMENSION(nlon,nlev),     INTENT(IN)    :: zdt, zdq, zdql, zdqs, zdqbs
    662671REAL, DIMENSION(nlon,nlev+1),   INTENT(IN)    :: paprs
    663672CHARACTER*(*),                  INTENT(IN)    :: text
    … …  
    666675!--------
    667676REAL, DIMENSION(nlon,nlev)      :: uu_n, vv_n
    668 REAL, DIMENSION(nlon,nlev)      :: temp_n, qv_n, ql_n, qs_n
     677REAL, DIMENSION(nlon,nlev)      :: temp_n, qv_n, ql_n, qs_n, qbs_n
    669678
    670679
    … …  
    689698! zqw_col------  total mass of watter vapour (kg/m2)
    690699! zql_col------  total mass of liquid watter (kg/m2)
    691 ! zqs_col------  total mass of solid watter (kg/m2)
     700! zqs_col------  total mass of cloud ice (kg/m2)
     701! zqbs_col------  total mass of blowing snow (kg/m2)
    692702! zek_col------  total kinetic energy (kg/m2)
    693703!
    … …  
    696706REAL zql_col(nlon,2)
    697707REAL zqs_col(nlon,2)
     708REAL zqbs_col(nlon,2)
    698709REAL zek_col(nlon,2)
    699710REAL zh_dair_col(nlon,2)
    700 REAL zh_qw_col(nlon,2), zh_ql_col(nlon,2), zh_qs_col(nlon,2)
     711REAL zh_qw_col(nlon,2), zh_ql_col(nlon,2), zh_qs_col(nlon,2), zh_qbs_col(nlon,2)
    701712REAL zh_col(nlon,2)
    702713
    … …  
    731742
    732743    CALL integr_v(nlon, nlev, rcpv, &
    733                   temp, qv, ql, qs, uu, vv, zairm, &
    734                     zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
    735                     zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n))
     744                  temp, qv, ql, qs, qbs, uu, vv, zairm, &
     745                    zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zqbs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
     746                    zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_qbs_col(:,n), zh_col(:,n))
    736747
    737748  end if ! end if (fl_ebil .GT. 0)
    … …  
    746757     ql_n(:,:)=ql(:,:)+zdql(:,:)
    747758     qs_n(:,:)=qs(:,:)+zdqs(:,:)
     759     qbs_n(:,:)=qbs(:,:)+zdqbs(:,:)
    748760     temp_n(:,:)=temp(:,:)+zdt(:,:)
    749761
    … …  
    762774
    763775    CALL integr_v(nlon, nlev, rcpv, &
    764                   temp_n, qv_n, ql_n, qs_n, uu_n, vv_n, zairm, &
    765                     zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
    766                     zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n))
     776                  temp_n, qv_n, ql_n, qs_n, qbs_n, uu_n, vv_n, zairm, &
     777                    zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zqbs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
     778                    zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_qbs_col(:,n), zh_col(:,n))
    767779
    768780    ! ------------------------------------------------
    … …  
    773785    d_ql_col(:) = (zql_col(:,2)-zql_col(:,1))/phys_tstep
    774786    d_qs_col(:) = (zqs_col(:,2)-zqs_col(:,1))/phys_tstep
    775     d_qt_col(:) = d_qw_col(:) + d_ql_col(:) + d_qs_col(:)
     787    d_qbs_col(:) = (zqbs_col(:,2)-zqbs_col(:,1))/phys_tstep
     788    d_qt_col(:) = d_qw_col(:) + d_ql_col(:) + d_qs_col(:) + d_qbs_col(:)
    776789
    777790    d_ek_col(:) = (zek_col(:,2)-zek_col(:,1))/phys_tstep
    … …  
    785798    d_h_ql_col(:) = (zh_ql_col(:,2)-zh_ql_col(:,1))/phys_tstep
    786799    d_h_qs_col(:) = (zh_qs_col(:,2)-zh_qs_col(:,1))/phys_tstep
     800    d_h_qbs_col(:) = (zh_qbs_col(:,2)-zh_qbs_col(:,1))/phys_tstep
    787801
    788802    d_h_col = (zh_col(:,2)-zh_col(:,1))/phys_tstep
    … …  
    795809
    796810SUBROUTINE integr_v(nlon, nlev, zcpvap, &
    797                     temp, qv, ql, qs, uu, vv, zairm,  &
    798                     zqw_col, zql_col, zqs_col, zek_col, zh_dair_col, &
    799                     zh_qw_col, zh_ql_col, zh_qs_col, zh_col)
     811                    temp, qv, ql, qs, qbs, uu, vv, zairm,  &
     812                    zqw_col, zql_col, zqs_col, zqbs_col, zek_col, zh_dair_col, &
     813                    zh_qw_col, zh_ql_col, zh_qs_col, zh_qbs_col, zh_col)
    800814
    801815IMPLICIT none
    … …  
    804818INTEGER,                    INTENT(IN)    :: nlon,nlev
    805819REAL,                       INTENT(IN)    :: zcpvap
    806 REAL, DIMENSION(nlon,nlev), INTENT(IN)    :: temp, qv, ql, qs, uu, vv
     820REAL, DIMENSION(nlon,nlev), INTENT(IN)    :: temp, qv, ql, qs, qbs, uu, vv
    807821REAL, DIMENSION(nlon,nlev), INTENT(IN)    :: zairm
    808822REAL, DIMENSION(nlon),      INTENT(OUT)   :: zqw_col
    809823REAL, DIMENSION(nlon),      INTENT(OUT)   :: zql_col
    810 REAL, DIMENSION(nlon),      INTENT(OUT)   :: zqs_col
     824REAL, DIMENSION(nlon),      INTENT(OUT)   :: zqs_col, zqbs_col
    811825REAL, DIMENSION(nlon),      INTENT(OUT)   :: zek_col
    812826REAL, DIMENSION(nlon),      INTENT(OUT)   :: zh_dair_col
    813827REAL, DIMENSION(nlon),      INTENT(OUT)   :: zh_qw_col
    814828REAL, DIMENSION(nlon),      INTENT(OUT)   :: zh_ql_col
    815 REAL, DIMENSION(nlon),      INTENT(OUT)   :: zh_qs_col
     829REAL, DIMENSION(nlon),      INTENT(OUT)   :: zh_qs_col, zh_qbs_col
    816830REAL, DIMENSION(nlon),      INTENT(OUT)   :: zh_col
    817831
    … …  
    823837    zql_col(:) = 0.
    824838    zqs_col(:) = 0.
     839    zqbs_col(:) = 0.
    825840    zek_col(:) = 0.
    826841    zh_dair_col(:) = 0.
    … …  
    828843    zh_ql_col(:) = 0.
    829844    zh_qs_col(:) = 0.
     845    zh_qbs_col(:) = 0.
    830846
    831847!JLD    write (*,*) "rcpd, zcpvap, zcwat, zcice ",rcpd, zcpvap, zcwat, zcice
    … …  
    840856        zql_col(i) = zql_col(i) + ql(i, k)*zairm(i, k)
    841857        zqs_col(i) = zqs_col(i) + qs(i, k)*zairm(i, k)
     858        zqbs_col(i)= zqbs_col(i) + qbs(i,k)*zairm(i,k)
    842859        ! Kinetic Energy
    843860        zek_col(i) = zek_col(i) + 0.5*(uu(i,k)**2+vv(i,k)**2)*zairm(i, k)
    … …  
    848865        zh_ql_col(i) = zh_ql_col(i) + (zcpvap*temp(i, k) - rlvtt)*ql(i, k)*zairm(i, k)   !jyg
    849866        zh_qs_col(i) = zh_qs_col(i) + (zcpvap*temp(i, k) - rlstt)*qs(i, k)*zairm(i, k)   !jyg
     867        zh_qbs_col(i) = zh_qbs_col(i) + (zcpvap*temp(i, k) - rlstt)*qbs(i, k)*zairm(i, k)   !jyg
    850868      END DO
    851869    END DO
    852870    ! compute total air enthalpy
    853     zh_col(:) = zh_dair_col(:) + zh_qw_col(:) + zh_ql_col(:) + zh_qs_col(:)
     871    zh_col(:) = zh_dair_col(:) + zh_qw_col(:) + zh_ql_col(:) + zh_qs_col(:) + zh_qbs_col(:)
    854872
    855873END SUBROUTINE integr_v
    … …  
    866884USE dimphy, ONLY: klon, klev
    867885USE phys_state_var_mod, ONLY : phys_tstep
    868 USE phys_state_var_mod, ONLY : topsw, toplw, solsw, sollw, rain_con, snow_con
     886USE phys_state_var_mod, ONLY : topsw, toplw, solsw, sollw, rain_con, snow_con, bs_fall
    869887USE geometry_mod, ONLY: longitude_deg, latitude_deg
    870888USE print_control_mod, ONLY: prt_level
    871889USE cmp_seri_mod
    872 USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qt_col, d_ek_col, d_h_dair_col &
    873   &           , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_col
     890USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qbs_col, d_qt_col, d_ek_col, d_h_dair_col &
     891  &           , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_qbs_col, d_h_col
    874892USE phys_local_var_mod, ONLY: evap, sens
    875 USE phys_local_var_mod, ONLY: u_seri, v_seri, ql_seri, qs_seri, q_seri, t_seri &
     893USE phys_local_var_mod, ONLY: u_seri, v_seri, ql_seri, qs_seri, qbs_seri, q_seri, t_seri &
    876894   &  , rain_lsc, snow_lsc
    877895USE climb_hq_mod, ONLY : d_h_col_vdf, f_h_bnd
    … …  
    910928      bilh_bnd = (-(rcw-rcpd)*t_seri(1,1) + rlvtt) * rain_lsc(1) &
    911929    &         + (-(rcs-rcpd)*t_seri(1,1) + rlstt) * snow_lsc(1)
     930  CASE("bs") param
     931      bilq_bnd = - bs_fall(1)
     932      bilh_bnd = (-(rcs-rcpd)*t_seri(1,1) + rlstt) * bs_fall(1)
    912933  CASE("convection") param
    913934      bilq_bnd = - rain_con(1) - snow_con(1)
    … …  
    946967  if ( prt_level .GE. 5) then
    947968    write(*,9000) text,"enerbil at boundaries: Q, H",bilq_bnd, bilh_bnd
    948     write(*,9000) text,"enerbil: water budget",d_qt_col(1),d_qw_col(1),d_ql_col(1),d_qs_col(1)
    949     write(*,9000) text,"enerbil: enthalpy budget",d_h_col(1),d_h_dair_col(1),d_h_qw_col(1),d_h_ql_col(1),d_h_qs_col(1)
     969    write(*,9000) text,"enerbil: water budget",d_qt_col(1),d_qw_col(1),d_ql_col(1),d_qs_col(1), d_qbs_col(1)
     970    write(*,9000) text,"enerbil: enthalpy budget",d_h_col(1),d_h_dair_col(1),d_h_qw_col(1),d_h_ql_col(1),d_h_qs_col(1),d_h_qbs_col(1)
    950971  end if
    951972

  • LMDZ6/branches/blowing_snow/libf/phylmdiso/pbl_surface_mod.F90

    r4478 r4506  
    2121  USE cpl_mod,             ONLY : gath2cpl
    2222  USE climb_hq_mod,        ONLY : climb_hq_down, climb_hq_up
     23  USE climb_qbs_mod,       ONLY : climb_qbs_down, climb_qbs_up
    2324  USE climb_wind_mod,      ONLY : climb_wind_down, climb_wind_up
    2425  USE coef_diff_turb_mod,  ONLY : coef_diff_turb
    … …  
    261262       rlon,      rlat,      rugoro,   rmu0,          &
    262263       lwdown_m,  cldt,          &
    263        rain_f,    snow_f,    solsw_m,  solswfdiff_m, sollw_m,       &
     264       rain_f,    snow_f,    bs_f, solsw_m,  solswfdiff_m, sollw_m,       &
    264265       gustiness,                                     &
    265        t,         q,        u,        v,             &
     266       t,         q,        qbs, u,        v,             &
    266267!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
    267268!!       t_x,       q_x,       t_w,      q_w,           &
    … …  
    275276       beta, &
    276277!>jyg
    277        alb_dir_m,    alb_dif_m,  zxsens,   zxevap,    &
     278       alb_dir_m,    alb_dif_m,  zxsens,   zxevap,  zxsnowerosion,  &
    278279       alb3_lic,  runoff,    snowhgt,   qsnow,     to_ice,    sissnow,  &
    279280       zxtsol,    zxfluxlat, zt2m,     qsat2m, zn2mout, &
    280        d_t,       d_q,       d_u,      d_v, d_t_diss, &
     281       d_t,       d_q,    d_qbs,    d_u,      d_v, d_t_diss, &
    281282!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
    282283       d_t_w,     d_q_w,                              &
    … …  
    301302       rh2m,      zxfluxu,  zxfluxv,               &
    302303       z0m, z0h,   agesno,  sollw,    solsw,         &
    303        d_ts,      evap,    fluxlat,  t2m,           &
     304       d_ts,      evap,    fluxlat,   t2m,           &
    304305       wfbils,    wfbilo, wfevap, wfrain, wfsnow,   &
    305306       flux_t,   flux_u, flux_v,                    &
    … …  
    307308!jyg<
    308309!!       zxfluxt,   zxfluxq,   q2m,      flux_q, tke,   &
    309        zxfluxt,   zxfluxq,   q2m,      flux_q, tke_x,  &
     310       zxfluxt,   zxfluxq, zxfluxqbs,   q2m, flux_q, flux_qbs, tke_x,  &
    310311!>jyg
    311312!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
    … …  
    410411         dser, dt_ds, zsig, zmea
    411412    use phys_output_var_mod, only: tkt, tks, taur, sss
     413    use blowing_snow_ini_mod, only : zeta_bs
    412414#ifdef CPP_XIOS
    413415    USE wxios, ONLY: missing_val
    … …  
    444446    REAL, DIMENSION(klon),        INTENT(IN)        :: rain_f  ! rain fall
    445447    REAL, DIMENSION(klon),        INTENT(IN)        :: snow_f  ! snow fall
     448    REAL, DIMENSION(klon),        INTENT(IN)        :: bs_f  ! blowing snow fall
    446449    REAL, DIMENSION(klon),        INTENT(IN)        :: solsw_m ! net shortwave radiation at mean surface
    447450    REAL, DIMENSION(klon),        INTENT(IN)        :: solswfdiff_m ! diffuse fraction fordownward shortwave radiation at mean surface
    … …  
    449452    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: t       ! temperature (K)
    450453    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: q       ! water vapour (kg/kg)
     454    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: qbs       ! blowing snow specific content (kg/kg)
    451455    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: u       ! u speed
    452456    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: v       ! v speed
    … …  
    521525                                                                  ! (=> positive sign upwards)
    522526    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxevap     ! water vapour flux at surface, positiv upwards
     527    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxsnowerosion     ! blowing snow flux at surface
    523528    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxtsol     ! temperature at surface, mean for each grid point
    524529!!! jyg le ???
    … …  
    537542    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_u        ! change in u speed
    538543    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_v        ! change in v speed
     544    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_qbs        ! change in blowing snow specific content
     545
    539546
    540547    REAL, INTENT(OUT):: zcoefh(:, :, :) ! (klon, klev, nbsrf + 1)
    … …  
    604611    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: sollw      ! net longwave radiation at surface
    605612    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: d_ts       ! change in temperature at surface
    606     REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: evap     ! evaporation at surface
     613    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: evap       ! evaporation at surface
    607614    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: fluxlat    ! latent flux
    608615    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: t2m        ! temperature at 2 meter height
    … …  
    631638    REAL, DIMENSION(klon, klev), INTENT(OUT)        :: zxfluxt    ! sensible heat flux, mean for each grid point
    632639    REAL, DIMENSION(klon, klev), INTENT(OUT)        :: zxfluxq    ! water vapour flux, mean for each grid point
     640    REAL, DIMENSION(klon, klev), INTENT(OUT)        :: zxfluxqbs    ! blowing snow flux, mean for each grid point
    633641    REAL, DIMENSION(klon, nbsrf),INTENT(OUT)        :: q2m        ! water vapour at 2 meter height
    634642    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_q     ! water vapour flux(latent flux) (kg/m**2/s)
     643    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_qbs   ! blowind snow vertical flux (kg/m**2
     644
    635645#ifdef ISO   
    636646    REAL, DIMENSION(ntraciso,klon),       INTENT(OUT)        :: dflux_xt    ! change of water vapour flux
    … …  
    683693    REAL, DIMENSION(klon)              :: yalb,yalb_vis
    684694!albedo SB <<<
    685     REAL, DIMENSION(klon)              :: yt1, yq1, yu1, yv1
     695    REAL, DIMENSION(klon)              :: yt1, yq1, yu1, yv1, yqbs1
    686696    REAL, DIMENSION(klon)              :: yqa
    687697    REAL, DIMENSION(klon)              :: ysnow, yqsurf, yagesno, yqsol
    688     REAL, DIMENSION(klon)              :: yrain_f, ysnow_f
     698    REAL, DIMENSION(klon)              :: yrain_f, ysnow_f, ybs_f
    689699#ifdef ISO
    690700    REAL, DIMENSION(ntraciso,klon)              :: yxt1
    … …  
    699709    REAL, DIMENSION(klon)              :: yrugoro
    700710    REAL, DIMENSION(klon)              :: yfluxlat
     711    REAL, DIMENSION(klon)              :: yfluxbs
    701712    REAL, DIMENSION(klon)              :: y_d_ts
    702713    REAL, DIMENSION(klon)              :: y_flux_t1, y_flux_q1
    … …  
    707718#endif
    708719    REAL, DIMENSION(klon)              :: y_flux_u1, y_flux_v1
     720    REAL, DIMENSION(klon)              :: y_flux_bs, y_flux0
    709721    REAL, DIMENSION(klon)              :: yt2m, yq2m, yu10m
    710722    INTEGER, DIMENSION(klon, nbsrf, 6) :: yn2mout, yn2mout_x, yn2mout_w
    … …  
    736748#endif
    737749    REAL, DIMENSION(klon)              :: AcoefU, AcoefV, BcoefU, BcoefV
     750    REAL, DIMENSION(klon)              :: AcoefQBS, BcoefQBS
    738751    REAL, DIMENSION(klon)              :: ypsref
    739752    REAL, DIMENSION(klon)              :: yevap, yevap_pot, ytsurf_new, yalb3_new
    … …  
    744757    REAL, DIMENSION(klon)              :: meansqT ! mean square deviation of subsurface temperatures
    745758    REAL, DIMENSION(klon)              :: alb_m  ! mean albedo for whole SW interval
    746     REAL, DIMENSION(klon,klev)         :: y_d_t, y_d_q, y_d_t_diss
     759    REAL, DIMENSION(klon,klev)         :: y_d_t, y_d_q, y_d_t_diss, y_d_qbs
    747760    REAL, DIMENSION(klon,klev)         :: y_d_u, y_d_v
    748     REAL, DIMENSION(klon,klev)         :: y_flux_t, y_flux_q
     761    REAL, DIMENSION(klon,klev)         :: y_flux_t, y_flux_q, y_flux_qbs
    749762    REAL, DIMENSION(klon,klev)         :: y_flux_u, y_flux_v
    750     REAL, DIMENSION(klon,klev)         :: ycoefh, ycoefm,ycoefq
     763    REAL, DIMENSION(klon,klev)         :: ycoefh, ycoefm,ycoefq, ycoefqbs
    751764    REAL, DIMENSION(klon)              :: ycdragh, ycdragq, ycdragm
    752765    REAL, DIMENSION(klon,klev)         :: yu, yv
    753     REAL, DIMENSION(klon,klev)         :: yt, yq
     766    REAL, DIMENSION(klon,klev)         :: yt, yq, yqbs
    754767#ifdef ISO
    755768    REAL, DIMENSION(ntraciso,klon)              :: yxtevap
    … …  
    819832    REAL, DIMENSION(klon,klev)         :: CcoefH, CcoefQ, DcoefH, DcoefQ
    820833    REAL, DIMENSION(klon,klev)         :: CcoefU, CcoefV, DcoefU, DcoefV
     834    REAL, DIMENSION(klon,klev)         :: CcoefQBS, DcoefQBS
    821835    REAL, DIMENSION(klon,klev)         :: CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x
    822836    REAL, DIMENSION(klon,klev)         :: CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w
    … …  
    824838    REAL, DIMENSION(klon,klev)         :: CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w
    825839    REAL, DIMENSION(klon,klev)         :: Kcoef_hq, Kcoef_m, gama_h, gama_q
     840    REAL, DIMENSION(klon,klev)         :: gama_qbs, Kcoef_qbs
    826841    REAL, DIMENSION(klon,klev)         :: Kcoef_hq_x, Kcoef_m_x, gama_h_x, gama_q_x
    827842    REAL, DIMENSION(klon,klev)         :: Kcoef_hq_w, Kcoef_m_w, gama_h_w, gama_q_w
    … …  
    10161031    REAL, DIMENSION(klon,nbsrf)        :: zx_t1
    10171032    REAL, DIMENSION(klon, nbsrf)       :: alb          ! mean albedo for whole SW interval
     1033    REAL, DIMENSION(klon,nbsrf)        :: snowerosion   
    10181034    REAL, DIMENSION(klon)              :: ylwdown      ! jg : temporary (ysollwdown)
    10191035    REAL, DIMENSION(klon)              :: ygustiness      ! jg : temporary (ysollwdown)
    … …  
    11891205  alb_dir_m=0. ; alb_dif_m=0. ; alb3_lic(:)=0.
    11901206!albedo SB <<<
    1191  zxsens(:)=0. ; zxevap(:)=0. ; zxtsol(:)=0.
     1207 zxsens(:)=0. ; zxevap(:)=0. ; zxtsol(:)=0. ; zxsnowerosion(:)=0.
    11921208 d_t_w(:,:)=0. ; d_q_w(:,:)=0. ; d_t_x(:,:)=0. ; d_q_x(:,:)=0.
    11931209 zxfluxlat(:)=0.
    11941210 zt2m(:)=0. ; zq2m(:)=0. ; qsat2m(:)=0. ; rh2m(:)=0.
    11951211 zn2mout(:,:)=0 ;
    1196  d_t(:,:)=0. ; d_t_diss(:,:)=0. ; d_q(:,:)=0. ; d_u(:,:)=0. ; d_v(:,:)=0.
     1212 d_t(:,:)=0. ; d_t_diss(:,:)=0. ; d_q(:,:)=0. ; d_qbs(:,:)=0. ; d_u(:,:)=0. ; d_v(:,:)=0.
    11971213 zcoefh(:,:,:)=0. ; zcoefm(:,:,:)=0.
    11981214 zxsens_x(:)=0. ; zxsens_w(:)=0. ; zxfluxlat_x(:)=0. ; zxfluxlat_w(:)=0.
    … …  
    12141230 d_ts(:,:)=0.
    12151231 evap(:,:)=0.
     1232 snowerosion(:,:)=0.
    12161233 fluxlat(:,:)=0.
    12171234 wfbils(:,:)=0. ; wfbilo(:,:)=0.
    12181235 wfevap(:,:)=0. ; wfrain(:,:)=0. ; wfsnow(:,:)=0.
    12191236 flux_t(:,:,:)=0. ; flux_q(:,:,:)=0. ; flux_u(:,:,:)=0. ; flux_v(:,:,:)=0.
     1237 flux_qbs(:,:,:)=0.
    12201238 dflux_t(:)=0. ; dflux_q(:)=0.
    12211239 zxsnow(:)=0.
    1222  zxfluxt(:,:)=0. ; zxfluxq(:,:)=0.
     1240 zxfluxt(:,:)=0. ; zxfluxq(:,:)=0.; zxfluxqbs(:,:)=0.
    12231241 qsnow(:)=0. ; snowhgt(:)=0. ; to_ice(:)=0. ; sissnow(:)=0.
    12241242 runoff(:)=0.
    … …  
    12661284    yqsurf = 0.0  ; yalb = 0.0 ; yalb_vis = 0.0
    12671285!albedo SB <<<
    1268     yrain_f = 0.0 ; ysnow_f = 0.0    ; yfder = 0.0     ; ysolsw = 0.0
     1286    yrain_f = 0.0 ; ysnow_f = 0.0  ; ybs_f=0.0  ; yfder = 0.0     ; ysolsw = 0.0
    12691287    ysollw = 0.0  ; yz0m = 0.0 ; yz0h = 0.0    ; yu1 = 0.0   
    1270     yv1 = 0.0     ; ypaprs = 0.0     ; ypplay = 0.0
     1288    yv1 = 0.0     ; ypaprs = 0.0     ; ypplay = 0.0     ; yqbs1 = 0.0
    12711289    ydelp = 0.0   ; yu = 0.0         ; yv = 0.0        ; yt = 0.0         
    12721290    yq = 0.0      ; y_dflux_t = 0.0  ; y_dflux_q = 0.0
     1291    yqbs(:,:)=0.0 
    12731292    yrugoro = 0.0 ; ywindsp = 0.0   
    12741293!!    d_ts = 0.0    ; yfluxlat=0.0     ; flux_t = 0.0    ; flux_q = 0.0     
    1275     yfluxlat=0.0
     1294    yfluxlat=0.0 ; y_flux0(:)=0.0
    12761295!!    flux_u = 0.0  ; flux_v = 0.0     ; d_t = 0.0       ; d_q = 0.0     
    12771296!!    d_t_diss= 0.0 ;d_u = 0.0     ; d_v = 0.0
    … …  
    12881307    ycldt = 0.0      ; yrmu0 = 0.0
    12891308    ! Martin
     1309    y_d_qbs(:,:)=0.0
    12901310
    12911311!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
    … …  
    16291649          yrain_f(j) = rain_f(i)
    16301650          ysnow_f(j) = snow_f(i)
     1651          ybs_f(j)   = bs_f(i)
    16311652          yagesno(j) = agesno(i,nsrf)
    16321653          yfder(j)   = fder(i)
    … …  
    16401661          yu1(j)     = u(i,1)
    16411662          yv1(j)     = v(i,1)
     1663          yqbs1(j)   = qbs(i,1)
    16421664          ypaprs(j,klev+1) = paprs(i,klev+1)
    16431665!jyg<
    … …  
    16531675!!! nrlmd le 13/06/2011
    16541676          y_delta_tsurf(j)=delta_tsurf(i,nsrf)
     1677          yfluxbs(j)=0.0
     1678          y_flux_bs(j) = 0.0
    16551679!!!
    16561680#ifdef ISO
    … …  
    17211745             yt(j,k) = t(i,k)
    17221746             yq(j,k) = q(i,k)
     1747             yqbs(j,k)=qbs(i,k)
    17231748#ifdef ISO
    17241749             do ixt=1,ntraciso   
    … …  
    19511976      print *,' args coef_diff_turb: ycdragh ', ycdragh
    19521977      print *,' args coef_diff_turb: ytke ', ytke
    1953 
    19541978       ENDIF
    19551979
    … …  
    19601984
    19611985        ELSE
    1962 
    19631986
    19641987        CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
    … …  
    19812004
    19822005        IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh ',ycoefh
    1983 !
     2006
     2007
    19842008       ELSE  !(iflag_split .eq.0)
     2009
    19852010      IF (prt_level >=10) THEN
    19862011      print *,' args coef_diff_turb: yu_x ',  yu_x 
    … …  
    20202045       ENDIF
    20212046
    2022        ENDIF ! iflag_pbl >= 50
     2047        ENDIF ! iflag_pbl >= 50
    20232048
    20242049        IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_x ',ycoefh_x
    … …  
    20342059      print *,' args coef_diff_turb: ycdragh_w ', ycdragh_w
    20352060      print *,' args coef_diff_turb: ytke_w ', ytke_w
    2036        ENDIF
     2061      ENDIF
    20372062
    20382063        IF (iflag_pbl>=50) THEN
    … …  
    20622087
    20632088        IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_w ',ycoefh_w
    2064 !
     2089
    20652090!!!jyg le 10/04/2013
    20662091!!   En attendant de traiter le transport des traceurs dans les poches froides, formule
    … …  
    20722097        ENDDO
    20732098      ENDDO
    2074 !!!
    20752099       ENDIF  ! (iflag_split .eq.0)
    2076 !!!
    20772100       
    20782101!****************************************************************************************
    … …  
    21712194       ENDIF  ! (iflag_split .eq.0)
    21722195!!!
     2196
     2197 ! For blowing snow:
     2198    IF (ok_bs) THEN
     2199     ! following Bintanja et al 2000, part II
     2200     ! we assume that the eddy diffuvisity coefficient for
     2201     ! suspended particles is larger than Km by a factor zeta_bs
     2202     ! which is equal to 3 by default
     2203     ycoefqbs=ycoefm*zeta_bs
     2204     CALL climb_qbs_down(knon, ycoefqbs, ypaprs, ypplay, &
     2205     ydelp, yt, yqbs, dtime, &
     2206     CcoefQBS, DcoefQBS, &
     2207     Kcoef_qbs, gama_qbs, &
     2208     AcoefQBS, BcoefQBS)
     2209    ENDIF
     2210
    21732211
    21742212!****************************************************************************************
    … …  
    23562394               debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, &
    23572395!!jyg               yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
    2358                yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&
     2396               yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt1, yq1,&
    23592397               AcoefH, AcoefQ, BcoefH, BcoefQ, &
    23602398               AcoefU, AcoefV, BcoefU, BcoefV, &
    … …  
    23622400               ylwdown, yq2m, yt2m, &
    23632401               ysnow, yqsol, yagesno, ytsoil, &
    2364                yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
     2402               yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,yfluxbs,&
    23652403               yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &
    23662404               y_flux_u1, y_flux_v1, &
    … …  
    24312469                  yrmu0, ylwdown, yalb, zgeo1, &
    24322470                  ysolsw, ysollw, yts, ypplay(:,1), &
    2433                   !!jyg               ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
    2434                   ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&
     2471                  ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt1, yq1,&
    24352472                  AcoefH, AcoefQ, BcoefH, BcoefQ, &
    24362473                  AcoefU, AcoefV, BcoefU, BcoefV, &
    24372474                  ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
    2438                   ysnow, yqsurf, yqsol, yagesno, &
     2475                  ysnow, yqsurf, yqsol,yqbs1, yagesno, &
    24392476                  ytsoil, yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap,yfluxsens,yfluxlat, &
    2440                   ytsurf_new, y_dflux_t, y_dflux_q, &
     2477                  yfluxbs, ytsurf_new, y_dflux_t, y_dflux_q, &
    24412478                  yzmea, yzsig, ycldt, &
    24422479                  ysnowhgt, yqsnow, ytoice, ysissnow, &
    … …  
    24952532               itap, dtime, jour, knon, ni, &
    24962533!!jyg               ypplay(:,1), zgeo1/RG, ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
    2497                ypplay(:,1), zgeo1(1:knon)/RG, ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&    ! ym missing init
     2534               ypplay(:,1), zgeo1(1:knon)/RG, ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt(:,1), yq(:,1),&    ! ym missing init
    24982535               AcoefH, AcoefQ, BcoefH, BcoefQ, &
    24992536               AcoefU, AcoefV, BcoefU, BcoefV, &
    … …  
    26582695          ENDDO
    26592696        ENDIF ! (ok_flux_surf)
     2697
     2698        ! flux of blowing snow at the first level
     2699        IF (ok_bs) THEN
     2700        DO j=1,knon
     2701        y_flux_bs(j)=yfluxbs(j)
     2702        ENDDO
     2703        ENDIF
    26602704!
    26612705! ------------------------------------------------------------------------------
    … …  
    29853029!
    29863030       ENDIF  ! (iflag_split .eq.0)
     3031
     3032       IF (ok_bs) THEN
     3033            CALL climb_qbs_up(knon, dtime, yqbs, &
     3034            y_flux_bs, ypaprs, ypplay, &
     3035            AcoefQBS, BcoefQBS, &
     3036            CcoefQBS, DcoefQBS, &
     3037            Kcoef_qbs, gama_qbs, &
     3038            y_flux_qbs(:,:), y_d_qbs(:,:))
     3039       ENDIF
     3040
    29873041!!!
    29883042!!
    … …  
    31513205!!!
    31523206
    3153 !      print*,'Dans pbl OK1'
    3154 
    3155 !jyg<
    3156 !!       evap(:,nsrf) = - flux_q(:,1,nsrf)
    3157 !>jyg
     3207       ! tendencies of blowing snow
     3208       IF (ok_bs) THEN
     3209           DO k = 1, klev   
     3210            DO j = 1, knon
     3211                i = ni(j)
     3212                y_d_qbs(j,k)=y_d_qbs(j,k) * ypct(j)
     3213                flux_qbs(i,k,nsrf) = y_flux_qbs(j,k)
     3214            ENDDO
     3215          ENDDO
     3216       ENDIF
     3217
     3218
    31583219       DO j = 1, knon
    31593220          i = ni(j)
    31603221          evap(i,nsrf) = - flux_q(i,1,nsrf)                  !jyg
     3222          if (ok_bs) then ; snowerosion(i,nsrf)=flux_qbs(i,1,nsrf); endif
    31613223          beta(i,nsrf) = ybeta(j)                             !jyg
    31623224          d_ts(i,nsrf) = y_d_ts(j)
    … …  
    33863448          ENDDO
    33873449       ENDDO
     3450
     3451
     3452       IF (ok_bs) THEN
     3453         DO k = 1, klev
     3454         DO j = 1, knon
     3455         i = ni(j)
     3456         d_qbs(i,k) = d_qbs(i,k) + y_d_qbs(j,k)
     3457         ENDDO
     3458         ENDDO
     3459        ENDIF
     3460
    33883461
    33893462#ifdef ISO
    … …  
    38883961       fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
    38893962    ENDDO
     3963
     3964    ! if blowing snow
     3965    if (ok_bs) then 
     3966       DO nsrf = 1, nbsrf
     3967       DO k = 1, klev
     3968       DO i = 1, klon
     3969         zxfluxqbs(i,k) = zxfluxqbs(i,k) + flux_qbs(i,k,nsrf) * pctsrf(i,nsrf)
     3970       ENDDO
     3971       ENDDO
     3972       ENDDO
     3973
     3974       DO i = 1, klon
     3975        zxsnowerosion(i)     = zxfluxqbs(i,1) ! blowings snow flux at the surface
     3976       END DO
     3977    endif
     3978
     3979   
     3980   
    38903981#ifdef ISO
    38913982    DO i = 1, klon

  • LMDZ6/branches/blowing_snow/libf/phylmdiso/phyredem.F90

    r4389 r4506  
    1515                                ftsol, beta_aridity, delta_tsurf, falb_dir,  &
    1616                                falb_dif, qsol, fevap, radsol, solsw, sollw, &
    17                                 sollwdown, rain_fall, snow_fall, z0m, z0h,  &
     17                                sollwdown, rain_fall, snow_fall, bs_fall, z0m, z0h, &
    1818                                agesno, zmea, zstd, zsig, zgam, zthe, zpic,  &
    1919                                zval, rugoro, t_ancien, q_ancien,            &
    20                                 prw_ancien, prlw_ancien, prsw_ancien,        &
    21                                 ql_ancien, qs_ancien,  u_ancien,            &
     20                                prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien,      &
     21                                ql_ancien, qs_ancien, qbs_ancien, rneb_ancien, u_ancien, &
    2222                                v_ancien, clwcon, rnebcon, ratqs, pbl_tke,   &
    2323                                wake_delta_pbl_tke, zmax0, f0, sig1, w01,    &
    … …  
    260260
    261261    CALL put_field(pass,"QSANCIEN", "QSANCIEN", qs_ancien)
     262
     263    IF (ok_bs) THEN
     264       CALL put_field(pass,"bs_f", "precipitation neige soufflee", bs_fall)
     265       CALL put_field(pass,"QBSANCIEN", "QBSANCIEN", qbs_ancien)
     266       CALL put_field(pass,"PRBSWANCIEN", "PRBSWANCIEN", prbsw_ancien)
     267    ENDIF
     268
     269    CALL put_field(pass,"RNEBANCIEN", "RNEBANCIEN", rneb_ancien)
    262270
    263271    CALL put_field(pass,"PRWANCIEN", "PRWANCIEN", prw_ancien)

  • LMDZ6/branches/blowing_snow/libf/phylmdiso/phys_local_var_mod.F90

    r4143 r4506  
    1414      REAL, SAVE, ALLOCATABLE :: ql_seri(:,:),qs_seri(:,:)
    1515      !$OMP THREADPRIVATE(ql_seri,qs_seri)
     16      REAL, SAVE, ALLOCATABLE :: qbs_seri(:,:)
     17      !$OMP THREADPRIVATE(qbs_seri)
    1618      REAL, SAVE, ALLOCATABLE :: u_seri(:,:), v_seri(:,:)
    1719      !$OMP THREADPRIVATE(u_seri, v_seri)
    … …  
    2628      REAL, SAVE, ALLOCATABLE :: d_t_dyn(:,:), d_q_dyn(:,:)
    2729      !$OMP THREADPRIVATE(d_t_dyn, d_q_dyn)
    28       REAL, SAVE, ALLOCATABLE :: d_ql_dyn(:,:), d_qs_dyn(:,:)
    29       !$OMP THREADPRIVATE(d_ql_dyn, d_qs_dyn)
    30       REAL, SAVE, ALLOCATABLE :: d_q_dyn2d(:), d_ql_dyn2d(:), d_qs_dyn2d(:)
    31       !$OMP THREADPRIVATE(d_q_dyn2d, d_ql_dyn2d, d_qs_dyn2d)
     30      REAL, SAVE, ALLOCATABLE :: d_ql_dyn(:,:), d_qs_dyn(:,:), d_qbs_dyn(:,:)
     31      !$OMP THREADPRIVATE(d_ql_dyn, d_qs_dyn, d_qbs_dyn)
     32      REAL, SAVE, ALLOCATABLE :: d_q_dyn2d(:), d_ql_dyn2d(:), d_qs_dyn2d(:), d_qbs_dyn2d(:)
     33      !$OMP THREADPRIVATE(d_q_dyn2d, d_ql_dyn2d, d_qs_dyn2d, d_qbs_dyn2d)
    3234      REAL, SAVE, ALLOCATABLE :: d_u_dyn(:,:), d_v_dyn(:,:)
    3335      !$OMP THREADPRIVATE(d_u_dyn, d_v_dyn)
    … …  
    6971      REAL, SAVE, ALLOCATABLE :: d_u_oli(:,:), d_v_oli(:,:)
    7072      !$OMP THREADPRIVATE(d_u_oli, d_v_oli)
    71       REAL, SAVE, ALLOCATABLE :: d_t_vdf(:,:), d_q_vdf(:,:), d_t_diss(:,:)
    72       !$OMP THREADPRIVATE( d_t_vdf, d_q_vdf,d_t_diss)
     73      REAL, SAVE, ALLOCATABLE :: d_t_vdf(:,:), d_q_vdf(:,:), d_qbs_vdf(:,:), d_t_diss(:,:)
     74      !$OMP THREADPRIVATE( d_t_vdf, d_q_vdf, d_qbs_vdf, d_t_diss)
    7375      REAL, SAVE, ALLOCATABLE :: d_u_vdf(:,:), d_v_vdf(:,:)
    7476      !$OMP THREADPRIVATE(d_u_vdf, d_v_vdf)
    … …  
    7880      REAL, SAVE, ALLOCATABLE :: d_t_vdf_x(:,:), d_q_vdf_x(:,:)
    7981      !$OMP THREADPRIVATE( d_t_vdf_x, d_q_vdf_x)
     82      REAL, SAVE, ALLOCATABLE :: d_t_bs(:,:), d_q_bs(:,:), d_qbs_bs(:,:)
     83      !$OMP THREADPRIVATE( d_t_bs,d_q_bs, d_qbs_bs)
    8084!>nrlmd+jyg
    8185      REAL, SAVE, ALLOCATABLE :: d_t_oro(:,:)
    … …  
    361365      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: JrNt
    362366!$OMP THREADPRIVATE(JrNt)
    363       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: dthmin, evap, fder, plcl, plfc, prw, prlw, prsw
    364 !$OMP THREADPRIVATE(dthmin, evap, fder, plcl, plfc, prw, prlw, prsw)
     367      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: dthmin, evap, snowerosion, fder, plcl, plfc, prw, prlw, prsw, prbsw
     368!$OMP THREADPRIVATE(dthmin, evap, snowerosion, fder, plcl, plfc, prw, prlw, prsw, prbsw)
    365369      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zustar, zu10m, zv10m, rh2m
    366370!$OMP THREADPRIVATE(zustar, zu10m, zv10m, rh2m)
    … …  
    379383      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: tpot, tpote, ue, uq, uwat, ve, vq, vwat, zxffonte
    380384!$OMP THREADPRIVATE(tpot, tpote, ue, uq, uwat, ve, vq, vwat, zxffonte)
     385      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zxustartlic, zxrhoslic
     386!$OMP THREADPRIVATE(zxustartlic, zxrhoslic)
    381387      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zxfqcalving
    382388!$OMP THREADPRIVATE(zxfqcalving)
    … …  
    567573      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: zx_rh, zx_rhl, zx_rhi
    568574!$OMP THREADPRIVATE(zx_rh, zx_rhl, zx_rhi)
    569       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: prfl, psfl, fraca
    570 !$OMP THREADPRIVATE(prfl, psfl, fraca)
     575      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: prfl, psfl, fraca, bsfl
     576!$OMP THREADPRIVATE(prfl, psfl, fraca, bsfl)
    571577      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: Vprecip, zw2
    572578!$OMP THREADPRIVATE(Vprecip, zw2)
    … …  
    734740
    735741IMPLICIT NONE
    736       ALLOCATE(t_seri(klon,klev),q_seri(klon,klev),ql_seri(klon,klev),qs_seri(klon,klev))
     742      ALLOCATE(t_seri(klon,klev),q_seri(klon,klev),ql_seri(klon,klev),qs_seri(klon,klev), qbs_seri(klon,klev))
    737743      ALLOCATE(u_seri(klon,klev),v_seri(klon,klev))
    738744      ALLOCATE(l_mixmin(klon,klev+1,nbsrf),l_mix(klon,klev+1,nbsrf),tke_dissip(klon,klev+1,nbsrf),wprime(klon,klev+1,nbsrf))
    … …  
    741747      ALLOCATE(tr_seri(klon,klev,nbtr))
    742748      ALLOCATE(d_t_dyn(klon,klev),d_q_dyn(klon,klev))
    743       ALLOCATE(d_ql_dyn(klon,klev),d_qs_dyn(klon,klev))
    744       ALLOCATE(d_q_dyn2d(klon),d_ql_dyn2d(klon),d_qs_dyn2d(klon))
     749      ALLOCATE(d_ql_dyn(klon,klev),d_qs_dyn(klon,klev), d_qbs_dyn(klon,klev))
     750      ALLOCATE(d_q_dyn2d(klon),d_ql_dyn2d(klon),d_qs_dyn2d(klon), d_qbs_dyn2d(klon))
    745751      ALLOCATE(d_u_dyn(klon,klev),d_v_dyn(klon,klev))
    746752      ALLOCATE(d_tr_dyn(klon,klev,nbtr))                   !RomP
    … …  
    765771      ALLOCATE(plul_st(klon),plul_th(klon))
    766772      ALLOCATE(d_t_vdf(klon,klev),d_q_vdf(klon,klev),d_t_diss(klon,klev))
    767 
     773      ALLOCATE (d_qbs_vdf(klon,klev))
     774      ALLOCATE(d_t_bs(klon,klev),d_q_bs(klon,klev),d_qbs_bs(klon,klev))
    768775      ALLOCATE(d_t_vdf_w(klon,klev),d_q_vdf_w(klon,klev))
    769776      ALLOCATE(d_t_vdf_x(klon,klev),d_q_vdf_x(klon,klev))
    … …  
    925932      ALLOCATE(cldm(klon), cldq(klon), cldt(klon), qsat2m(klon))
    926933      ALLOCATE(JrNt(klon))
    927       ALLOCATE(dthmin(klon), evap(klon), fder(klon), plcl(klon), plfc(klon))
    928       ALLOCATE(prw(klon), prlw(klon), prsw(klon), zustar(klon), zu10m(klon), zv10m(klon), rh2m(klon))
     934      ALLOCATE(dthmin(klon), evap(klon), snowerosion(klon), fder(klon), plcl(klon), plfc(klon))
     935      ALLOCATE(prw(klon), prlw(klon), prsw(klon), prbsw(klon), zustar(klon), zu10m(klon), zv10m(klon), rh2m(klon))
    929936      ALLOCATE(s_lcl(klon))
    930937      ALLOCATE(s_pblh(klon), s_pblt(klon), s_therm(klon))
    … …  
    10551062      ALLOCATE(prfl(klon, klev+1))
    10561063      ALLOCATE(psfl(klon, klev+1), fraca(klon, klev+1), Vprecip(klon, klev+1))
     1064      ALLOCATE(bsfl(klon,klev+1))
    10571065      ALLOCATE(zw2(klon, klev+1))
    10581066
    … …  
    11411149USE indice_sol_mod
    11421150IMPLICIT NONE
    1143       DEALLOCATE(t_seri,q_seri,ql_seri,qs_seri)
     1151      DEALLOCATE(t_seri,q_seri,ql_seri,qs_seri, qbs_seri)
    11441152      DEALLOCATE(u_seri,v_seri)
    11451153      DEALLOCATE(l_mixmin,l_mix, tke_dissip,wprime)
    … …  
    11471155      DEALLOCATE(tr_seri)
    11481156      DEALLOCATE(d_t_dyn,d_q_dyn)
    1149       DEALLOCATE(d_ql_dyn,d_qs_dyn)
    1150       DEALLOCATE(d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d)
     1157      DEALLOCATE(d_ql_dyn,d_qs_dyn, d_qbs_dyn)
     1158      DEALLOCATE(d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d, d_qbs_dyn2d)
    11511159      DEALLOCATE(d_u_dyn,d_v_dyn)
    11521160      DEALLOCATE(d_tr_dyn)                      !RomP
    … …  
    11711179      DEALLOCATE(plul_st,plul_th)
    11721180      DEALLOCATE(d_t_vdf,d_q_vdf,d_t_diss)
     1181      DEALLOCATE(d_qbs_vdf)
     1182      DEALLOCATE(d_t_bs,d_q_bs,d_qbs_bs)
    11731183#ifdef ISO
    11741184      deallocate(xt_seri,xtl_seri,xts_seri)
    … …  
    13081318      DEALLOCATE(cldm, cldq, cldt, qsat2m)
    13091319      DEALLOCATE(JrNt)
    1310       DEALLOCATE(dthmin, evap, fder, plcl, plfc)
    1311       DEALLOCATE(prw, prlw, prsw, zustar, zu10m, zv10m, rh2m, s_lcl)
     1320      DEALLOCATE(dthmin, evap, snowerosion, fder, plcl, plfc)
     1321      DEALLOCATE(prw, prlw, prsw, prbsw, zustar, zu10m, zv10m, rh2m, s_lcl)
    13121322      DEALLOCATE(s_pblh, s_pblt, s_therm)
    13131323!
    … …  
    13221332      DEALLOCATE(zxfqcalving, zxfluxlat)
    13231333      DEALLOCATE(zxrunofflic)
     1334      DEALLOCATE(zxustartlic, zxrhoslic)
    13241335      DEALLOCATE(zxtsol, snow_lsc, zxfqfonte, zxqsurf)
    13251336      DEALLOCATE(rain_lsc)
    … …  
    14231434
    14241435
    1425       DEALLOCATE(prfl, psfl, fraca, Vprecip)
     1436      DEALLOCATE(prfl, psfl, bsfl, fraca, Vprecip)
    14261437      DEALLOCATE(zw2)
    14271438

  • LMDZ6/branches/blowing_snow/libf/phylmdiso/phys_output_ctrlout_mod.F90

    r4065 r4506  
    378378  TYPE(ctrl_out), SAVE :: o_snow = ctrl_out((/ 1, 1, 10, 10, 5, 10, 11, 11, 11, 11/), &
    379379    'snow', 'Snow fall', 'kg/(s*m2)', (/ ('', i=1, 10) /))
     380  TYPE(ctrl_out), SAVE :: o_bsfall = ctrl_out((/ 10, 10, 10, 10, 5, 10, 11, 11, 11, 11/), &
     381    'bsfall', 'Blowing Snow fall', 'kg/(s*m2)', (/ ('', i=1, 10) /))
    380382  TYPE(ctrl_out), SAVE :: o_evap = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11, 11/), &
    381383    'evap', 'Evaporat', 'kg/(s*m2)', (/ ('', i=1, 10) /))
    382 
     384  TYPE(ctrl_out), SAVE :: o_snowerosion = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     385   'snowerosion', 'blowing snow flux', 'kg/(s*m2)', (/ ('', i=1, 10) /))
     386  TYPE(ctrl_out), SAVE :: o_ustart_lic = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     387    'ustart_lic', 'threshold velocity', 'm/s', (/ ('', i=1, 10) /))
     388  TYPE(ctrl_out), SAVE :: o_rhosnow_lic = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     389    'rhosnow_lic', 'snow density lic', 'kg/m3', (/ ('', i=1, 10) /))
    383390  TYPE(ctrl_out), SAVE :: o_sens_prec_liq_oce = ctrl_out((/ 5, 5, 10, 10, 5, 10, 11, 11, 11, 11/), &
    384391    'sens_rain_oce', 'Sensible heat flux of liquid prec. over ocean', 'W/m2', (/ ('', i=1, 10) /))
    … …  
    743750  TYPE(ctrl_out), SAVE :: o_prsw = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11, 11/), &
    744751    'prsw', 'Precipitable solid water', 'kg/m2', (/ ('', i=1, 10) /))
     752  TYPE(ctrl_out), SAVE :: o_prbsw = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     753    'prbsw', 'Precipitable blowing snow', 'kg/m2', (/ ('', i=1, 10) /))
    745754  TYPE(ctrl_out), SAVE :: o_s_pblh = ctrl_out((/ 1, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
    746755    's_pblh', 'Boundary Layer Height', 'm', (/ ('', i=1, 10) /))
    … …  
    14131422    'ec550aer', 'Extinction at 550nm', 'm^-1', (/ ('', i=1, 10) /))
    14141423  TYPE(ctrl_out), SAVE :: o_lwcon = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11, 11, 11/), &
    1415     'lwcon', 'Cloud liquid water content', 'kg/kg', (/ ('', i=1, 10) /))
     1424    'lwcon', 'Cloud liquid water content seen by radiation', 'kg/kg', (/ ('', i=1, 10) /))
    14161425  TYPE(ctrl_out), SAVE :: o_iwcon = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11, 11, 11/), &
    1417     'iwcon', 'Cloud ice water content', 'kg/kg', (/ ('', i=1, 10) /))
     1426    'iwcon', 'Cloud ice water content seen by radiation', 'kg/kg', (/ ('', i=1, 10) /))
    14181427  TYPE(ctrl_out), SAVE :: o_temp = ctrl_out((/ 2, 3, 4, 10, 10, 10, 11, 11, 11, 11/), &
    14191428    'temp', 'Air temperature', 'K', (/ ('', i=1, 10) /))
    … …  
    14321441  TYPE(ctrl_out), SAVE :: o_ocond = ctrl_out((/ 2, 3, 4, 10, 10, 10, 11, 11, 11, 11/), &
    14331442    'ocond', 'Condensed water', 'kg/kg', (/ ('', i=1, 10) /))
     1443  TYPE(ctrl_out), SAVE :: o_qbs = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     1444    'qbs', 'Specific content of blowing snow', 'kg/kg', (/ ('', i=1, 10) /))
    14341445  TYPE(ctrl_out), SAVE :: o_wvapp = ctrl_out((/ 2, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
    14351446    'wvapp', '', '', (/ ('', i=1, 10) /))
    … …  
    14941505  TYPE(ctrl_out), SAVE :: o_dqsphy2d = ctrl_out((/ 2, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
    14951506    'dqsphy2d', 'Physics dQS', '(kg/m2)/s', (/ ('', i=1, 10) /))
     1507  TYPE(ctrl_out), SAVE :: o_dqbsphy = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     1508    'dqbsphy', 'Physics dQBS', '(kg/kg)/s', (/ ('', i=1, 10) /))
     1509  TYPE(ctrl_out), SAVE :: o_dqbsphy2d = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     1510    'dqbsphy2d', 'Physics dQBS', '(kg/m2)/s', (/ ('', i=1, 10) /))
    14961511  TYPE(ctrl_out), SAVE :: o_pr_con_l = ctrl_out((/ 2, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
    14971512    'pr_con_l', 'Convective precipitation lic', ' ', (/ ('', i=1, 10) /))
    … …  
    15021517  TYPE(ctrl_out), SAVE :: o_pr_lsc_i = ctrl_out((/ 2, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
    15031518    'pr_lsc_i', 'Large scale precipitation ice', ' ', (/ ('', i=1, 10) /))
     1519  TYPE(ctrl_out), SAVE :: o_pr_bs = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     1520    'pr_bs', 'profile of blowing snow flux', ' ', (/ ('', i=1, 10) /))
    15041521  TYPE(ctrl_out), SAVE :: o_re = ctrl_out((/ 5, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
    15051522    're', 'Cloud droplet effective radius', 'um', (/ ('', i=1, 10) /))
    … …  
    15991616  TYPE(ctrl_out), SAVE :: o_dqsdyn2d = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
    16001617    'dqsdyn2d', 'Dynamics dQS', '(kg/m2)/s', (/ ('', i=1, 10) /))
     1618  TYPE(ctrl_out), SAVE :: o_dqbsdyn = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     1619    'dqbsdyn', 'Dynamics dQBS', '(kg/kg)/s', (/ ('', i=1, 10) /))
     1620  TYPE(ctrl_out), SAVE :: o_dqbsdyn2d = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     1621    'dqbsdyn2d', 'Dynamics dQBS', '(kg/m2)/s', (/ ('', i=1, 10) /))
    16011622  TYPE(ctrl_out), SAVE :: o_dudyn = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
    16021623    'dudyn', 'Dynamics dU', 'm/s2', (/ ('', i=1, 10) /))
    … …  
    16731694  TYPE(ctrl_out), SAVE :: o_dqeva2d = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
    16741695    'dqeva2d', 'Reevaporation dQ', '(kg/m2)/s', (/ ('', i=1, 10) /))
     1696  TYPE(ctrl_out), SAVE :: o_dqbsvdf = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     1697    'dqbsvdf', 'Boundary-layer dQBS', '(kg/kg)/s', (/ ('', i=1, 10) /))
     1698  TYPE(ctrl_out), SAVE :: o_dqbsbs = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     1699    'dqbsbs', 'Blowing snow dQBS', '(kg/kg)/s', (/ ('', i=1, 10) /))
     1700  TYPE(ctrl_out), SAVE :: o_dtbs = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     1701    'dtbs', 'Blowing snow dT', '(K)/s', (/ ('', i=1, 10) /))
     1702  TYPE(ctrl_out), SAVE :: o_dqbs = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     1703    'dqbs', 'Blowing snow dQ', '(kg/kg)/s', (/ ('', i=1, 10) /))
    16751704
    16761705!!!!!!!!!!!!!!!! Specifique thermiques

  • LMDZ6/branches/blowing_snow/libf/phylmdiso/phys_output_var_mod.F90

    r4374 r4506  
    2929  REAL, SAVE, ALLOCATABLE :: d_qw_col(:)      ! watter vapour mass budget for each column (kg/m2/s)
    3030  REAL, SAVE, ALLOCATABLE :: d_ql_col(:)      ! liquid watter mass budget for each column (kg/m2/s)
    31   REAL, SAVE, ALLOCATABLE :: d_qs_col(:)      ! solid watter mass budget for each column (kg/m2/s)
     31  REAL, SAVE, ALLOCATABLE :: d_qs_col(:)      ! cloud ice mass budget for each column (kg/m2/s)
     32  REAL, SAVE, ALLOCATABLE :: d_qbs_col(:)     ! blowing snow mass budget for each column (kg/m2/s)
    3233  REAL, SAVE, ALLOCATABLE :: d_qt_col(:)      ! total watter mass budget for each column (kg/m2/s)
    3334  REAL, SAVE, ALLOCATABLE :: d_ek_col(:)      ! kinetic energy budget for each column (W/m2)
    … …  
    3536  REAL, SAVE, ALLOCATABLE :: d_h_qw_col(:)    ! enthalpy budget of watter vapour for each column (W/m2)
    3637  REAL, SAVE, ALLOCATABLE :: d_h_ql_col(:)    ! enthalpy budget of liquid watter for each column (W/m2)
    37   REAL, SAVE, ALLOCATABLE :: d_h_qs_col(:)    ! enthalpy budget of solid watter  for each column (W/m2)
     38  REAL, SAVE, ALLOCATABLE :: d_h_qs_col(:)    ! enthalpy budget of cloud ice  for each column (W/m2)
     39  REAL, SAVE, ALLOCATABLE :: d_h_qbs_col(:)    ! enthalpy budget of blowing snow for each column (W/m2)
    3840  REAL, SAVE, ALLOCATABLE :: d_h_col(:)       ! total enthalpy budget for each column (W/m2)
    39   !$OMP THREADPRIVATE(d_qw_col, d_ql_col, d_qs_col, d_qt_col, d_ek_col, d_h_dair_col)
    40   !$OMP THREADPRIVATE(d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_col)
     41  !$OMP THREADPRIVATE(d_qw_col, d_ql_col, d_qs_col, d_qbs_col, d_qt_col, d_ek_col, d_h_dair_col)
     42  !$OMP THREADPRIVATE(d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_qbs_col, d_h_col)
    4143
    4244  ! Outputs used in cloudth_vert to extract the moments of the horizontal and
    … …  
    173175
    174176    allocate (bils_ec(klon),bils_ech(klon),bils_tke(klon),bils_diss(klon),bils_kinetic(klon),bils_enthalp(klon),bils_latent(klon))
    175     allocate (d_qw_col(klon), d_ql_col(klon), d_qs_col(klon), d_qt_col(klon), d_ek_col(klon), d_h_dair_col(klon) &
    176   &         , d_h_qw_col(klon), d_h_ql_col(klon), d_h_qs_col(klon), d_h_col(klon))
    177     d_qw_col=0. ; d_ql_col=0. ; d_qs_col=0. ; d_qt_col=0. ; d_ek_col=0. ; d_h_dair_col =0.
    178     d_h_qw_col=0. ; d_h_ql_col=0. ; d_h_qs_col=0. ; d_h_col=0.
     177    allocate (d_qw_col(klon), d_ql_col(klon), d_qs_col(klon), d_qbs_col(klon), d_qt_col(klon), d_ek_col(klon), d_h_dair_col(klon) &
     178  &         , d_h_qw_col(klon), d_h_ql_col(klon), d_h_qs_col(klon), d_h_qbs_col(klon), d_h_col(klon))
     179    d_qw_col=0. ; d_ql_col=0. ; d_qs_col=0. ; d_qbs_col=0. ; d_qt_col=0. ; d_ek_col=0. ; d_h_dair_col =0.
     180    d_h_qw_col=0. ; d_h_ql_col=0. ; d_h_qs_col=0. ; d_h_qbs_col=0. ; d_h_col=0.
    179181
    180182    ! Outputs used in cloudth_vert
    … …  
    223225    deallocate(sza_o)
    224226    deallocate (bils_ec,bils_ech,bils_tke,bils_diss,bils_kinetic,bils_enthalp,bils_latent)
    225     deallocate (d_qw_col, d_ql_col, d_qs_col, d_qt_col, d_ek_col, d_h_dair_col &
    226   &           , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_col)
     227    deallocate (d_qw_col, d_ql_col, d_qs_col, d_qbs_col, d_qt_col, d_ek_col, d_h_dair_col &
     228  &           , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_qbs_col, d_h_col)
    227229
    228230    ! Outputs used in cloudth_vert

  • LMDZ6/branches/blowing_snow/libf/phylmdiso/physiq_mod.F90

    r4504 r4506  
    7979    USE wxios, ONLY: g_ctx, wxios_set_context
    8080#endif
    81     USE atke_turbulence_ini_mod, ONLY : atke_ini
    82     USE lscp_ini_mod, ONLY : lscp_ini
    8381    USE lscp_mod, ONLY : lscp
    8482    USE wake_ini_mod, ONLY : wake_ini
    8583    USE yamada_ini_mod, ONLY : yamada_ini
     84    USE atke_turbulence_ini_mod, ONLY : atke_ini
    8685    USE thermcell_ini_mod, ONLY : thermcell_ini
     86    USE blowing_snow_ini_mod, ONLY : blowing_snow_ini , qbst_bs
     87    USE lscp_ini_mod, ONLY : lscp_ini
    8788
    8889    !USE cmp_seri_mod
    … …  
    183184       ! [Variables internes non sauvegardees de la physique]
    184185       ! Variables locales pour effectuer les appels en serie
    185        t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,tr_seri,rneb_seri, &
     186       t_seri,q_seri,ql_seri,qs_seri,qbs_seri,u_seri,v_seri,tr_seri,rneb_seri, &
    186187       ! Dynamic tendencies (diagnostics)
    187        d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_u_dyn,d_v_dyn,d_tr_dyn,d_rneb_dyn, &
    188        d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d, &
     188       d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_qbs_dyn,d_u_dyn,d_v_dyn,d_tr_dyn,d_rneb_dyn, &
     189       d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d,d_qbs_dyn2d, &
    189190       ! Physic tendencies
    190191       d_t_con,d_q_con,d_u_con,d_v_con, &
    … …  
    203204       plul_st,plul_th, &
    204205       !
    205        d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_t_diss, &
     206       d_t_vdf,d_q_vdf, d_qbs_vdf, d_u_vdf,d_v_vdf,d_t_diss, &
    206207       d_t_vdf_x, d_t_vdf_w, &
    207208       d_q_vdf_x, d_q_vdf_w, &
    208209       d_ts, &
     210       !
     211       d_t_bs,d_q_bs,d_qbs_bs, &
    209212       !
    210213!       d_t_oli,d_u_oli,d_v_oli, &
    … …  
    254257       cldh, cldl,cldm, cldq, cldt,      &
    255258       JrNt,                             &
    256        dthmin, evap, fder, plcl, plfc,   &
    257        prw, prlw, prsw,                  &
     259       dthmin, evap, snowerosion,fder, plcl, plfc,   &
     260       prw, prlw, prsw, prbsw,                  &
    258261       s_lcl, s_pblh, s_pblt, s_therm,   &
    259262       cdragm, cdragh,                   &
    … …  
    343346       fsolsw, wfbils, wfbilo,  &
    344347       wfevap, wfrain, wfsnow,  & 
    345        prfl, psfl, fraca, Vprecip,  &
     348       prfl, psfl,bsfl, fraca, Vprecip,  &
    346349       zw2,  &
    347350       !
    … …  
    526529    !======================================================================
    527530    !
    528     ! indices de traceurs eau vapeur, liquide, glace, fraction nuageuse LS (optional)
    529     INTEGER,SAVE :: ivap, iliq, isol, irneb
    530 !$OMP THREADPRIVATE(ivap, iliq, isol, irneb)
     531    ! indices de traceurs eau vapeur, liquide, glace, fraction nuageuse LS (optional), blowing snow (optional)
     532    INTEGER,SAVE :: ivap, iliq, isol, irneb, ibs
     533!$OMP THREADPRIVATE(ivap, iliq, isol, irneb, ibs)
    531534    !
    532535    !
    … …  
    905908    REAL dialiq(klon,klev)  ! eau liquide nuageuse
    906909    REAL diafra(klon,klev)  ! fraction nuageuse
    907     REAL cldliq(klon,klev)  ! eau liquide nuageuse
     910    REAL radocond(klon,klev)  ! eau condensee nuageuse
    908911    !
    909912    !XXX PB
    910913    REAL fluxq(klon,klev, nbsrf)   ! flux turbulent d'humidite
     914    REAL fluxqbs(klon,klev, nbsrf)   ! flux turbulent de neige soufflee
    911915    !
    912916    REAL zxfluxt(klon, klev)
    913917    REAL zxfluxq(klon, klev)
     918    REAL zxfluxqbs(klon,klev)
    914919    REAL zxfluxu(klon, klev)
    915920    REAL zxfluxv(klon, klev)
    … …  
    10091014    !
    10101015    ! tendance nulles
    1011     REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0, dqi0
     1016    REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0, dqi0, dqbs0
    10121017#ifdef ISO
    10131018    REAL, dimension(ntraciso,klon,klev):: dxt0, dxtl0, dxti0
    … …  
    11561161    REAL ztsol(klon)
    11571162    REAL q2m(klon,nbsrf)  ! humidite a 2m
     1163    REAL fsnowerosion(klon,nbsrf) ! blowing snow flux at surface
     1164    REAL qbsfra  ! blowing snow fraction
    11581165#ifdef ISO
    11591166    REAL d_xtw(ntraciso),d_xtl(ntraciso), d_xts(ntraciso)
    … …  
    13721379       isol = strIdx(tracers(:)%name, addPhase('H2O', 's'))
    13731380       irneb= strIdx(tracers(:)%name, addPhase('H2O', 'r'))
     1381       ibs  = strIdx(tracers(:)%name, addPhase('H2O', 'b'))
    13741382       CALL init_etat0_limit_unstruct
    13751383       IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
    … …  
    14151423       ENDIF
    14161424
    1417 
    14181425       IF (ok_ice_sursat.AND.(iflag_ice_thermo.EQ.0)) THEN
    14191426          WRITE (lunout, *) ' ok_ice_sursat=y requires iflag_ice_thermo=1 as well'
    … …  
    14221429       ENDIF
    14231430
    1424        IF (ok_ice_sursat.AND.(nqo.NE.4)) THEN
     1431       IF (ok_ice_sursat.AND.(nqo.LT.4)) THEN
    14251432          WRITE (lunout, *) ' ok_ice_sursat=y requires 4 H2O tracers ', &
    14261433               '(H2O_g, H2O_l, H2O_s, H2O_r) but nqo=', nqo, '. Might as well stop here.'
    … …  
    14411448       ENDIF
    14421449
    1443        IF (ok_bs) THEN
     1450        IF (ok_bs) THEN
    14441451          abort_message='blowing snow cannot be activated with water isotopes yet'
    14451452          CALL abort_physic(modname,abort_message, 1)
    1446        ENDIF
     1453         IF ((ok_ice_sursat.AND.nqo .LT.5).OR.(.NOT.ok_ice_sursat.AND.nqo.LT.4)) THEN
     1454             WRITE (lunout, *) 'activation of blowing snow needs a specific H2O tracer', &
     1455                               'but nqo=', nqo
     1456             abort_message='see above'
     1457             CALL abort_physic(modname,abort_message, 1)
     1458         ENDIF
     1459        ENDIF
    14471460       Ncvpaseq1 = 0
    14481461       dnwd0=0.0
    … …  
    18981911           CALL lscp_ini(pdtphys,ok_ice_sursat)
    18991912       endif
     1913       CALL blowing_snow_ini(prt_level,lunout, &
     1914                             RCPD, RLSTT, RLVTT, RLMLT, &
     1915                             RVTMP2, RTT,RD,RG)
    19001916
    19011917!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    19241940       CALL phys_output_write(itap, pdtphys, paprs, pphis,                    &
    19251941                              pplay, lmax_th, aerosol_couple,                 &
    1926                               ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, ok_sync,&
     1942                              ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, ibs,  ok_sync,&
    19271943                              ptconv, read_climoz, clevSTD,                   &
    19281944                              ptconvth, d_u, d_t, qx, d_qx, zmasse,           &
    … …  
    24112427    dql0(:,:)=0.
    24122428    dqi0(:,:)=0.
     2429    dqbs0(:,:)=0.
    24132430#ifdef ISO
    24142431      dxt0(:,:,:)=0.
    … …  
    24682485          q_seri(i,k)  = qx(i,k,ivap)
    24692486          ql_seri(i,k) = qx(i,k,iliq)
     2487          qbs_seri(i,k) = 0.
    24702488          !CR: ATTENTION, on rajoute la variable glace
    24712489          IF (nqo.EQ.2) THEN             !--vapour and liquid only
    … …  
    24752493             qs_seri(i,k) = qx(i,k,isol)
    24762494             rneb_seri(i,k) = 0.
    2477           ELSE IF (nqo.EQ.4) THEN        !--vapour, liquid, ice and rneb
     2495          ELSE IF (nqo.GE.4) THEN        !--vapour, liquid, ice and rneb and blowing snow
    24782496             qs_seri(i,k) = qx(i,k,isol)
     2497             IF (ok_ice_sursat) THEN
    24792498             rneb_seri(i,k) = qx(i,k,irneb)
     2499             ENDIF
     2500             IF (ok_bs) THEN
     2501             qbs_seri(i,k)= qx(i,k,ibs)
     2502             ENDIF
     2503
    24802504          ENDIF
     2505
     2506
    24812507       ENDDO
    24822508    ENDDO
    … …  
    25192545    qql1(:)=0.0
    25202546    DO k = 1, klev
    2521       qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
     2547      qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k)+qbs_seri(:,k))*zmasse(:,k)
    25222548    ENDDO
    25232549#ifdef ISO
    … …  
    26262652       d_ql_dyn(:,:) = (ql_seri(:,:)-ql_ancien(:,:))/phys_tstep
    26272653       d_qs_dyn(:,:) = (qs_seri(:,:)-qs_ancien(:,:))/phys_tstep
     2654       d_qbs_dyn(:,:) = (qbs_seri(:,:)-qbs_ancien(:,:))/phys_tstep
    26282655       CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
    26292656       d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/phys_tstep
    … …  
    26322659       CALL water_int(klon,klev,qs_seri,zmasse,zx_tmp_fi2d)
    26332660       d_qs_dyn2d(:)=(zx_tmp_fi2d(:)-prsw_ancien(:))/phys_tstep
     2661       CALL water_int(klon,klev,qbs_seri,zmasse,zx_tmp_fi2d)
     2662       d_qbs_dyn2d(:)=(zx_tmp_fi2d(:)-prbsw_ancien(:))/phys_tstep
    26342663       ! !! RomP >>>   td dyn traceur
    26352664       IF (nqtot > nqo) d_tr_dyn(:,:,:)=(tr_seri(:,:,:)-tr_ancien(:,:,:))/phys_tstep
    … …  
    27182747       d_ql_dyn2d(:) = 0.0
    27192748       d_qs_dyn2d(:) = 0.0
     2749       d_qbs_dyn2d(:)= 0.0
    27202750
    27212751#ifdef ISO
    … …  
    27412771       ! !! RomP <<<
    27422772       d_rneb_dyn(:,:)=0.0
     2773       d_qbs_dyn(:,:)=0.0
    27432774       ancien_ok = .TRUE.
    27442775    ENDIF
    … …  
    28602891
    28612892     CALL add_phys_tend &
    2862             (du0,dv0,d_t_eva,d_q_eva,d_ql_eva,d_qi_eva,paprs,&
     2893            (du0,dv0,d_t_eva,d_q_eva,d_ql_eva,d_qi_eva,dqbs0,paprs,&
    28632894               'eva',abortphy,flag_inhib_tend,itap,0 &
    28642895#ifdef ISO
    … …  
    30293060            longitude_deg, latitude_deg, rugoro,  zrmu0,      &
    30303061            sollwdown,    cldt,      &
    3031             rain_fall, snow_fall, solsw,   solswfdiff, sollw,     &
     3062            rain_fall, snow_fall, bs_fall, solsw,   solswfdiff, sollw,     &
    30323063            gustiness,                                &
    3033             t_seri,    q_seri,    u_seri,  v_seri,    &
     3064            t_seri,    q_seri,   qbs_seri, u_seri,  v_seri,    &
    30343065                                !nrlmd+jyg<
    30353066            wake_deltat, wake_deltaq, wake_cstar, wake_s, &
    … …  
    30423073                                !albedo SB >>>
    30433074                                ! albsol1,   albsol2,   sens,    evap,      &
    3044             albsol_dir,   albsol_dif,   sens,    evap,   & 
     3075            albsol_dir,   albsol_dif,   sens,    evap, snowerosion, & 
    30453076                                !albedo SB <<<
    30463077            albsol3_lic,runoff,   snowhgt,   qsnow, to_ice, sissnow, &
    30473078            zxtsol,    zxfluxlat, zt2m,    qsat2m,  zn2mout, &
    3048             d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf, d_t_diss, &
     3079            d_t_vdf,   d_q_vdf, d_qbs_vdf,  d_u_vdf, d_v_vdf, d_t_diss, &
    30493080                                !nrlmd<
    30503081                                !jyg<
    … …  
    30723103            fluxt,   fluxu,  fluxv, &
    30733104            dsens,     devap,     zxsnow, &
    3074             zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke, &
     3105            zxfluxt,   zxfluxq,  zxfluxqbs,  q2m, fluxq, fluxqbs, pbl_tke, &
    30753106                                !nrlmd+jyg<
    30763107            wake_delta_pbl_TKE, &
    … …  
    31843215       IF (klon_glo==1) THEN
    31853216          CALL add_pbl_tend &
    3186                (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
     3217               (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,d_qbs_vdf,paprs,&
    31873218               'vdf',abortphy,flag_inhib_tend,itap &
    31883219#ifdef ISO
    … …  
    31923223       ELSE
    31933224          CALL add_phys_tend &
    3194                (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
     3225               (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,d_qbs_vdf,paprs,&
    31953226               'vdf',abortphy,flag_inhib_tend,itap,0 &
    31963227#ifdef ISO
    … …  
    32543285
    32553286    ENDIF
     3287
     3288    ! ==================================================================
     3289    ! Blowing snow sublimation and sedimentation
     3290
     3291    d_t_bs(:,:)=0.
     3292    d_q_bs(:,:)=0.
     3293    d_qbs_bs(:,:)=0.
     3294    bsfl(:,:)=0.
     3295    bs_fall(:)=0.
     3296    IF (ok_bs) THEN
     3297
     3298     CALL call_blowing_snow_sublim_sedim(klon,klev,phys_tstep,t_seri,q_seri,qbs_seri,pplay,paprs, &
     3299                                        d_t_bs,d_q_bs,d_qbs_bs,bsfl,bs_fall)
     3300
     3301     CALL add_phys_tend &
     3302               (du0,dv0,d_t_bs,d_q_bs,dql0,dqi0,d_qbs_bs,paprs,&
     3303               'bs',abortphy,flag_inhib_tend,itap,0  &
     3304#ifdef ISO                                       
     3305       &    ,dxt0,dxtl0,dxti0 &                   
     3306#endif                                       
     3307       &   )
     3308
     3309    ENDIF
     3310
    32563311    ! =================================================================== c
    32573312    !   Calcul de Qsat
    … …  
    39253980!!
    39263981!!
    3927     CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, paprs, &
     3982    CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, dqbs0, paprs, &
    39283983         'convection',abortphy,flag_inhib_tend,itap,0 &
    39293984#ifdef ISO
    … …  
    41574212       !-----------------------------------------------------------------------
    41584213       ! ajout des tendances des poches froides
    4159        CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake', &
     4214       CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,dqbs0,paprs,'wake', &
    41604215            abortphy,flag_inhib_tend,itap,0 &
    41614216#ifdef ISO
    … …  
    44364491          !
    44374492          CALL add_phys_tend(d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs,  &
    4438                              dql0,dqi0,paprs,'thermals', abortphy,flag_inhib_tend,itap,0 &
     4493                             dql0,dqi0,dqbs0,paprs,'thermals', abortphy,flag_inhib_tend,itap,0 &
    44394494#ifdef ISO
    44404495     &    ,d_xt_ajs,dxtl0,dxti0 &
    … …  
    45524607          !--------------------------------------------------------------------
    45534608          ! ajout des tendances de l'ajustement sec ou des thermiques
    4554           CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs, &
     4609          CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,dqbs0,paprs, &
    45554610               'ajsb',abortphy,flag_inhib_tend,itap,0 &
    45564611#ifdef ISO
    … …  
    46964751    CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay, &
    46974752         t_seri, q_seri,ptconv,ratqs, &
    4698          d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, rneb_seri, &
    4699          cldliq, picefra, rain_lsc, snow_lsc, &
     4753         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, rneb_seri, & 
     4754         radocond, picefra, rain_lsc, snow_lsc, &
    47004755         frac_impa, frac_nucl, beta_prec_fisrt, &
    47014756         prfl, psfl, rhcl,  &
    … …  
    47094764    CALL fisrtilp(phys_tstep,paprs,pplay, &
    47104765         t_seri, q_seri,ptconv,ratqs, &
    4711          d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, &
     4766         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, radocond, &
    47124767         rain_lsc, snow_lsc, &
    47134768         pfrac_impa, pfrac_nucl, pfrac_1nucl, &
    … …  
    47574812!    write(*,9000) "rcpv","rcw",rcpv,rcw,rcs,t_seri(1,1)
    47584813!-JLD
    4759     CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs, &
     4814    CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,dqbs0,paprs, &
    47604815         'lsc',abortphy,flag_inhib_tend,itap,0 &
    47614816#ifdef ISO
    … …  
    47944849    ENDIF
    47954850
    4796     !---------------------------------------------------------------------------
     4851
     4852!---------------------------------------------------------------------------
    47974853    DO k = 1, klev
    47984854       DO i = 1, klon
    47994855          cldfra(i,k) = rneb(i,k)
    48004856          !CR: a quoi ca sert? Faut-il ajouter qs_seri?
    4801           IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
     4857          !EV: en effet etrange, j'ajouterais aussi qs_seri
     4858          !    plus largement, je nettoierais (enleverrais) ces lignes
     4859          IF (.NOT.new_oliq) radocond(i,k) = ql_seri(i,k)
    48024860       ENDDO
    48034861    ENDDO
    48044862
    48054863
    4806 
     4864    ! Option to activate the radiative effect of blowing snow (ok_rad_bs)
     4865    ! makes sense only if the new large scale condensation scheme is active
     4866    ! with the ok_icefra_lscp flag active as well
     4867
     4868    IF (ok_bs .AND. ok_rad_bs) THEN
     4869       IF (ok_new_lscp .AND. ok_icefra_lscp) THEN
     4870           DO k=1,klev
     4871             DO i=1,klon
     4872                radocond(i,k)=radocond(i,k)+qbs_seri(i,k)
     4873                picefra(i,k)=(radocond(i,k)*picefra(i,k)+qbs_seri(i,k))/(radocond(i,k))
     4874                qbsfra=min(qbs_seri(i,k)/qbst_bs,1.0)
     4875                cldfra(i,k)=max(cldfra(i,k),qbsfra)
     4876             ENDDO
     4877           ENDDO
     4878       ELSE
     4879          WRITE(lunout,*)"PAY ATTENTION, you try to activate the radiative effect of blowing snow"
     4880          WRITE(lunout,*)"with ok_new_lscp=false and/or ok_icefra_lscp=false"
     4881          abort_message='inconsistency in cloud phase for blowing snow'
     4882          CALL abort_physic(modname,abort_message,1)
     4883       ENDIF
     4884
     4885    ENDIF
    48074886#ifdef ISO     
    48084887!#ifdef ISOVERIF
    … …  
    49675046          DO i = 1, klon
    49685047             IF (diafra(i,k).GT.cldfra(i,k)) THEN
    4969                 cldliq(i,k) = dialiq(i,k)
     5048                radocond(i,k) = dialiq(i,k)
    49705049                cldfra(i,k) = diafra(i,k)
    49715050             ENDIF
    … …  
    50045083                DO i=1,klon
    50055084                   IF (ptconv(i,k).AND.ptconvth(i,k)) THEN
    5006                       cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
     5085                      radocond(i,k)=radocond(i,k)+rnebcon(i,k)*clwcon(i,k)
    50075086                      cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
    50085087                   ELSE IF (ptconv(i,k)) THEN
    50095088                      cldfra(i,k)=rnebcon(i,k)
    5010                       cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
     5089                      radocond(i,k)=rnebcon(i,k)*clwcon(i,k)
    50115090                   ENDIF
    50125091                ENDDO
    … …  
    50175096                DO i=1,klon
    50185097                   cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
    5019                    cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
     5098                   radocond(i,k)=radocond(i,k)+rnebcon(i,k)*clwcon(i,k)
    50205099                ENDDO
    50215100             ENDDO
    … …  
    50355114                   IF (ptconv(i,k).AND. .NOT.ptconvth(i,k)) THEN
    50365115                      cldfra(i,k)=rnebcon(i,k)
    5037                       cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
     5116                      radocond(i,k)=rnebcon(i,k)*clwcon(i,k)
    50385117                   ENDIF
    50395118                ENDDO
    … …  
    50465125          ! Ancienne version
    50475126          cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
    5048           cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
     5127          radocond(:,:)=radocond(:,:)+rnebcon(:,:)*clwcon(:,:)
    50495128       ENDIF
    50505129
    … …  
    50665145          DO i = 1, klon
    50675146             IF (diafra(i,k).GT.cldfra(i,k)) THEN
    5068                 cldliq(i,k) = dialiq(i,k)
     5147                radocond(i,k) = dialiq(i,k)
    50695148                cldfra(i,k) = diafra(i,k)
    50705149             ENDIF
    … …  
    54415520          ENDIF
    54425521          CALL newmicro (flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, &
    5443                paprs, pplay, t_seri, cldliq, picefra, cldfra, &
     5522               paprs, pplay, t_seri, radocond, picefra, cldfra, &
    54445523               cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
    54455524               flwp, fiwp, flwc, fiwc, &
    … …  
    54495528       ELSE
    54505529          CALL nuage (paprs, pplay, &
    5451                t_seri, cldliq, picefra, cldfra, cldtau, cldemi, &
     5530               t_seri, radocond, picefra, cldfra, cldtau, cldemi, &
    54525531               cldh, cldl, cldm, cldt, cldq, &
    54535532               ok_aie, &
    … …  
    57585837    ENDDO
    57595838
    5760     CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy,flag_inhib_tend,itap,0 &
     5839    CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,dqbs0,paprs,'SW',abortphy,flag_inhib_tend,itap,0 &
    57615840#ifdef ISO
    57625841     &    ,dxt0,dxtl0,dxti0 &
    … …  
    57645843     &  )
    57655844    CALL prt_enerbil('SW',itap)
    5766     CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy,flag_inhib_tend,itap,0 &
     5845    CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,dqbs0,paprs,'LW',abortphy,flag_inhib_tend,itap,0 &
    57675846#ifdef ISO
    57685847     &    ,dxt0,dxtl0,dxti0 &
    … …  
    58125891          ! -> condition on zrel_oro can deactivate the drag on tilted planar terrains
    58135892          !    such as ice sheets (work by V. Wiener)
    5814           ! zpmm_orodr_t and zstd_orodr_t are activation thresholds set by F. Lott to 
     5893          ! zpmm_orodr_t and zstd_orodr_t are activation thresholds set by F. Lott to
    58155894          ! earn computation time but they are not physical.
    58165895          IF (((zpic(i)-zmea(i)).GT.zpmm_orodr_t).AND.(zstd(i).GT.zstd_orodr_t).AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
    … …  
    58435922       !-----------------------------------------------------------------------
    58445923       ! ajout des tendances de la trainee de l'orographie
    5845        CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro', &
     5924       CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,dqbs0,paprs,'oro', &
    58465925            abortphy,flag_inhib_tend,itap,0 &
    58475926#ifdef ISO
    … …  
    58985977
    58995978       ! ajout des tendances de la portance de l'orographie
    5900        CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, paprs, &
     5979       CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, dqbs0,paprs, &
    59015980            'lif', abortphy,flag_inhib_tend,itap,0 &
    59025981#ifdef ISO
    … …  
    59276006       d_t_hin(:, :)=0.
    59286007       CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
    5929             dqi0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0 &
     6008            dqi0, dqbs0,paprs, 'hin', abortphy,flag_inhib_tend,itap,0 &
    59306009#ifdef ISO
    59316010     &    ,dxt0,dxtl0,dxti0 &
    … …  
    59496028       ENDDO
    59506029
    5951        CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, &
     6030       CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, dqbs0, &
    59526031            paprs, 'front_gwd_rando', abortphy,flag_inhib_tend,itap,0 &
    59536032#ifdef ISO
    … …  
    59626041            rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
    59636042            du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
    5964        CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, &
     6043       CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, dqbs0, &
    59656044            paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend,itap,0 &
    59666045#ifdef ISO
    … …  
    60216100       d_xt_ch4_dtime(:,:,:) = d_xt_ch4(:,:,:)*phys_tstep
    60226101#endif
    6023        CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, paprs, &
     6102       CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, dqbs0, paprs, &
    60246103            'q_ch4', abortphy,flag_inhib_tend,itap,0 &
    60256104#ifdef ISO
    … …  
    60876166! car on peut s'attendre a ce que les petites echelles produisent aussi de la TKE
    60886167! Mais attention, cela ne va pas dans le sens de la conservation de l'energie!
    6089           IF ((zstd(i).GT.1.0).AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
     6168          IF ((zstd(i).GT.1.0) .AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
    60906169             itest(i)=1
    60916170             igwd=igwd+1
    … …  
    63476426    !
    63486427
    6349     IF (type_trac=='repr') THEN
     6428    IF (type_trac == 'repr') THEN
    63506429!MM pas d'impact, car on recupere q_seri,tr_seri,t_seri via phys_local_var_mod
    63516430!MM                               dans Reprobus
    … …  
    63976476         presnivs, pphis,     pphi,     albsol1, &
    63986477         sh_in,   ch_in,    rhcl,      cldfra,   rneb, &
    6399          diafra,   cldliq,    itop_con, ibas_con, &
     6478         diafra,   radocond,    itop_con, ibas_con, &
    64006479         pmflxr,   pmflxs,    prfl,     psfl, &
    64016480         da,       phi,       mp,       upwd, &
    … …  
    64146493
    64156494     CALL add_phys_tend &
    6416             (du0,dv0,dt0,d_q_rep,d_ql_rep,d_qi_rep,paprs,&
     6495            (du0,dv0,dt0,d_q_rep,d_ql_rep,d_qi_rep,dqbs0,paprs,&
    64176496             'rep',abortphy,flag_inhib_tend,itap,0)
    64186497        IF (abortphy==1) Print*,'ERROR ABORT REP'
    … …  
    64916570    !   prlw = colonne eau liquide
    64926571    !   prlw = colonne eau solide
     6572    !   prbsw = colonne neige soufflee
    64936573    prw(:) = 0.
    64946574    prlw(:) = 0.
    64956575    prsw(:) = 0.
     6576    prbsw(:) = 0.
    64966577    DO k = 1, klev
    64976578       prw(:)  = prw(:)  + q_seri(:,k)*zmasse(:,k)
    64986579       prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
    64996580       prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
     6581       prbsw(:)= prbsw(:) + qbs_seri(:,k)*zmasse(:,k)
    65006582    ENDDO
    65016583
    … …  
    65686650          ENDIF
    65696651          !--ice_sursat: nqo=4, on ajoute rneb
    6570           IF (nqo == 4) THEN
     6652          IF (nqo.ge.4 .and. ok_ice_sursat) THEN
    65716653             d_qx(i,k,irneb) = ( rneb_seri(i,k) - qx(i,k,irneb) ) / phys_tstep
    65726654          ENDIF
     6655
     6656           IF (nqo.ge.4 .and. ok_bs) THEN
     6657             d_qx(i,k,ibs) = ( qbs_seri(i,k) - qx(i,k,ibs) ) / phys_tstep
     6658          ENDIF
     6659
     6660
    65736661       ENDDO
    65746662    ENDDO
    … …  
    66506738    ql_ancien(:,:) = ql_seri(:,:)
    66516739    qs_ancien(:,:) = qs_seri(:,:)
     6740    qbs_ancien(:,:) = qbs_seri(:,:)
    66526741    rneb_ancien(:,:) = rneb_seri(:,:)
    66536742#ifdef ISO
    … …  
    66596748    CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
    66606749    CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
     6750    CALL water_int(klon,klev,qbs_ancien,zmasse,prbsw_ancien)
    66616751    ! !! RomP >>>
    66626752    IF (nqtot > nqo) tr_ancien(:,:,:) = tr_seri(:,:,:)
    … …  
    67876877    CALL phys_output_write(itap, pdtphys, paprs, pphis,  &
    67886878         pplay, lmax_th, aerosol_couple,                 &
    6789          ok_ade, ok_aie, ok_volcan, ivap, iliq, isol,    &
     6879         ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, ibs,   &
    67906880         ok_sync, ptconv, read_climoz, clevSTD,          &
    67916881         ptconvth, d_u, d_t, qx, d_qx, zmasse,           &

  • LMDZ6/branches/blowing_snow/libf/phylmdiso/surf_land_mod.F90

    r4285 r4506  
    1111       rlon, rlat, yrmu0, &
    1212       debut, lafin, zlev, ccanopy, swnet, lwnet, albedo, &
    13        tsurf, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
     13       tsurf, p1lay, cdragh, cdragm, precip_rain, precip_snow, precip_bs, temp_air, spechum, &
    1414       AcoefH, AcoefQ, BcoefH, BcoefQ, &
    1515       AcoefU, AcoefV, BcoefU, BcoefV, &
    … …  
    1717       lwdown_m, q2m, t2m, &
    1818       snow, qsol, agesno, tsoil, &
    19        z0m, z0h, SFRWL, alb_dir_new, alb_dif_new, evap, fluxsens, fluxlat, &   
     19       z0m, z0h, SFRWL, alb_dir_new, alb_dif_new, evap, fluxsens, fluxlat, fluxbs, &   
    2020       qsurf, tsurf_new, dflux_s, dflux_l, &
    2121       flux_u1, flux_v1 , &
    … …  
    9494    REAL, DIMENSION(klon), INTENT(IN)       :: p1lay
    9595    REAL, DIMENSION(klon), INTENT(IN)       :: cdragh, cdragm
    96     REAL, DIMENSION(klon), INTENT(IN)       :: precip_rain, precip_snow
     96    REAL, DIMENSION(klon), INTENT(IN)       :: precip_rain, precip_snow, precip_bs
    9797    REAL, DIMENSION(klon), INTENT(IN)       :: temp_air, spechum
    9898    REAL, DIMENSION(klon), INTENT(IN)       :: AcoefH, AcoefQ, BcoefH, BcoefQ
    … …  
    129129!albedo SB <<<
    130130    REAL, DIMENSION(klon), INTENT(OUT)       :: evap
    131     REAL, DIMENSION(klon), INTENT(OUT)       :: fluxsens, fluxlat
     131    REAL, DIMENSION(klon), INTENT(OUT)       :: fluxsens, fluxlat, fluxbs
    132132    REAL, DIMENSION(klon), INTENT(OUT)       :: qsurf
    133133    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
    … …  
    152152    REAL, DIMENSION(klon) :: tsol_rad, emis_new ! output from interfsol not used
    153153    REAL, DIMENSION(klon) :: u0, v0     ! surface speed
     154    REAL, DIMENSION(klon) :: precip_totsnow     ! total solid precip
    154155    INTEGER               :: i
    155156
    … …  
    166167#endif 
    167168
     169!****************************************************************************************
     170!Total solid precip
     171
     172IF (ok_bs) THEN
     173precip_totsnow=precip_snow+precip_bs
     174ELSE
     175precip_totsnow=precip_snow
     176ENDIF
     177!****************************************************************************************
    168178#ifdef ISO
    169179#ifdef ISOVERIF
    … …  
    228238            zlev,  u1, v1, gustiness, temp_air, spechum, epot_air, ccanopy, &
    229239            cdragh, AcoefH, AcoefQ, BcoefH, BcoefQ, &
    230             precip_rain, precip_snow, lwdown_m, swnet, swdown, &
     240            precip_rain, precip_totsnow, lwdown_m, swnet, swdown, &
    231241            pref_tmp, q2m, t2m, &
    232             evap, fluxsens, fluxlat, &             
     242            evap, fluxsens, fluxlat,fluxbs, &             
    233243            tsol_rad, tsurf_new, alb1_new, alb2_new, &
    234244            emis_new, z0m, z0h, qsurf, &
    … …  
    277287#endif
    278288       CALL surf_land_bucket(itime, jour, knon, knindex, debut, dtime,&
    279             tsurf, p1lay, cdragh, precip_rain, precip_snow, temp_air, &
     289            tsurf, p1lay, cdragh, precip_rain, precip_totsnow, temp_air, &
    280290            spechum, AcoefH, AcoefQ, BcoefH, BcoefQ, pref, &
    281291            u1, v1, gustiness, rugoro, swnet, lwnet, &
    … …  
    293303    ENDIF ! ok_veget
    294304
     305        ! blowing snow not treated yet over land
     306        fluxbs(:)=0.
    295307!****************************************************************************************
    296308! Calculation for all land models

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

    r4285 r4506  
    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, &
    1717       ps, u1, v1, gustiness, rugoro, pctsrf, &
    18        snow, qsurf, qsol, agesno, &
    19        tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, fluxsens, fluxlat, &
     18       snow, qsurf, qsol, qbs1, agesno, &
     19       tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, fluxsens, fluxlat, fluxbs, &
    2020       tsurf_new, dflux_s, dflux_l, &
    2121       alt, slope, cloudf, &
    … …  
    3434    USE cpl_mod,          ONLY : cpl_send_landice_fields
    3535    USE calcul_fluxs_mod
    36     USE phys_output_var_mod
     36    USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic
     37    USE phys_output_var_mod, ONLY : snow_o,zfra_o
    3738#ifdef ISO   
    3839    USE fonte_neige_mod,  ONLY : xtrun_off_lic
    … …  
    4849!FC
    4950    USE ioipsl_getin_p_mod, ONLY : getin_p
    50 
     51    USE blowing_snow_ini_mod, ONLY : zeta_bs, pbst_bs, prt_bs
    5152
    5253#ifdef CPP_INLANDSIS
    … …  
    7273    REAL, DIMENSION(klon), INTENT(IN)             :: p1lay
    7374    REAL, DIMENSION(klon), INTENT(IN)             :: cdragh, cdragm
    74     REAL, DIMENSION(klon), INTENT(IN)             :: precip_rain, precip_snow
     75    REAL, DIMENSION(klon), INTENT(IN)             :: precip_rain, precip_snow, precip_bs
    7576    REAL, DIMENSION(klon), INTENT(IN)             :: temp_air, spechum
    7677    REAL, DIMENSION(klon), INTENT(IN)             :: AcoefH, AcoefQ
    … …  
    7879    REAL, DIMENSION(klon), INTENT(IN)             :: AcoefU, AcoefV, BcoefU, BcoefV
    7980    REAL, DIMENSION(klon), INTENT(IN)             :: ps
    80     REAL, DIMENSION(klon), INTENT(IN)             :: u1, v1, gustiness
     81    REAL, DIMENSION(klon), INTENT(IN)             :: u1, v1, gustiness, qbs1
    8182    REAL, DIMENSION(klon), INTENT(IN)             :: rugoro
    8283    REAL, DIMENSION(klon,nbsrf), INTENT(IN)       :: pctsrf
    … …  
    118119!albedo SB <<<
    119120    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
     121    REAL, DIMENSION(klon), INTENT(OUT)            :: fluxbs
    120122    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
    121123    REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l     
    … …  
    179181
    180182    REAL,DIMENSION(klon) :: alb1,alb2
     183    REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow
    181184    REAL, DIMENSION (klon,6) :: alb6
     185    REAL                   :: rho0, rhoice, ustart0, hsalt, esalt, qsalt
     186    REAL                   :: tau_dens, tau_dens0, tau_densmin, rhomax, rhohard
     187    REAL, DIMENSION(klon)  :: ws1, rhos, ustart
    182188! End definition
    183189!****************************************************************************************
    … …  
    214220           PRINT*, 'alb_nir_sno_lic',alb_nir_sno_lic
    215221 
    216 !  z0m=1.e-3
    217 !  z0h = z0m
    218222  firstcall=.false.
    219223  ENDIF
    220224!******************************************************************************************
    221 !
     225
    222226! Initialize output variables
    223227    alb3(:) = 999999.
    224228    alb2(:) = 999999.
    225229    alb1(:) = 999999.
    226    
     230    fluxbs(:)=0. 
    227231    runoff(:) = 0.
    228232!****************************************************************************************
    … …  
    232236    radsol(:) = 0.0
    233237    radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
     238
     239!****************************************************************************************
    234240
    235241!****************************************************************************************
    … …  
    327333
    328334
    329 
    330335    ELSE
    331336
    … …  
    342347    IF (soil_model) THEN
    343348       CALL soil(dtime, is_lic, knon, snow, tsurf, qsol, &
    344     & longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux)
    345 
     349        & longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux)
    346350       cal(1:knon) = RCPD / soilcap(1:knon)
    347351       radsol(1:knon)  = radsol(1:knon) + soilflux(1:knon)
    … …  
    406410         flux_u1, flux_v1)
    407411
    408 !****************************************************************************************
    409 ! Calculate snow height, age, run-off,..
     412
     413!****************************************************************************************
     414! Calculate albedo
     415!
     416!****************************************************************************************
     417 
     418    CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)) 
     419
     420
     421! EV: following lines are obsolete since we set alb1 and alb2 to constant values
     422! I therefore comment them
     423!    alb1(1:knon) = alb_neig(1:knon)*zfra(1:knon) + &
     424!         0.6 * (1.0-zfra(1:knon))
     425!
     426!IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
     427!       alb1(1 : knon)  = 0.6 !IM cf FH/GK
     428!       alb1(1 : knon)  = 0.82
     429!       alb1(1 : knon)  = 0.77 !211003 Ksta0.77
     430!       alb1(1 : knon)  = 0.8 !KstaTER0.8 & LMD_ARMIP5
     431!IM: KstaTER0.77 & LMD_ARMIP6   
     432
     433! Attantion: alb1 and alb2 are not the same!
     434    alb1(1:knon)  = alb_vis_sno_lic
     435    alb2(1:knon)  = alb_nir_sno_lic
     436
     437
     438!****************************************************************************************
     439! Rugosity
     440!
     441!****************************************************************************************
     442    z0m = z0m_landice
     443    z0h = z0h_landice
     444    !z0m = SQRT(z0m**2+rugoro**2)
     445
     446
     447
     448  ! Simple blowing snow param
     449  if (ok_bs) then
     450       ustart0 = 0.211
     451       rhoice = 920.0
     452       rho0 = 200.0
     453       rhomax=450.0
     454       rhohard=400.0
     455       tau_dens0=86400.0*10.  ! 10 days by default, in s
     456       tau_densmin=86400.0 ! 1 days according to in situ obs by C. Amory
     457       do i = 1, knon
     458           ! estimation of snow density
     459           ! snow density increases with snow age and
     460           ! increases even faster in case of sedimentation of blowing snow
     461           tau_dens=max(tau_densmin, tau_dens0*exp(-abs(precip_bs(i))/pbst_bs-abs(precip_rain(i))/prt_bs))
     462           rhos(i)=rho0+(rhohard-rho0)*(1.-exp(-agesno(i)*86400.0/tau_dens))
     463           ! blowing snow flux formula used in MAR
     464           ws1(i)=(u1(i)**2+v1(i)**2)**0.5
     465           ustar(i)=(cdragm(i)*(u1(i)**2+v1(i)**2))**0.5
     466           ustart(i)=ustart0*exp(max(rhoice/rho0-rhoice/rhos(i),0.))*exp(max(0.,rhos(i)-rhomax))
     467           ! we have multiplied by exp to prevent erosion when rhos>rhomax (usefull till
     468           ! rhohard<450)
     469           esalt=1./(3.25*max(ustar(i),0.001))
     470           hsalt=0.08436*ustar(i)**1.27
     471           qsalt=(max(ustar(i)**2-ustart(i)**2,0.))/(RG*hsalt)*esalt
     472           !ep=qsalt*cdragm(i)*sqrt(u1(i)**2+v1(i)**2)
     473           fluxbs(i)= zeta_bs*p1lay(i)/RD/temp_air(i)*ws1(i)*cdragm(i)*(qbs1(i)-qsalt)
     474       enddo
     475
     476       ! for outputs
     477       do j = 1, knon
     478          i = knindex(j)
     479          zxustartlic(i) = ustart(j)
     480          zxrhoslic(i) = rhos(j)
     481       enddo
     482
     483  endif
     484
     485
     486
     487!****************************************************************************************
     488! Calculate surface snow amount
    410489!   
    411490!****************************************************************************************
     491    IF (ok_bs) THEN
     492      precip_totsnow=precip_snow+precip_bs
     493      evap_totsnow=evap-fluxbs ! flux bs is positive towards the surface (snow erosion)
     494    ELSE
     495      precip_totsnow=precip_snow
     496      evap_totsnow=evap
     497    ENDIF
     498
    412499    CALL fonte_neige(knon, is_lic, knindex, dtime, &
    413          tsurf, precip_rain, precip_snow, &
    414          snow, qsol, tsurf_new, evap &
     500         tsurf, precip_rain, precip_totsnow, &
     501         snow, qsol, tsurf_new, evap_totsnow &
    415502#ifdef ISO   
    416503     & ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag   &
    … …  
    446533
    447534
    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 
     535    WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.                                         
     536    zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0))) 
    479537
    480538
    … …  
    494552          run_off_lic_frac(j) = pctsrf(i,is_lic)
    495553       ENDDO
    496        
     554
    497555       CALL cpl_send_landice_fields(itime, knon, knindex, run_off_lic, run_off_lic_frac)
    498556    ENDIF
    … …  
    501559    runoff(1:knon)=run_off_lic(1:knon)/dtime
    502560
    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 !****************************************************************************************
    515561       snow_o=0.
    516562       zfra_o = 0.
    … …  
    553599!albedo SB <<<
    554600
    555  
    556  
    557601
    558602  END SUBROUTINE surf_landice

  • LMDZ6/branches/blowing_snow/libf/phylmdiso/surf_ocean_mod.F90

    r4374 r4506  
    1313       windsp, rmu0, fder, tsurf_in, &
    1414       itime, dtime, jour, knon, knindex, &
    15        p1lay, z1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
     15       p1lay, z1lay, cdragh, cdragm, precip_rain, precip_snow, precip_bs, temp_air, spechum, &
    1616       AcoefH, AcoefQ, BcoefH, BcoefQ, &
    1717       AcoefU, AcoefV, BcoefU, BcoefV, &
    … …  
    7272    REAL, DIMENSION(klon), INTENT(IN)        :: cdragh
    7373    REAL, DIMENSION(klon), INTENT(IN)        :: cdragm
    74     REAL, DIMENSION(klon), INTENT(IN)        :: precip_rain, precip_snow
     74    REAL, DIMENSION(klon), INTENT(IN)        :: precip_rain, precip_snow, precip_bs
    7575    REAL, DIMENSION(klon), INTENT(IN)        :: temp_air, spechum
    7676    REAL, DIMENSION(klon), INTENT(IN)        :: AcoefH, AcoefQ, BcoefH, BcoefQ
    … …  
    169169    REAL, DIMENSION(klon) :: radsol
    170170    REAL, DIMENSION(klon) :: cdragq ! Cdrag pour l'evaporation
     171    REAL, DIMENSION(klon) :: precip_totsnow
    171172    CHARACTER(len=20),PARAMETER :: modname="surf_ocean"
    172173    real rhoa(knon) ! density of moist air  (kg / m3)
    … …  
    200201    radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
    201202
     203
     204    !****************************************************************************************
     205    !Total solid precip
     206
     207    IF (ok_bs) THEN
     208       precip_totsnow=precip_snow+precip_bs
     209    ELSE
     210       precip_totsnow=precip_snow
     211    ENDIF
     212
    202213    !******************************************************************************
    203214    ! Cdragq computed from cdrag
    … …  
    227238            windsp, fder, &
    228239            itime, dtime, knon, knindex, &
    229             p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow,temp_air,spechum,&
     240            p1lay, cdragh, cdragq, cdragm, precip_rain, precip_totsnow,temp_air,spechum,&
    230241            AcoefH, AcoefQ, BcoefH, BcoefQ, &
    231242            AcoefU, AcoefV, BcoefU, BcoefV, &
    … …  
    239250       CALL ocean_slab_noice( &
    240251            itime, dtime, jour, knon, knindex, &
    241             p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, temp_air, spechum,&
     252            p1lay, cdragh, cdragq, cdragm, precip_rain, precip_totsnow, temp_air, spechum,&
    242253            AcoefH, AcoefQ, BcoefH, BcoefQ, &
    243254            AcoefU, AcoefV, BcoefU, BcoefV, &
    … …  
    250261       CALL ocean_forced_noice( &
    251262            itime, dtime, jour, knon, knindex, &
    252             p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, &
     263            p1lay, cdragh, cdragq, cdragm, precip_rain, precip_totsnow, &
    253264            temp_air, spechum, &
    254265            AcoefH, AcoefQ, BcoefH, BcoefQ, &
    … …  
    370381       call bulk_flux(tkt, tks, taur, dter, dser, t_int, s_int, ds_ns, dt_ns, &
    371382            u = windsp(:knon), t_ocean_1 = tsurf_new(:knon), s1 = sss(:knon), &
    372             rain = precip_rain(:knon) + precip_snow(:knon), &
     383            rain = precip_rain(:knon) + precip_totsnow(:knon), &
    373384            hf = - fluxsens(:knon), hlb = - fluxlat(:knon), &
    374385            rnl = - lwnet(:knon), &
Note: See TracChangeset for help on using the changeset viewer.