Ignore:
Timestamp:
Sep 4, 2023, 10:17:16 AM (15 months ago)
Author:
Laurent Fairhead
Message:

Merged with trunk revision 4586 corresponding to june 2023 testing

Location:
LMDZ6/branches/LMDZ_cdrag_LSCE
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/LMDZ_cdrag_LSCE

  • LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmd/physiq_mod.F90

    r4661 r4669  
    8484    USE atke_turbulence_ini_mod, ONLY : atke_ini
    8585    USE thermcell_ini_mod, ONLY : thermcell_ini
     86    USE blowing_snow_ini_mod, ONLY : blowing_snow_ini , qbst_bs
    8687    USE lscp_ini_mod, ONLY : lscp_ini
    8788
    … …  
    141142!!!!!!!!!!!!!!!!!!  END "USE" for CPP keys !!!!!!!!!!!!!!!!!!!!!!
    142143
     144USE physiqex_mod, ONLY : physiqex
    143145USE phys_local_var_mod, ONLY: phys_local_var_init, phys_local_var_end, &
    144146       ! [Variables internes non sauvegardees de la physique]
    145147       ! Variables locales pour effectuer les appels en serie
    146        t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,tr_seri,rneb_seri, &
     148       t_seri,q_seri,ql_seri,qs_seri,qbs_seri,u_seri,v_seri,tr_seri,rneb_seri, &
    147149       rhcl, &       
    148150       ! Dynamic tendencies (diagnostics)
    149        d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_u_dyn,d_v_dyn,d_tr_dyn,d_rneb_dyn, &
    150        d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d, &
     151       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, &
     152       d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d,d_qbs_dyn2d, &
    151153       ! Physic tendencies
    152154       d_t_con,d_q_con,d_u_con,d_v_con, &
    … …  
    165167       plul_st,plul_th, &
    166168       !
    167        d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_t_diss, &
     169       d_t_vdf,d_q_vdf, d_qbs_vdf, d_u_vdf,d_v_vdf,d_t_diss, &
    168170       d_t_vdf_x, d_t_vdf_w, &
    169171       d_q_vdf_x, d_q_vdf_w, &
    170172       d_ts, &
     173       !
     174       d_t_bs,d_q_bs,d_qbs_bs, &
    171175       !
    172176!       d_t_oli,d_u_oli,d_v_oli, &
    … …  
    216220       cldh, cldl,cldm, cldq, cldt,      &
    217221       JrNt,                             &
    218        dthmin, evap, fder, plcl, plfc,   &
    219        prw, prlw, prsw,                  &
     222       dthmin, evap, snowerosion,fder, plcl, plfc,   &
     223       prw, prlw, prsw, prbsw,                  &
    220224       s_lcl, s_pblh, s_pblt, s_therm,   &
    221225       cdragm, cdragh,                   &
    … …  
    290294       !
    291295       rneblsvol, &
     296       pfraclr,pfracld, &
     297       distcltop, &
    292298       zqsatl, zqsats, &
    293299       qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
    … …  
    305311       fsolsw, wfbils, wfbilo,  &
    306312       wfevap, wfrain, wfsnow,  & 
    307        prfl, psfl, fraca, Vprecip,  &
     313       prfl, psfl,bsfl, fraca, Vprecip,  &
    308314       zw2,  &
    309315       !
    … …  
    317323       beta_prec,  &
    318324       rneb,  &
    319        zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic
    320        !
     325       zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic, &
     326       zxfluxt,zxfluxq
     327       !
     328      USE output_physiqex_mod, ONLY: output_physiqex
     329
    321330
    322331    IMPLICIT NONE
    … …  
    461470    !======================================================================
    462471    !
    463     ! indices de traceurs eau vapeur, liquide, glace, fraction nuageuse LS (optional)
    464     INTEGER,SAVE :: ivap, iliq, isol, irneb
    465 !$OMP THREADPRIVATE(ivap, iliq, isol, irneb)
     472    ! indices de traceurs eau vapeur, liquide, glace, fraction nuageuse LS (optional), blowing snow (optional)
     473    INTEGER,SAVE :: ivap, iliq, isol, irneb, ibs
     474!$OMP THREADPRIVATE(ivap, iliq, isol, irneb, ibs)
    466475    !
    467476    !
    … …  
    834843    !XXX PB
    835844    REAL fluxq(klon,klev, nbsrf)   ! flux turbulent d'humidite
    836     !
    837     REAL zxfluxt(klon, klev)
    838     REAL zxfluxq(klon, klev)
     845    REAL fluxqbs(klon,klev, nbsrf)   ! flux turbulent de neige soufflee
     846    !
     847    !FC    REAL zxfluxt(klon, klev)
     848    !FC    REAL zxfluxq(klon, klev)
     849    REAL zxfluxqbs(klon,klev)
    839850    REAL zxfluxu(klon, klev)
    840851    REAL zxfluxv(klon, klev)
    … …  
    920931    !
    921932    ! tendance nulles
    922     REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0, dqi0
     933    REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0, dqi0, dqbs0
    923934    REAL, dimension(klon)     :: dsig0, ddens0
    924935    INTEGER, dimension(klon)  :: wkoccur1
    … …  
    944955    LOGICAL, SAVE :: ok_bug_split_th = .TRUE.
    945956    !$OMP THREADPRIVATE(ok_bug_split_th)
     957
     958    ! Logical switch to a bug : modifying directly wake_deltat  by adding
     959    ! the (w) dry adjustment tendency to wake_deltat
     960    LOGICAL, SAVE :: ok_bug_ajs_cv = .TRUE.
     961    !$OMP THREADPRIVATE(ok_bug_ajs_cv)
    946962
    947963    !
    … …  
    10611077    REAL ztsol(klon)
    10621078    REAL q2m(klon,nbsrf)  ! humidite a 2m
    1063 
     1079    REAL fsnowerosion(klon,nbsrf) ! blowing snow flux at surface
     1080    REAL qbsfra  ! blowing snow fraction
    10641081    !IM: t2m, q2m, ustar, u10m, v10m et t2mincels, t2maxcels
    10651082    CHARACTER*40 t2mincels, t2maxcels       !t2m min., t2m max
    … …  
    11361153    !IM 100106 BEG : pouvoir sortir les ctes de la physique
    11371154    include "conema3.h"
    1138     include "fisrtilp.h"
    11391155    include "nuage.h"
    11401156    include "compbl.h"
    … …  
    11441160    ! Declarations pour Simulateur COSP
    11451161    !============================================================
     1162    ! AI 10-22
     1163#ifdef CPP_COSP
     1164    include "ini_COSP.h"
     1165#endif
    11461166    real :: mr_ozone(klon,klev), phicosp(klon,klev)
    11471167
    … …  
    12181238    REAL pi
    12191239
     1240
     1241    !======================================================================!
     1242    ! Bifurcation vers un nouveau moniteur physique pour experimenter      !
     1243    ! des solutions et préparer le couplage avec la physique de MesoNH     !
     1244    ! 14 mai 2023                                                          !
     1245    !======================================================================!
     1246    if (debut) then                                                        !
     1247       iflag_physiq=0
     1248       call getin_p('iflag_physiq', iflag_physiq)                          !
     1249    endif                                                                  !
     1250    if ( iflag_physiq == 2 ) then                                          !
     1251       call physiqex (nlon,nlev, &                                         !
     1252       debut,lafin,pdtphys_, &                                             !
     1253       paprs,pplay,pphi,pphis,presnivs, &                                  !
     1254       u,v,rot,t,qx, &                                                     !
     1255       flxmass_w, &                                                        !
     1256       d_u, d_v, d_t, d_qx, d_ps)                                          !
     1257       return                                                              !
     1258    endif                                                                  !
     1259    !======================================================================!
     1260
     1261
    12201262    pi = 4. * ATAN(1.)
    12211263
    … …  
    12441286    ! Utilise notamment en 1D mais peut etre active egalement en 3D
    12451287    ! en imposant la valeur de igout.
    1246     !======================================================================d
     1288    !======================================================================
    12471289    IF (prt_level.ge.1) THEN
    12481290       igout=klon/2+1/klon
    … …  
    12751317       isol = strIdx(tracers(:)%name, addPhase('H2O', 's'))
    12761318       irneb= strIdx(tracers(:)%name, addPhase('H2O', 'r'))
     1319       ibs  = strIdx(tracers(:)%name, addPhase('H2O', 'b'))
    12771320       CALL init_etat0_limit_unstruct
    12781321       IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
    … …  
    13241367       ENDIF
    13251368
    1326        IF (ok_ice_sursat.AND.(nqo.NE.4)) THEN
     1369       IF (ok_ice_sursat.AND.(nqo.LT.4)) THEN
    13271370          WRITE (lunout, *) ' ok_ice_sursat=y requires 4 H2O tracers ', &
    13281371               '(H2O_g, H2O_l, H2O_s, H2O_r) but nqo=', nqo, '. Might as well stop here.'
    … …  
    13421385          CALL abort_physic(modname,abort_message,1)
    13431386       ENDIF
     1387
     1388        IF (ok_bs) THEN
     1389         IF ((ok_ice_sursat.AND.nqo .LT.5).OR.(.NOT.ok_ice_sursat.AND.nqo.LT.4)) THEN
     1390             WRITE (lunout, *) 'activation of blowing snow needs a specific H2O tracer', &
     1391                               'but nqo=', nqo
     1392             abort_message='see above'
     1393             CALL abort_physic(modname,abort_message, 1)
     1394         ENDIF
     1395        ENDIF
    13441396
    13451397       Ncvpaseq1 = 0
    … …  
    14041456       CALL getin_p('ok_bug_cv_trac',ok_bug_cv_trac)
    14051457       CALL getin_p('ok_bug_split_th',ok_bug_split_th)
     1458       CALL getin_p('ok_bug_ajs_cv',ok_bug_ajs_cv)
    14061459       fl_ebil = 0 ! by default, conservation diagnostics are desactivated
    14071460       CALL getin_p('fl_ebil',fl_ebil)
    … …  
    16761729      IF (.NOT. create_etat0_limit) CALL init_readaerosolstrato(flag_aerosol_strat)  !! initialise aero strato from file for XIOS interpolation (unstructured_grid)
    16771730
     1731      if (ok_cosp) then
    16781732#ifdef CPP_COSP
    1679       IF (ok_cosp) THEN
    1680 !           DO k = 1, klev
    1681 !             DO i = 1, klon
    1682 !               phicosp(i,k) = pphi(i,k) + pphis(i)
    1683 !             ENDDO
    1684 !           ENDDO
     1733        ! A.I : Initialisations pour le 1er passage a Cosp
     1734        CALL ini_COSP(ref_liq_cosp0,ref_ice_cosp0,pctsrf_cosp0,zu10m_cosp0,zv10m_cosp0, &
     1735               zxtsol_cosp0,zx_rh_cosp0,cldfra_cosp0,rnebcon_cosp0,flwc_cosp0, &
     1736               fiwc_cosp0,prfl_cosp0,psfl_cosp0,pmflxr_cosp0,pmflxs_cosp0, &
     1737               mr_ozone_cosp0,cldtau_cosp0,cldemi_cosp0,JrNt_cosp0)
     1738
    16851739        CALL phys_cosp(itap,phys_tstep,freq_cosp, &
    16861740               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
    … …  
    16951749               pmflxr(:,1:klev),pmflxs(:,1:klev), &
    16961750               mr_ozone,cldtau, cldemi)
    1697       ENDIF
    16981751#endif
    16991752
    17001753#ifdef CPP_COSP2
    1701         IF (ok_cosp) THEN
    1702 !           DO k = 1, klev
    1703 !             DO i = 1, klon
    1704 !               phicosp(i,k) = pphi(i,k) + pphis(i)
    1705 !             ENDDO
    1706 !           ENDDO
    17071754          CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
    17081755               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
    … …  
    17171764               pmflxr(:,1:klev),pmflxs(:,1:klev), &
    17181765               mr_ozone,cldtau, cldemi)
    1719        ENDIF
    17201766#endif
    17211767
    17221768#ifdef CPP_COSPV2
    1723         IF (ok_cosp) THEN
    1724            DO k = 1, klev
    1725              DO i = 1, klon
    1726                phicosp(i,k) = pphi(i,k) + pphis(i)
    1727              ENDDO
    1728            ENDDO
    17291769          CALL lmdz_cosp_interface(itap,phys_tstep,freq_cosp, &
    17301770               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
    … …  
    17391779               pmflxr(:,1:klev),pmflxs(:,1:klev), &
    17401780               mr_ozone,cldtau, cldemi)
    1741        ENDIF
    17421781#endif
     1782      ENDIF
    17431783
    17441784       !
    … …  
    17561796       CALL thermcell_ini(iflag_thermals,prt_level,tau_thermals,lunout, &
    17571797   &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
    1758        IF (ok_new_lscp) then
    1759            CALL lscp_ini(pdtphys,ok_ice_sursat)
    1760        endif
     1798       CALL lscp_ini(pdtphys,ok_ice_sursat, RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG)
     1799       CALL blowing_snow_ini(prt_level,lunout, &
     1800                             RCPD, RLSTT, RLVTT, RLMLT, &
     1801                             RVTMP2, RTT,RD,RG)
     1802
    17611803!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    17621804
    … …  
    17841826       CALL phys_output_write(itap, pdtphys, paprs, pphis,                    &
    17851827                              pplay, lmax_th, aerosol_couple,                 &
    1786                               ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, ok_sync,&
     1828                              ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, ibs,  ok_sync,&
    17871829                              ptconv, read_climoz, clevSTD,                   &
    17881830                              ptconvth, d_u, d_t, qx, d_qx, zmasse,           &
    … …  
    22552297    dql0(:,:)=0.
    22562298    dqi0(:,:)=0.
     2299    dqbs0(:,:)=0.
    22572300    dsig0(:) = 0.
    22582301    ddens0(:) = 0.
    … …  
    23072350          q_seri(i,k)  = qx(i,k,ivap)
    23082351          ql_seri(i,k) = qx(i,k,iliq)
     2352          qbs_seri(i,k) = 0.
    23092353          !CR: ATTENTION, on rajoute la variable glace
    23102354          IF (nqo.EQ.2) THEN             !--vapour and liquid only
    … …  
    23142358             qs_seri(i,k) = qx(i,k,isol)
    23152359             rneb_seri(i,k) = 0.
    2316           ELSE IF (nqo.EQ.4) THEN        !--vapour, liquid, ice and rneb
     2360          ELSE IF (nqo.GE.4) THEN        !--vapour, liquid, ice and rneb and blowing snow
    23172361             qs_seri(i,k) = qx(i,k,isol)
     2362             IF (ok_ice_sursat) THEN
    23182363             rneb_seri(i,k) = qx(i,k,irneb)
     2364             ENDIF
     2365             IF (ok_bs) THEN
     2366             qbs_seri(i,k)= qx(i,k,ibs)
     2367             ENDIF
     2368
    23192369          ENDIF
     2370
     2371
    23202372       ENDDO
    23212373    ENDDO
    … …  
    23262378    qql1(:)=0.0
    23272379    DO k = 1, klev
    2328       qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
     2380      qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k)+qbs_seri(:,k))*zmasse(:,k)
    23292381    ENDDO
    23302382    ENDIF
    … …  
    23822434       d_ql_dyn(:,:) = (ql_seri(:,:)-ql_ancien(:,:))/phys_tstep
    23832435       d_qs_dyn(:,:) = (qs_seri(:,:)-qs_ancien(:,:))/phys_tstep
     2436       d_qbs_dyn(:,:) = (qbs_seri(:,:)-qbs_ancien(:,:))/phys_tstep
    23842437       CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
    23852438       d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/phys_tstep
    … …  
    23882441       CALL water_int(klon,klev,qs_seri,zmasse,zx_tmp_fi2d)
    23892442       d_qs_dyn2d(:)=(zx_tmp_fi2d(:)-prsw_ancien(:))/phys_tstep
     2443       CALL water_int(klon,klev,qbs_seri,zmasse,zx_tmp_fi2d)
     2444       d_qbs_dyn2d(:)=(zx_tmp_fi2d(:)-prbsw_ancien(:))/phys_tstep
    23902445       ! !! RomP >>>   td dyn traceur
    23912446       IF (nqtot > nqo) d_tr_dyn(:,:,:)=(tr_seri(:,:,:)-tr_ancien(:,:,:))/phys_tstep
    … …  
    24032458       d_ql_dyn2d(:) = 0.0
    24042459       d_qs_dyn2d(:) = 0.0
     2460       d_qbs_dyn2d(:)= 0.0
    24052461       ! !! RomP >>>   td dyn traceur
    24062462       IF (nqtot > nqo) d_tr_dyn(:,:,:)= 0.0
    24072463       ! !! RomP <<<
    24082464       d_rneb_dyn(:,:)=0.0
     2465       d_qbs_dyn(:,:)=0.0
    24092466       ancien_ok = .TRUE.
    24102467    ENDIF
    … …  
    25222579
    25232580     CALL add_phys_tend &
    2524             (du0,dv0,d_t_eva,d_q_eva,d_ql_eva,d_qi_eva,paprs,&
     2581            (du0,dv0,d_t_eva,d_q_eva,d_ql_eva,d_qi_eva,dqbs0,paprs,&
    25252582               'eva',abortphy,flag_inhib_tend,itap,0)
    25262583    CALL prt_enerbil('eva',itap)
    … …  
    26982755            longitude_deg, latitude_deg, rugoro,  zrmu0,      &
    26992756            sollwdown,    cldt,      &
    2700             rain_fall, snow_fall, solsw,   solswfdiff, sollw,     &
     2757            rain_fall, snow_fall, bs_fall, solsw,   solswfdiff, sollw,     &
    27012758            gustiness,                                &
    2702             t_seri,    q_seri,    u_seri,  v_seri,    &
     2759            t_seri,    q_seri,   qbs_seri, u_seri,  v_seri,    &
    27032760                                !nrlmd+jyg<
    27042761            wake_deltat, wake_deltaq, wake_cstar, wake_s, &
    … …  
    27112768                                !albedo SB >>>
    27122769                                ! albsol1,   albsol2,   sens,    evap,      &
    2713             albsol_dir,   albsol_dif,   sens,    evap,   & 
     2770            albsol_dir,   albsol_dif,   sens,    evap, snowerosion, & 
    27142771                                !albedo SB <<<
    27152772            albsol3_lic,runoff,   snowhgt,   qsnow, to_ice, sissnow, &
    27162773            zxtsol,    zxfluxlat, zt2m,    qsat2m,  zn2mout, &
    2717             d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf, d_t_diss, &
     2774            d_t_vdf,   d_q_vdf, d_qbs_vdf,  d_u_vdf, d_v_vdf, d_t_diss, &
    27182775                                !nrlmd<
    27192776                                !jyg<
    … …  
    27412798            fluxt,   fluxu,  fluxv, &
    27422799            dsens,     devap,     zxsnow, &
    2743             zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke, &
     2800            zxfluxt,   zxfluxq,  zxfluxqbs,  q2m, fluxq, fluxqbs, pbl_tke, &
    27442801                                !nrlmd+jyg<
    27452802            wake_delta_pbl_TKE, &
    … …  
    27662823       IF (klon_glo==1) THEN
    27672824          CALL add_pbl_tend &
    2768                (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
     2825               (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,d_qbs_vdf,paprs,&
    27692826               'vdf',abortphy,flag_inhib_tend,itap)
    27702827       ELSE
    27712828          CALL add_phys_tend &
    2772                (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
     2829               (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,d_qbs_vdf,paprs,&
    27732830               'vdf',abortphy,flag_inhib_tend,itap,0)
    27742831       ENDIF
    27752832       CALL prt_enerbil('vdf',itap)
     2833
    27762834       !--------------------------------------------------------------------
    27772835
    … …  
    28242882
    28252883    ENDIF
     2884
     2885    ! ==================================================================
     2886    ! Blowing snow sublimation and sedimentation
     2887
     2888    d_t_bs(:,:)=0.
     2889    d_q_bs(:,:)=0.
     2890    d_qbs_bs(:,:)=0.
     2891    bsfl(:,:)=0.
     2892    bs_fall(:)=0.
     2893    IF (ok_bs) THEN
     2894
     2895     CALL call_blowing_snow_sublim_sedim(klon,klev,phys_tstep,t_seri,q_seri,qbs_seri,pplay,paprs, &
     2896                                        d_t_bs,d_q_bs,d_qbs_bs,bsfl,bs_fall)
     2897
     2898     CALL add_phys_tend &
     2899               (du0,dv0,d_t_bs,d_q_bs,dql0,dqi0,d_qbs_bs,paprs,&
     2900               'bs',abortphy,flag_inhib_tend,itap,0)
     2901
     2902    ENDIF
     2903
    28262904    ! =================================================================== c
    28272905    !   Calcul de Qsat
    … …  
    29943072                ENDDO
    29953073             ENDDO
    2996              IF (iflag_adjwk == 2) THEN
     3074             IF (iflag_adjwk == 2 .AND. OK_bug_ajs_cv) THEN
    29973075               CALL add_wake_tend &
    29983076                 (d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, ddens0, ddens0, wkoccur1, 'ajs_cv', abortphy)
    2999              ENDIF  ! (iflag_adjwk == 2)
     3077             ENDIF  ! (iflag_adjwk == 2 .AND. OK_bug_ajs_cv)
    30003078          ENDIF  ! (iflag_adjwk >= 1)
    30013079       ENDIF ! (iflag_wake>=1)
    … …  
    32483326!!
    32493327!!
    3250     CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, paprs, &
     3328    CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, dqbs0, paprs, &
    32513329         'convection',abortphy,flag_inhib_tend,itap,0)
    32523330    CALL prt_enerbil('convection',itap)
    … …  
    33793457       !-----------------------------------------------------------------------
    33803458       ! ajout des tendances des poches froides
    3381        CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake', &
     3459       CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,dqbs0,paprs,'wake', &
    33823460            abortphy,flag_inhib_tend,itap,0)
    33833461       CALL prt_enerbil('wake',itap)
    … …  
    35373615          !
    35383616          CALL add_phys_tend(d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs,  &
    3539                              dql0,dqi0,paprs,'thermals', abortphy,flag_inhib_tend,itap,0)
     3617                             dql0,dqi0,dqbs0,paprs,'thermals', abortphy,flag_inhib_tend,itap,0)
    35403618          CALL prt_enerbil('thermals',itap)
    35413619          !
    … …  
    35993677          !--------------------------------------------------------------------
    36003678          ! ajout des tendances de l'ajustement sec ou des thermiques
    3601           CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs, &
     3679          CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,dqbs0,paprs, &
    36023680               'ajsb',abortphy,flag_inhib_tend,itap,0)
    36033681          CALL prt_enerbil('ajsb',itap)
    … …  
    36333711         ptconv,ptconvth,clwcon0th, rnebcon0th,     &
    36343712         paprs,pplay,t_seri,q_seri, qtc_cv, sigt_cv, zqsat, &
    3635          pbl_tke(:,:,is_ave),tke_dissip_ave,l_mix_ave,wprime_ave,t2m,q2m,fm_therm, &
     3713         pbl_tke(:,:,is_ave),tke_dissip_ave,l_mix_ave,wprime_ave,t2m,q2m,fm_therm,cell_area, &
    36363714         ratqs,ratqsc,ratqs_inter)
    36373715
    … …  
    36573735         t_seri, q_seri,ptconv,ratqs, &
    36583736         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, rneb_seri, &
     3737         pfraclr,pfracld, &
    36593738         radocond, picefra, rain_lsc, snow_lsc, &
    36603739         frac_impa, frac_nucl, beta_prec_fisrt, &
    36613740         prfl, psfl, rhcl,  &
    36623741         zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
    3663          iflag_ice_thermo, ok_ice_sursat, zqsatl, zqsats, &
     3742         iflag_ice_thermo, ok_ice_sursat, zqsatl, zqsats, distcltop, &
    36643743         qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
    36653744         Tcontr, qcontr, qcontr2, fcontrN, fcontrP )
    … …  
    36883767!    write(*,9000) "rcpv","rcw",rcpv,rcw,rcs,t_seri(1,1)
    36893768!-JLD
    3690     CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs, &
     3769    CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,dqbs0,paprs, &
    36913770         'lsc',abortphy,flag_inhib_tend,itap,0)
    36923771    CALL prt_enerbil('lsc',itap)
    … …  
    37113790    ENDIF
    37123791
    3713     !---------------------------------------------------------------------------
     3792
     3793!---------------------------------------------------------------------------
    37143794    DO k = 1, klev
    37153795       DO i = 1, klon
    37163796          cldfra(i,k) = rneb(i,k)
    37173797          !CR: a quoi ca sert? Faut-il ajouter qs_seri?
     3798          !EV: en effet etrange, j'ajouterais aussi qs_seri
     3799          !    plus largement, je nettoierais (enleverrais) ces lignes
    37183800          IF (.NOT.new_oliq) radocond(i,k) = ql_seri(i,k)
    37193801       ENDDO
    37203802    ENDDO
     3803
     3804
     3805    ! Option to activate the radiative effect of blowing snow (ok_rad_bs)
     3806    ! makes sense only if the new large scale condensation scheme is active
     3807    ! with the ok_icefra_lscp flag active as well
     3808
     3809    IF (ok_bs .AND. ok_rad_bs) THEN
     3810       IF (ok_new_lscp .AND. ok_icefra_lscp) THEN
     3811           DO k=1,klev
     3812             DO i=1,klon
     3813                radocond(i,k)=radocond(i,k)+qbs_seri(i,k)
     3814                picefra(i,k)=(radocond(i,k)*picefra(i,k)+qbs_seri(i,k))/(radocond(i,k))
     3815                qbsfra=min(qbs_seri(i,k)/qbst_bs,1.0)
     3816                cldfra(i,k)=max(cldfra(i,k),qbsfra)
     3817             ENDDO
     3818           ENDDO
     3819       ELSE
     3820          WRITE(lunout,*)"PAY ATTENTION, you try to activate the radiative effect of blowing snow"
     3821          WRITE(lunout,*)"with ok_new_lscp=false and/or ok_icefra_lscp=false"
     3822          abort_message='inconsistency in cloud phase for blowing snow'
     3823          CALL abort_physic(modname,abort_message,1)
     3824       ENDIF
     3825
     3826    ENDIF
     3827
    37213828    IF (check) THEN
    37223829       za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
    … …  
    42064313               flwp, fiwp, flwc, fiwc, &
    42074314               mass_solu_aero, mass_solu_aero_pi, &
    4208                cldtaupi, latitude_deg, re, fl, ref_liq, ref_ice, &
     4315               cldtaupi, latitude_deg, distcltop, re, fl, ref_liq, ref_ice, &
    42094316               ref_liq_pi, ref_ice_pi)
    42104317       ELSE
    … …  
    42144321               ok_aie, &
    42154322               mass_solu_aero, mass_solu_aero_pi, &
    4216                bl95_b0, bl95_b1, &
     4323               bl95_b0, bl95_b1, distcltop, &
    42174324               cldtaupi, re, fl)
    42184325       ENDIF
    … …  
    45194626    ENDDO
    45204627
    4521     CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy,flag_inhib_tend,itap,0)
     4628    CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,dqbs0,paprs,'SW',abortphy,flag_inhib_tend,itap,0)
    45224629    CALL prt_enerbil('SW',itap)
    4523     CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy,flag_inhib_tend,itap,0)
     4630    CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,dqbs0,paprs,'LW',abortphy,flag_inhib_tend,itap,0)
    45244631    CALL prt_enerbil('LW',itap)
    45254632
    … …  
    45964703       !-----------------------------------------------------------------------
    45974704       ! ajout des tendances de la trainee de l'orographie
    4598        CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro', &
     4705       CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,dqbs0,paprs,'oro', &
    45994706            abortphy,flag_inhib_tend,itap,0)
    46004707       CALL prt_enerbil('oro',itap)
    … …  
    46474754
    46484755       ! ajout des tendances de la portance de l'orographie
    4649        CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, paprs, &
     4756       CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, dqbs0,paprs, &
    46504757            'lif', abortphy,flag_inhib_tend,itap,0)
    46514758       CALL prt_enerbil('lif',itap)
    … …  
    46724779       d_t_hin(:, :)=0.
    46734780       CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
    4674             dqi0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0)
     4781            dqi0, dqbs0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0)
    46754782       CALL prt_enerbil('hin',itap)
    46764783    ENDIF
    … …  
    46904797       ENDDO
    46914798
    4692        CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, &
     4799       CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, dqbs0, &
    46934800            paprs, 'front_gwd_rando', abortphy,flag_inhib_tend,itap,0)
    46944801       CALL prt_enerbil('front_gwd_rando',itap)
    … …  
    46994806            rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
    47004807            du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
    4701        CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, &
     4808       CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, dqbs0, &
    47024809            paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend,itap,0)
    47034810       CALL prt_enerbil('flott_gwd_rando',itap)
    … …  
    47514858       ! ajout de la tendance d'humidite due au methane
    47524859       d_q_ch4_dtime(:,:) = d_q_ch4(:,:)*phys_tstep
    4753        CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, paprs, &
     4860       CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, dqbs0, paprs, &
    47544861            'q_ch4', abortphy,flag_inhib_tend,itap,0)
    47554862       d_q_ch4(:,:) = d_q_ch4_dtime(:,:)/phys_tstep
    … …  
    50565163
    50575164     CALL add_phys_tend &
    5058             (du0,dv0,dt0,d_q_rep,d_ql_rep,d_qi_rep,paprs,&
     5165            (du0,dv0,dt0,d_q_rep,d_ql_rep,d_qi_rep,dqbs0,paprs,&
    50595166             'rep',abortphy,flag_inhib_tend,itap,0)
    50605167        IF (abortphy==1) Print*,'ERROR ABORT REP'
    … …  
    51335240    !   prlw = colonne eau liquide
    51345241    !   prlw = colonne eau solide
     5242    !   prbsw = colonne neige soufflee
    51355243    prw(:) = 0.
    51365244    prlw(:) = 0.
    51375245    prsw(:) = 0.
     5246    prbsw(:) = 0.
    51385247    DO k = 1, klev
    51395248       prw(:)  = prw(:)  + q_seri(:,k)*zmasse(:,k)
    51405249       prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
    51415250       prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
     5251       prbsw(:)= prbsw(:) + qbs_seri(:,k)*zmasse(:,k)
    51425252    ENDDO
    51435253    !
    … …  
    51985308          ENDIF
    51995309          !--ice_sursat: nqo=4, on ajoute rneb
    5200           IF (nqo == 4) THEN
     5310          IF (nqo.ge.4 .and. ok_ice_sursat) THEN
    52015311             d_qx(i,k,irneb) = ( rneb_seri(i,k) - qx(i,k,irneb) ) / phys_tstep
    52025312          ENDIF
     5313
     5314           IF (nqo.ge.4 .and. ok_bs) THEN
     5315             d_qx(i,k,ibs) = ( qbs_seri(i,k) - qx(i,k,ibs) ) / phys_tstep
     5316          ENDIF
     5317
    52035318       ENDDO
    52045319    ENDDO
    … …  
    52475362    ql_ancien(:,:) = ql_seri(:,:)
    52485363    qs_ancien(:,:) = qs_seri(:,:)
     5364    qbs_ancien(:,:) = qbs_seri(:,:)
    52495365    rneb_ancien(:,:) = rneb_seri(:,:)
    52505366    CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
    52515367    CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
    52525368    CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
     5369    CALL water_int(klon,klev,qbs_ancien,zmasse,prbsw_ancien)
    52535370    ! !! RomP >>>
    52545371    IF (nqtot > nqo) tr_ancien(:,:,:) = tr_seri(:,:,:)
    … …  
    53755492    CALL phys_output_write(itap, pdtphys, paprs, pphis,  &
    53765493         pplay, lmax_th, aerosol_couple,                 &
    5377          ok_ade, ok_aie, ok_volcan, ivap, iliq, isol,    &
     5494         ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, ibs,   &
    53785495         ok_sync, ptconv, read_climoz, clevSTD,          &
    53795496         ptconvth, d_u, d_t, qx, d_qx, zmasse,           &
    … …  
    53865503
    53875504#endif
     5505    ! Petit appelle de sorties pour accompagner le travail sur phyex
     5506    if ( iflag_physiq == 1 ) then
     5507        call output_physiqex(debut,jD_eq,pdtphys,presnivs,paprs,u,v,t,qx,cldfra,0.*t,0.*t,0.*t,pbl_tke,theta)
     5508    endif
    53885509
    53895510    !====================================================================
    … …  
    54065527    ! Disabling calls to the prt_alerte function
    54075528    alert_first_call = .FALSE.
     5529
    54085530   
    54095531    IF (lafin) THEN
Note: See TracChangeset for help on using the changeset viewer.