Ignore:
Timestamp:
Nov 19, 2021, 4:58:59 PM (3 years ago)
Author:
lguez
Message:

Sync latest trunk changes to Ocean_skin

Location:
LMDZ6/branches/Ocean_skin
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • LMDZ6/branches/Ocean_skin

  • LMDZ6/branches/Ocean_skin/libf/phylmd/pbl_surface_mod.F90

    r3798 r4013  
    2323  USE climb_wind_mod,      ONLY : climb_wind_down, climb_wind_up
    2424  USE coef_diff_turb_mod,  ONLY : coef_diff_turb
    25   USE wx_pbl_mod,          ONLY : wx_pbl_init, wx_pbl_final, &
    26 !!                                  wx_pbl_fuse_no_dts, wx_pbl_split_no_dts, &
    27 !!                                  wx_pbl_fuse, wx_pbl_split
    28                                   wx_pbl0_fuse, wx_pbl0_split
     25  USE ioipsl_getin_p_mod,  ONLY : getin_p
     26  USE cdrag_mod
     27  USE stdlevvar_mod
     28  USE wx_pbl_var_mod,      ONLY : wx_pbl_init, wx_pbl_final, &
     29                                  wx_pbl_prelim_0, wx_pbl_prelim_beta
     30  USE wx_pbl_mod,          ONLY : wx_pbl0_merge, wx_pbl_split, wx_pbl_dts_merge, &
     31                                  wx_pbl_check, wx_pbl_dts_check, wx_evappot
    2932  use config_ocean_skin_m, only: activate_ocean_skin
    3033
     
    3437  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: fder   ! flux drift
    3538  !$OMP THREADPRIVATE(fder)
    36   REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: snow   ! snow at surface
     39  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE    :: snow   ! snow at surface
    3740  !$OMP THREADPRIVATE(snow)
    3841  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: qsurf  ! humidity at surface
    3942  !$OMP THREADPRIVATE(qsurf)
    40   REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: ftsoil ! soil temperature
     43  REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE          :: ftsoil ! soil temperature
    4144  !$OMP THREADPRIVATE(ftsoil)
     45  REAL, ALLOCATABLE, DIMENSION(:), SAVE              :: ydTs0, ydqs0 
     46                                                     ! nul forced temperature and humidity differences
     47  !$OMP THREADPRIVATE(ydTs0, ydqs0)
    4248
    4349  INTEGER, SAVE :: iflag_pbl_surface_t2m_bug
    4450  !$OMP THREADPRIVATE(iflag_pbl_surface_t2m_bug)
     51  INTEGER, SAVE :: iflag_new_t2mq2m
     52  !$OMP THREADPRIVATE(iflag_new_t2mq2m)
     53
    4554!FC
    4655!  integer, save :: iflag_frein
     
    93102    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
    94103
     104    ALLOCATE(ydTs0(klon), stat=ierr)
     105    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
     106
     107    ALLOCATE(ydqs0(klon), stat=ierr)
     108    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
     109
    95110    fder(:)       = fder_rst(:)
    96111    snow(:,:)     = snow_rst(:,:)
    97112    qsurf(:,:)    = qsurf_rst(:,:)
    98113    ftsoil(:,:,:) = ftsoil_rst(:,:,:)
     114    ydTs0(:) = 0.
     115    ydqs0(:) = 0.
    99116
    100117!****************************************************************************************
     
    142159    iflag_pbl_surface_t2m_bug=0
    143160    CALL getin_p('iflag_pbl_surface_t2m_bug',iflag_pbl_surface_t2m_bug)
     161    WRITE(lunout,*) 'iflag_pbl_surface_t2m_bug=',iflag_pbl_surface_t2m_bug
    144162!FC
    145163!    iflag_frein = 0
     
    164182       debut,     lafin,                              &
    165183       rlon,      rlat,      rugoro,   rmu0,          &
    166        zsig,      lwdown_m,  pphi,     cldt,          &
     184       lwdown_m,  cldt,          &
    167185       rain_f,    snow_f,    solsw_m,  solswfdiff_m, sollw_m,       &
    168186       gustiness,                                     &
     
    176194       ts,SFRWL,   alb_dir, alb_dif,ustar, u10m, v10m,wstar, &
    177195       cdragh,    cdragm,   zu1,    zv1,              &
     196!jyg<   (26/09/2019)
     197       beta, &
     198!>jyg
    178199       alb_dir_m,    alb_dif_m,  zxsens,   zxevap,    &
    179200       alb3_lic,  runoff,    snowhgt,   qsnow,     to_ice,    sissnow,  &
    180        zxtsol,    zxfluxlat, zt2m,     qsat2m,       &
     201       zxtsol,    zxfluxlat, zt2m,     qsat2m, zn2mout, &
    181202       d_t,       d_q,       d_u,      d_v, d_t_diss, &
    182203!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
     
    199220       s_therm,   s_trmb1,   s_trmb2,  s_trmb3,       &
    200221       zustar,zu10m,  zv10m,    fder_print,    &
    201        zxqsurf,   rh2m,      zxfluxu,  zxfluxv,       &
     222       zxqsurf, delta_qsurf,                       &
     223       rh2m,      zxfluxu,  zxfluxv,               &
    202224       z0m, z0h,   agesno,  sollw,    solsw,         &
    203225       d_ts,      evap,    fluxlat,  t2m,           &
     
    255277! z0m, z0h ----input-R- longeur de rugosite (en m)
    256278! Martin
    257 ! zsig-----input-R- slope
    258279! cldt-----input-R- total cloud fraction
    259 ! pphi-----input-R- geopotentiel de chaque couche (g z) (reference sol)
    260280! Martin
    261281!
     
    293313    USE print_control_mod,  ONLY : prt_level,lunout
    294314    USE ioipsl_getin_p_mod, ONLY : getin_p
    295     use phys_state_var_mod, only: ds_ns, dt_ns, delta_sst, delta_sal
     315    use phys_state_var_mod, only: ds_ns, dt_ns, delta_sst, delta_sal, zsig, zmea
    296316    use phys_output_var_mod, only: dter, dser, tkt, tks, taur, sss
    297317#ifdef CPP_XIOS
     
    300320    use netcdf, only: missing_val => nf90_fill_real
    301321#endif
     322
     323     
     324
    302325
    303326    IMPLICIT NONE
     
    337360    REAL, DIMENSION(klon, nbsrf), INTENT(IN)        :: pctsrf  ! sub-surface fraction
    338361! Martin
    339     REAL, DIMENSION(klon),        INTENT(IN)        :: zsig    ! slope
    340362    REAL, DIMENSION(klon),        INTENT(IN)        :: lwdown_m ! downward longwave radiation at mean s   
    341363    REAL, DIMENSION(klon),        INTENT(IN)        :: gustiness ! gustiness
    342364
    343365    REAL, DIMENSION(klon),        INTENT(IN)        :: cldt    ! total cloud fraction
    344     REAL, DIMENSION(klon,klev),   INTENT(IN)        :: pphi    ! geopotential (m2/s2)
    345 ! Martin
    346366
    347367!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
     
    359379! Input/Output variables
    360380!****************************************************************************************
     381!jyg<
     382    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: beta    ! Aridity factor
     383!>jyg
    361384    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: ts      ! temperature at surface (K)
    362385    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: delta_tsurf !surface temperature difference between
     
    404427    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxfluxlat  ! latent flux, mean for each grid point
    405428    REAL, DIMENSION(klon),        INTENT(OUT)       :: zt2m       ! temperature at 2m, mean for each grid point
     429    INTEGER, DIMENSION(klon, 6),  INTENT(OUT)       :: zn2mout    ! number of times the 2m temperature is out of the [tsol,temp]
    406430    REAL, DIMENSION(klon),        INTENT(OUT)       :: qsat2m
    407431    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_t        ! change in temperature
     
    460484    REAL, DIMENSION(klon),        INTENT(OUT)       :: fder_print ! fder for printing (=fder(i) + dflux_t(i) + dflux_q(i))
    461485    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxqsurf    ! humidity at surface, mean for each grid point
     486    REAL, DIMENSION(klon),        INTENT(OUT)       :: delta_qsurf! humidity difference at surface, mean for each grid point
    462487    REAL, DIMENSION(klon),        INTENT(OUT)       :: rh2m       ! relative humidity at 2m
    463488    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zxfluxu    ! u wind tension, mean for each grid point
     
    494519
    495520! Martin
    496 ! sisvat
     521! inlandsis
    497522    REAL, DIMENSION(klon),       INTENT(OUT)        :: qsnow      ! snow water content
    498523    REAL, DIMENSION(klon),       INTENT(OUT)        :: snowhgt    ! snow height
     
    521546    INTEGER                            :: n
    522547! << PC
    523     INTEGER                            :: iflag_split
     548    INTEGER                            :: iflag_split, iflag_split_ref
    524549    INTEGER                            :: i, k, nsrf
    525550    INTEGER                            :: knon, j
     
    532557    REAL, DIMENSION(klon)              :: r_co2_ppm     ! taux CO2 atmosphere
    533558    REAL, DIMENSION(klon)              :: yts, yz0m, yz0h, ypct
     559    REAL, DIMENSION(klon)              :: yz0h_old
    534560!albedo SB >>>
    535561    REAL, DIMENSION(klon)              :: yalb,yalb_vis
    536562!albedo SB <<<
    537563    REAL, DIMENSION(klon)              :: yt1, yq1, yu1, yv1
     564    REAL, DIMENSION(klon)              :: yqa
    538565    REAL, DIMENSION(klon)              :: ysnow, yqsurf, yagesno, yqsol
    539566    REAL, DIMENSION(klon)              :: yrain_f, ysnow_f
     
    547574    REAL, DIMENSION(klon)              :: y_flux_u1, y_flux_v1
    548575    REAL, DIMENSION(klon)              :: yt2m, yq2m, yu10m
     576    INTEGER, DIMENSION(klon, nbsrf, 6) :: yn2mout, yn2mout_x, yn2mout_w
     577    INTEGER, DIMENSION(klon, nbsrf, 6) :: n2mout, n2mout_x, n2mout_w
    549578    REAL, DIMENSION(klon)              :: yustar
    550579    REAL, DIMENSION(klon)              :: ywstar
     
    567596    REAL, DIMENSION(klon)              :: yz0h_oupas
    568597    REAL, DIMENSION(klon)              :: yfluxsens
     598    REAL, DIMENSION(klon)              :: AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0
    569599    REAL, DIMENSION(klon)              :: AcoefH, AcoefQ, BcoefH, BcoefQ
    570600    REAL, DIMENSION(klon)              :: AcoefU, AcoefV, BcoefU, BcoefV
    571601    REAL, DIMENSION(klon)              :: ypsref
    572     REAL, DIMENSION(klon)              :: yevap, ytsurf_new, yalb3_new
     602    REAL, DIMENSION(klon)              :: yevap, yevap_pot, ytsurf_new, yalb3_new
    573603!albedo SB >>>
    574604    REAL, DIMENSION(klon,nsw)          :: yalb_dir_new, yalb_dif_new
     
    582612    REAL, DIMENSION(klon,klev)         :: y_flux_u, y_flux_v
    583613    REAL, DIMENSION(klon,klev)         :: ycoefh, ycoefm,ycoefq
    584     REAL, DIMENSION(klon)              :: ycdragh, ycdragm
     614    REAL, DIMENSION(klon)              :: ycdragh, ycdragq, ycdragm
    585615    REAL, DIMENSION(klon,klev)         :: yu, yv
    586616    REAL, DIMENSION(klon,klev)         :: yt, yq
     
    614644    REAL, DIMENSION(klon,klev)         :: ycoefh_x, ycoefm_x, ycoefh_w, ycoefm_w
    615645    REAL, DIMENSION(klon,klev)         :: ycoefq_x, ycoefq_w
    616     REAL, DIMENSION(klon)              :: ycdragh_x, ycdragm_x, ycdragh_w, ycdragm_w
     646    REAL, DIMENSION(klon)              :: ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w
     647    REAL, DIMENSION(klon)              :: ycdragm_x, ycdragm_w
    617648    REAL, DIMENSION(klon)              :: AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x
    618649    REAL, DIMENSION(klon)              :: AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w
     
    634665    REAL, DIMENSION(klon, klev)        :: zxfluxu_x, zxfluxv_x, zxfluxu_w, zxfluxv_w
    635666    REAL                               :: zx_qs_surf, zcor_surf, zdelta_surf
    636     REAL, DIMENSION(klon)              :: ytsurf_th, yqsatsurf
     667!jyg<
    637668    REAL, DIMENSION(klon)              :: ybeta
     669    REAL, DIMENSION(klon)              :: ybeta_prev
     670!>jyg
    638671    REAL, DIMENSION(klon, klev)        :: d_u_x
    639672    REAL, DIMENSION(klon, klev)        :: d_u_w
     
    770803!!! nrlmd le 13/06/2011
    771804    REAL, DIMENSION(klon)              :: y_delta_flux_t1, y_delta_flux_q1, y_delta_flux_u1, y_delta_flux_v1
    772     REAL, DIMENSION(klon)              :: y_delta_tsurf,delta_coef,tau_eq
     805    REAL, DIMENSION(klon)              :: y_delta_tsurf, y_delta_tsurf_new
     806    REAL, DIMENSION(klon)              :: delta_coef, tau_eq
     807    REAL, DIMENSION(klon)              :: HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn
     808    REAL, DIMENSION(klon)              :: phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0
     809    REAL, DIMENSION(klon)              :: y_delta_qsurf
     810    REAL, DIMENSION(klon)              :: y_delta_qsats
     811    REAL, DIMENSION(klon)              :: yg_T, yg_Q
     812    REAL, DIMENSION(klon)              :: yGamma_dTs_phiT, yGamma_dQs_phiQ
     813    REAL, DIMENSION(klon)              :: ydTs_ins, ydqs_ins
     814!
    773815    REAL, PARAMETER                    :: facteur=2./sqrt(3.14)
    774816    REAL, PARAMETER                    :: inertia=2000.
    775     REAL, DIMENSION(klon)              :: ytsurf_th_x,ytsurf_th_w,yqsatsurf_x,yqsatsurf_w
    776817    REAL, DIMENSION(klon)              :: ydtsurf_th
    777818    REAL                               :: zdelta_surf_x,zdelta_surf_w,zx_qs_surf_x,zx_qs_surf_w
     
    783824    REAL, DIMENSION(klon)              :: Kech_m
    784825    REAL, DIMENSION(klon)              :: Kech_m_x, Kech_m_w
    785     REAL, DIMENSION(klon)              :: yts_x,yts_w
     826    REAL, DIMENSION(klon)              :: yts_x, yts_w
     827    REAL, DIMENSION(klon)              :: yqsatsrf0_x, yqsatsrf0_w
     828    REAL, DIMENSION(klon)              :: yqsurf_x, yqsurf_w
    786829!jyg<
    787830!!    REAL, DIMENSION(klon)              :: Kech_Hp, Kech_H_xp, Kech_H_wp
     
    790833!!    REAL, DIMENSION(klon)              :: Kech_Vp, Kech_V_xp, Kech_V_wp
    791834!>jyg
    792 !jyg<
    793     REAL, DIMENSION(klon)              :: ah, bh     ! coefficients of the delta_Tsurf equation
    794 !>jyg
     835
     836    REAL                               :: fact_cdrag
     837    REAL                               :: z1lay
    795838
    796839    REAL                               :: vent
     
    826869    REAL, DIMENSION(klon)              :: ytoice
    827870    REAL, DIMENSION(klon)              :: ysnowhgt, yqsnow, ysissnow, yrunoff
     871    REAL, DIMENSION(klon)              :: yzmea
    828872    REAL, DIMENSION(klon)              :: yzsig
    829     REAL, DIMENSION(klon,klev)         :: ypphi
    830873    REAL, DIMENSION(klon)              :: ycldt
    831874    REAL, DIMENSION(klon)              :: yrmu0
    832875    ! Martin
    833876
    834     real, DIMENSION(klon):: ydelta_sst, ydelta_sal, yds_ns, ydt_ns, ydter, ydser, &
     877    REAL, DIMENSION(klon):: ydelta_sst, ydelta_sal, yds_ns, ydt_ns, ydter, ydser, &
    835878         ytkt, ytks, ytaur, ysss
    836879    ! compression of delta_sst, delta_sal, ds_ns, dt_ns, dter, dser, tkt, tks,
     
    844887!
    845888!!jyg      iflag_split = mod(iflag_pbl_split,2)
    846       iflag_split = mod(iflag_pbl_split,10)
     889!!jyg      iflag_split = mod(iflag_pbl_split,10)
     890!
     891! Flags controlling the splitting of the turbulent boundary layer:
     892!   iflag_split_ref = 0  ==> no splitting
     893!                   = 1  ==> splitting without coupling with surface temperature
     894!                   = 2  ==> splitting with coupling with surface temperature over land
     895!                   = 3  ==> splitting over ocean; no splitting over land
     896!   iflag_split: actual flag controlling the splitting.
     897!   iflag_split = iflag_split_ref outside the sub-surface loop
     898!               = iflag_split_ref if iflag_split_ref = 0, 1, or 2
     899!               = 0 over land  if iflga_split_ref = 3
     900!               = 1 over ocean if iflga_split_ref = 3
     901
     902      iflag_split_ref = mod(iflag_pbl_split,10)
     903      iflag_split = iflag_split_ref
    847904
    848905!****************************************************************************************
     
    853910
    854911    IF (first_call) THEN
     912
     913       iflag_new_t2mq2m=1
     914       CALL getin_p('iflag_new_t2mq2m',iflag_new_t2mq2m)
     915       WRITE(lunout,*) 'pbl_iflag_new_t2mq2m=',iflag_new_t2mq2m
     916
    855917       print*,'PBL SURFACE AVEC GUSTINESS'
    856918       first_call=.FALSE.
    857919     
    858920       ! Initialize ok_flux_surf (for 1D model)
    859        if (klon_glo>1) ok_flux_surf=.FALSE.
    860        if (klon_glo>1) ok_forc_tsurf=.FALSE.
     921       IF (klon_glo>1) ok_flux_surf=.FALSE.
     922       IF (klon_glo>1) ok_forc_tsurf=.FALSE.
    861923
    862924       ! intialize beta_land
     
    919981 zxfluxlat(:)=0.
    920982 zt2m(:)=0. ; zq2m(:)=0. ; qsat2m(:)=0. ; rh2m(:)=0.
     983 zn2mout(:,:)=0 ;
    921984 d_t(:,:)=0. ; d_t_diss(:,:)=0. ; d_q(:,:)=0. ; d_u(:,:)=0. ; d_v(:,:)=0.
    922985 zcoefh(:,:,:)=0. ; zcoefm(:,:,:)=0.
     
    934997 fder_print(:)=0.
    935998 zxqsurf(:)=0.
     999 delta_qsurf(:) = 0.
    9361000 zxfluxu(:,:)=0. ; zxfluxv(:,:)=0.
    9371001 solsw(:,:)=0. ; sollw(:,:)=0.
     
    10001064    ysnowhgt = 0.0; yqsnow = 0.0     ; yrunoff = 0.0   ; ytoice =0.0
    10011065    yalb3_new = 0.0  ; ysissnow = 0.0
    1002     ypphi = 0.0   ; ycldt = 0.0      ; yrmu0 = 0.0
     1066    ycldt = 0.0      ; yrmu0 = 0.0
    10031067    ! Martin
    10041068
     
    10161080    y_delta_flux_t1=0.
    10171081    ydtsurf_th=0.
    1018     yts_x=0.      ; yts_w=0.
    1019     y_delta_tsurf=0.
     1082    yts_x(:)=0.      ; yts_w(:)=0.
     1083    y_delta_tsurf(:)=0. ; y_delta_qsurf(:)=0.
     1084    yqsurf_x(:)=0.      ; yqsurf_w(:)=0.
     1085    yg_T(:) = 0. ;        yg_Q(:) = 0.
     1086    yGamma_dTs_phiT(:) = 0. ; yGamma_dQs_phiQ(:) = 0.
     1087    ydTs_ins(:) = 0. ; ydqs_ins(:) = 0.
     1088
    10201089!!!
    10211090    ytsoil = 999999.
     
    11921261       DO i = 1, klon
    11931262          sollw(i,nsrf) = sollw_m(i) + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ts(i,nsrf))
    1194 
    1195 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1196 !         ! Martin
    1197 ! Apparently introduced for sisvat but not used
    1198 !         sollwd(i,nsrf)= sollwd_m(i)
    1199 !         ! Martin
    1200 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1201 
    12021263!--OB this line is not satisfactory because alb is the direct albedo not total albedo
    12031264          solsw(i,nsrf) = solsw_m(i) * (1.-alb(i,nsrf)) / (1.-alb_m(i))
     
    12481309!
    12491310!****************************************************************************************
    1250    
    1251     loop_nbsrf: DO nsrf = 1, nbsrf
     1311                                                                          !<<<<<<<<<<<<<
     1312    loop_nbsrf: DO nsrf = 1, nbsrf                                        !<<<<<<<<<<<<<
     1313                                                                          !<<<<<<<<<<<<<
    12521314       IF (prt_level >=10) print *,' Loop nsrf ',nsrf
     1315!
     1316       IF (iflag_split_ref == 3) THEN
     1317         IF (nsrf == is_oce) THEN
     1318            iflag_split = 1
     1319         ELSE
     1320            iflag_split=0
     1321         ENDIF   !! (nsrf == is_oce)
     1322       ELSE                     
     1323         iflag_split = iflag_split_ref
     1324       ENDIF   !! (iflag_split_ref == 3)
    12531325
    12541326! Search for index(ni) and size(knon) of domaine to treat
     
    12861358!****************************************************************************************
    12871359
     1360!
     1361!jyg<    (20190926)
     1362!   Provisional : set ybeta to standard values
     1363       IF (nsrf .NE. is_ter) THEN
     1364           ybeta(:) = 1.
     1365       ELSE
     1366           IF (iflag_split .EQ. 0) THEN
     1367              ybeta(:) = 1.
     1368           ELSE
     1369             DO j = 1, knon
     1370                i = ni(j)
     1371                ybeta(j)   = beta(i,nsrf)
     1372             ENDDO
     1373           ENDIF  ! (iflag_split .LE.1)
     1374       ENDIF !  (nsrf .NE. is_ter)
     1375!>jyg
     1376!
    12881377       DO j = 1, knon
    12891378          i = ni(j)
     
    13181407          ywindsp(j) = windsp(i,nsrf)
    13191408!>jyg
    1320           ! Martin
     1409          ! Martin and Etienne
     1410          yzmea(j)   = zmea(i)
    13211411          yzsig(j)   = zsig(i)
    13221412          ycldt(j)   = cldt(i)
     
    14531543!
    14541544!****************************************************************************************
     1545
    14551546
    14561547!!! jyg le 07/02/2012
     
    15031594           speed_x(i) = SQRT(yu_x(i,1)**2+yv_x(i,1)**2)
    15041595        ENDDO
    1505         CALL cdrag(knon, nsrf, &
     1596
     1597
     1598            CALL cdrag(knon, nsrf, &
    15061599            speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),&
    1507             yts_x, yqsurf, yz0m, yz0h, &
     1600            yts_x, yqsurf_x, yz0m, yz0h, &
    15081601            ycdragm_x, ycdragh_x, zri1_x, pref_x )
    15091602
     
    15321625        CALL cdrag(knon, nsrf, &
    15331626            speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),&
    1534             yts_w, yqsurf, yz0m, yz0h, &
     1627            yts_w, yqsurf_w, yz0m, yz0h, &
    15351628            ycdragm_w, ycdragh_w, zri1_w, pref_w )
    15361629!
     
    16051698       ENDIF
    16061699        CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
    1607             ypaprs, ypplay, yu_x, yv_x, yq_x, yt_x, yts_x, yqsurf, ycdragm_x, &
     1700            ypaprs, ypplay, yu_x, yv_x, yq_x, yt_x, yts_x, yqsurf_x, ycdragm_x, &
    16081701            ycoefm_x, ycoefh_x, ytke_x,y_treedrg)
    16091702!            ycoefm_x, ycoefh_x, ytke_x)
     
    16331726       ENDIF
    16341727        CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
    1635             ypaprs, ypplay, yu_w, yv_w, yq_w, yt_w, yts_w, yqsurf, ycdragm_w, &
     1728            ypaprs, ypplay, yu_w, yv_w, yq_w, yt_w, yts_w, yqsurf_w, ycdragm_w, &
    16361729            ycoefm_w, ycoefh_w, ytke_w,y_treedrg)
    16371730!            ycoefm_w, ycoefh_w, ytke_w)
     
    17701863         yt1(:) = yt(:,1)
    17711864         yq1(:) = yq(:,1)
    1772 !!       ELSE IF (iflag_split .eq. 1) THEN
    1773 !!!
    1774 !jyg<
    1775 !!         CALL wx_pbl_fuse_no_dts(knon, dtime, ypplay, ywake_s, &
    1776 !!                                 yt_x, yt_w, yq_x, yq_w, &
    1777 !!                                 yu_x, yu_w, yv_x, yv_w, &
    1778 !!                                 ycdragh_x, ycdragh_w, ycdragm_x, ycdragm_w, &
    1779 !!                                 AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
    1780 !!                                 AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
    1781 !!                                 BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
    1782 !!                                 BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
    1783 !!                                 AcoefH, AcoefQ, AcoefU, AcoefV, &
    1784 !!                                 BcoefH, BcoefQ, BcoefU, BcoefV, &
    1785 !!                                 ycdragh, ycdragm, &
    1786 !!                                 yt1, yq1, yu1, yv1 &
    1787 !!                                 )
    17881865       ELSE IF (iflag_split .ge. 1) THEN
    1789          CALL wx_pbl0_fuse(knon, dtime, ypplay, ywake_s, &
     1866!
     1867! Cdragq computation
     1868! ------------------
     1869    !******************************************************************************
     1870    ! Cdragq computed from cdrag
     1871    ! The difference comes only from a factor (f_z0qh_oce) on z0, so that
     1872    ! it can be computed inside wx_pbl0_merge
     1873    ! More complicated appraches may require the propagation through
     1874    ! pbl_surface of an independant cdragq variable.
     1875    !******************************************************************************
     1876!
     1877    IF ( f_z0qh_oce .ne. 1. .and. nsrf .eq.is_oce) THEN
     1878       ! Si on suit les formulations par exemple de Tessel, on
     1879       ! a z0h=0.4*nu/u*, z0q=0.62*nu/u*, d'ou f_z0qh_oce=0.62/0.4=1.55
     1880!!       ycdragq_x(1:knon)=ycdragh_x(1:knon)*                                      &
     1881!!            log(z1lay(1:knon)/yz0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*yz0h(1:knon)))
     1882!!       ycdragq_w(1:knon)=ycdragh_w(1:knon)*                                      &
     1883!!            log(z1lay(1:knon)/yz0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*yz0h(1:knon)))
     1884!
     1885       DO j = 1,knon
     1886         z1lay = zgeo1(j)/RG
     1887         fact_cdrag = log(z1lay/yz0h(j))/log(z1lay/(f_z0qh_oce*yz0h(j)))
     1888         ycdragq_x(j)=ycdragh_x(j)*fact_cdrag
     1889         ycdragq_w(j)=ycdragh_w(j)*fact_cdrag
     1890!!     Print *,'YYYYpbl0: fact_cdrag ', fact_cdrag
     1891       ENDDO  ! j = 1,knon
     1892!
     1893!!  Print *,'YYYYpbl0: z1lay, yz0h, f_z0qh_oce, ycdragh_w, ycdragq_w ', &
     1894!!                z1lay, yz0h(1:knon), f_z0qh_oce, ycdragh_w(1:knon), ycdragq_w(1:knon)
     1895    ELSE
     1896       ycdragq_x(1:knon)=ycdragh_x(1:knon)
     1897       ycdragq_w(1:knon)=ycdragh_w(1:knon)
     1898    ENDIF  ! ( f_z0qh_oce .ne. 1. .and. nsrf .eq.is_oce)
     1899!
     1900         CALL wx_pbl_prelim_0(knon, nsrf, dtime, ypplay, ypaprs, ywake_s,  &
     1901                         yts, y_delta_tsurf, ygustiness, &
    17901902                         yt_x, yt_w, yq_x, yq_w, &
    17911903                         yu_x, yu_w, yv_x, yv_w, &
    1792                          ycdragh_x, ycdragh_w, ycdragm_x, ycdragm_w, &
     1904                         ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
     1905                         ycdragm_x, ycdragm_w, &
    17931906                         AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
    17941907                         AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
    17951908                         BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
    17961909                         BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
    1797                          AcoefH, AcoefQ, AcoefU, AcoefV, &
    1798                          BcoefH, BcoefQ, BcoefU, BcoefV, &
    1799                          ycdragh, ycdragm, &
     1910                         Kech_h_x, Kech_h_w, Kech_h  &
     1911                         )
     1912         CALL wx_pbl_prelim_beta(knon, dtime, ywake_s, ybeta,  &
     1913                         BcoefQ_x, BcoefQ_w  &
     1914                         )
     1915         CALL wx_pbl0_merge(knon, ypplay, ypaprs,  &
     1916                         ywake_s, ydTs0, ydqs0, &
     1917                         yt_x, yt_w, yq_x, yq_w, &
     1918                         yu_x, yu_w, yv_x, yv_w, &
     1919                         ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
     1920                         ycdragm_x, ycdragm_w, &
     1921                         AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
     1922                         AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
     1923                         BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
     1924                         BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
     1925                         AcoefH_0, AcoefQ_0, AcoefU, AcoefV, &
     1926                         BcoefH_0, BcoefQ_0, BcoefU, BcoefV, &
     1927                         ycdragh, ycdragq, ycdragm, &
    18001928                         yt1, yq1, yu1, yv1 &
    18011929                         )
    1802 !!       ELSE IF (iflag_split .ge.2) THEN
    1803 !!!    Provisoire
    1804 !!         ah(:) = 0.
    1805 !!         bh(:) = 0.
    1806 !!         IF (nsrf == is_oce) THEN
    1807 !!           ybeta(:) = 1.
    1808 !!         ELSE
    1809 !!           ybeta(:) = beta_land
    1810 !!         ENDIF
    1811 !!         ycdragh(:) = ywake_s(:)*ycdragh_w(:) + (1.-ywake_s(:))*ycdragh_x(:)
    1812 !!         CALL wx_dts(knon, nsrf, ywake_cstar, ywake_s, ywake_dens, &
    1813 !!                     yts, ypplay(:,1), ybeta, ycdragh , ypaprs(:,1), &
    1814 !!                     yq(:,1), yt(:,1), yu(:,1), yv(:,1), ygustiness, &
    1815 !!                     ah, bh &
    1816 !!                     )
    1817 !!!
    1818 !!         CALL wx_pbl_fuse(knon, dtime, ypplay, ywake_s, &
    1819 !!                         yt_x, yt_w, yq_x, yq_w, &
    1820 !!                         yu_x, yu_w, yv_x, yv_w, &
    1821 !!                         ycdragh_x, ycdragh_w, ycdragm_x, ycdragm_w, &
    1822 !!                         AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
    1823 !!                         AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
    1824 !!                         BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
    1825 !!                         BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
    1826 !!                         ah, bh, &
    1827 !!                         AcoefH, AcoefQ, AcoefU, AcoefV, &
    1828 !!                         BcoefH, BcoefQ, BcoefU, BcoefV, &
    1829 !!                         ycdragh, ycdragm, &
    1830 !!                         yt1, yq1, yu1, yv1 &
    1831 !!                         )
    1832 !>jyg
    1833 !!!
    1834          ENDIF  ! (iflag_split .eq.0)
     1930         IF (iflag_split .eq. 2 .AND. nsrf .ne. is_oce) THEN
     1931           CALL wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, &
     1932                           ywake_s, ybeta, ywake_cstar, ywake_dens, &
     1933                           AcoefH_x, AcoefH_w, &
     1934                           BcoefH_x, BcoefH_w, &
     1935                           AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0,  &
     1936                           AcoefH, AcoefQ, BcoefH, BcoefQ,  &
     1937                           HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
     1938                           phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
     1939                           yg_T, yg_Q, &
     1940                           yGamma_dTs_phiT, yGamma_dQs_phiQ, &
     1941                           ydTs_ins, ydqs_ins &
     1942                           )
     1943         ELSE !
     1944           AcoefH(:) = AcoefH_0(:)
     1945           AcoefQ(:) = AcoefQ_0(:)
     1946           BcoefH(:) = BcoefH_0(:)
     1947           BcoefQ(:) = BcoefQ_0(:)
     1948           yg_T(:) = 0.
     1949           yg_Q(:) = 0.
     1950           yGamma_dTs_phiT(:) = 0.
     1951           yGamma_dQs_phiQ(:) = 0.
     1952           ydTs_ins(:) = 0.
     1953           ydqs_ins(:) = 0.
     1954         ENDIF   ! (iflag_split .eq. 2)
     1955       ENDIF  ! (iflag_split .eq.0)
    18351956!!!
    18361957       IF (prt_level >=10) THEN
    1837          PRINT *,'pbl_surface (fuse->): yt(1,:) ',yt(1,:)
    1838          PRINT *,'pbl_surface (fuse->): yq(1,:) ',yq(1,:)
    1839          PRINT *,'pbl_surface (fuse->): yu(1,:) ',yu(1,:)
    1840          PRINT *,'pbl_surface (fuse->): yv(1,:) ',yv(1,:)
    1841          PRINT *,'pbl_surface (fuse->): AcoefH(1) ',AcoefH(1)
    1842          PRINT *,'pbl_surface (fuse->): BcoefH(1) ',BcoefH(1)
     1958         PRINT *,'pbl_surface (merge->): yt(1,:) ',yt(1,:)
     1959         PRINT *,'pbl_surface (merge->): yq(1,:) ',yq(1,:)
     1960         PRINT *,'pbl_surface (merge->): yu(1,:) ',yu(1,:)
     1961         PRINT *,'pbl_surface (merge->): yv(1,:) ',yv(1,:)
     1962         PRINT *,'pbl_surface (merge->): AcoefH(1), AcoefQ(1), AcoefU(1), AcoefV(1) ', &
     1963                                         AcoefH(1), AcoefQ(1), AcoefU(1), AcoefV(1)
     1964         PRINT *,'pbl_surface (merge->): BcoefH(1), BcoefQ(1), BcoefU(1), BcoefV(1) ', &
     1965                                         BcoefH(1), BcoefQ(1), BcoefU(1), BcoefV(1)
     1966
    18431967       ENDIF
    18441968
     1969!  Save initial value of z0h for use in evappot (z0h wiil be computed again in the surface models)
     1970          yz0h_old(1:knon) = yz0h(1:knon)
     1971!
    18451972!****************************************************************************************
    18461973!
     
    18571984
    18581985          ! Calculate the temperature et relative humidity at 2m and the wind at 10m
     1986          IF (iflag_new_t2mq2m==1) THEN
     1987           CALL stdlevvarn(klon, knon, is_ter, zxli, &
     1988               yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
     1989               yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
     1990               yt2m, yq2m, yt10m, yq10m, yu10m, yustar, &
     1991               yn2mout(:, nsrf, :))
     1992          ELSE
    18591993          CALL stdlevvar(klon, knon, is_ter, zxli, &
    18601994               yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
    18611995               yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
    18621996               yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
     1997          ENDIF
    18631998         
    18641999       ENDIF
     
    19232058          CALL surf_landice(itap, dtime, knon, ni, &
    19242059               rlon, rlat, debut, lafin, &
    1925                yrmu0, ylwdown, yalb, ypphi(:,1), &
     2060               yrmu0, ylwdown, yalb, zgeo1, &
    19262061               ysolsw, ysollw, yts, ypplay(:,1), &
    19272062!!jyg               ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
     
    19332068               ytsoil, yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap,yfluxsens,yfluxlat, &
    19342069               ytsurf_new, y_dflux_t, y_dflux_q, &
    1935                yzsig, ycldt, &
     2070               yzmea, yzsig, ycldt, &
    19362071               ysnowhgt, yqsnow, ytoice, ysissnow, &
    19372072               yalb3_new, yrunoff, &
     
    20932228          y_flux_q1(j) = -yevap(j)
    20942229          ENDDO
    2095         ENDIF
    2096 
    2097        IF (prt_level >=10) THEN
    2098         DO j=1,knon
    2099          print*,'y_flux_t1,yfluxlat,wakes' &
    2100  &             ,  y_flux_t1(j), yfluxlat(j), ywake_s(j)
    2101          print*,'beta,ytsurf_new', ybeta(j), ytsurf_new(j)
    2102          print*,'inertia,facteur,cstar', inertia, facteur,wake_cstar(j)
    2103         ENDDO
    2104        ENDIF
    2105 
    2106 !!! jyg le 07/02/2012 puis le 10/04/2013
    2107 !!       IF (iflag_split .eq.1) THEN
    2108 !!!!!
    2109 !!!jyg<
    2110 !!         CALL wx_pbl_split_no_dts(knon, ywake_s, &
    2111 !!                                AcoefH_x, AcoefH_w, &
    2112 !!                                AcoefQ_x, AcoefQ_w, &
    2113 !!                                AcoefU_x, AcoefU_w, &
    2114 !!                                AcoefV_x, AcoefV_w, &
    2115 !!                                y_flux_t1, y_flux_q1, y_flux_u1, y_flux_v1, &
    2116 !!                                y_flux_t1_x, y_flux_t1_w, &
    2117 !!                                y_flux_q1_x, y_flux_q1_w, &
    2118 !!                                y_flux_u1_x, y_flux_u1_w, &
    2119 !!                                y_flux_v1_x, y_flux_v1_w, &
    2120 !!                                yfluxlat_x, yfluxlat_w &
    2121 !!                                )
    2122 !!       ELSE IF (iflag_split .ge. 2) THEN
     2230        ENDIF ! (ok_flux_surf)
     2231!
     2232! ------------------------------------------------------------------------------
     2233! 12a)  Splitting
     2234! ------------------------------------------------------------------------------
     2235
    21232236       IF (iflag_split .GE. 1) THEN
    2124          CALL wx_pbl0_split(knon, dtime, ywake_s, &
     2237!
     2238         IF (nsrf .ne. is_oce) THEN
     2239!
     2240!         Compute potential evaporation and aridity factor  (jyg, 20200328)
     2241          ybeta_prev(:) = ybeta(:)
     2242             DO j = 1, knon
     2243               yqa(j) = AcoefQ(j) - BcoefQ(j)*yevap(j)*dtime
     2244             ENDDO
     2245!
     2246          CALL wx_evappot(knon, yqa, yTsurf_new, yevap_pot)
     2247!
     2248          ybeta(1:knon) = min(yevap(1:knon)/yevap_pot(1:knon), 1.)
     2249         
     2250          IF (prt_level >=10) THEN
     2251           DO j=1,knon
     2252            print*,'y_flux_t1,yfluxlat,wakes' &
     2253 &                ,  y_flux_t1(j), yfluxlat(j), ywake_s(j)
     2254            print*,'beta_prev, beta, ytsurf_new', ybeta_prev(j), ybeta(j), ytsurf_new(j)
     2255            print*,'inertia,facteur,cstar', inertia, facteur,wake_cstar(j)
     2256           ENDDO
     2257          ENDIF  ! (prt_level >=10)
     2258!
     2259! Second call to wx_pbl0_merge and wx_pbl_dts_merge in order to take into account
     2260! the update of the aridity coeficient beta.
     2261!
     2262        CALL wx_pbl_prelim_beta(knon, dtime, ywake_s, ybeta,  &
     2263                        BcoefQ_x, BcoefQ_w  &
     2264                        )
     2265        CALL wx_pbl0_merge(knon, ypplay, ypaprs,  &
     2266                          ywake_s, ydTs0, ydqs0, &
     2267                          yt_x, yt_w, yq_x, yq_w, &
     2268                          yu_x, yu_w, yv_x, yv_w, &
     2269                          ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
     2270                          ycdragm_x, ycdragm_w, &
     2271                          AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
     2272                          AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
     2273                          BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
     2274                          BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
     2275                          AcoefH_0, AcoefQ_0, AcoefU, AcoefV, &
     2276                          BcoefH_0, BcoefQ_0, BcoefU, BcoefV, &
     2277                          ycdragh, ycdragq, ycdragm, &
     2278                          yt1, yq1, yu1, yv1 &
     2279                          )
     2280          IF (iflag_split .eq. 2) THEN
     2281            CALL wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, &
     2282                            ywake_s, ybeta, ywake_cstar, ywake_dens, &
     2283                            AcoefH_x, AcoefH_w, &
     2284                            BcoefH_x, BcoefH_w, &
     2285                            AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0,  &
     2286                            AcoefH, AcoefQ, BcoefH, BcoefQ,  &
     2287                            HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
     2288                            phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
     2289                            yg_T, yg_Q, &
     2290                            yGamma_dTs_phiT, yGamma_dQs_phiQ, &
     2291                            ydTs_ins, ydqs_ins &
     2292                            )
     2293          ELSE !
     2294            AcoefH(:) = AcoefH_0(:)
     2295            AcoefQ(:) = AcoefQ_0(:)
     2296            BcoefH(:) = BcoefH_0(:)
     2297            BcoefQ(:) = BcoefQ_0(:)
     2298            yg_T(:) = 0.
     2299            yg_Q(:) = 0.
     2300            yGamma_dTs_phiT(:) = 0.
     2301            yGamma_dQs_phiQ(:) = 0.
     2302            ydTs_ins(:) = 0.
     2303            ydqs_ins(:) = 0.
     2304          ENDIF   ! (iflag_split .eq. 2)
     2305!
     2306        ELSE    ! (nsrf .ne. is_oce)
     2307          ybeta(1:knon) = 1.
     2308          yevap_pot(1:knon) = yevap(1:knon)
     2309          AcoefH(:) = AcoefH_0(:)
     2310          AcoefQ(:) = AcoefQ_0(:)
     2311          BcoefH(:) = BcoefH_0(:)
     2312          BcoefQ(:) = BcoefQ_0(:)
     2313          yg_T(:) = 0.
     2314          yg_Q(:) = 0.
     2315          yGamma_dTs_phiT(:) = 0.
     2316          yGamma_dQs_phiQ(:) = 0.
     2317          ydTs_ins(:) = 0.
     2318          ydqs_ins(:) = 0.
     2319        ENDIF   ! (nsrf .ne. is_oce)
     2320!
     2321        CALL wx_pbl_split(knon, nsrf, dtime, ywake_s, ybeta, iflag_split, &
     2322                       yg_T, yg_Q, &
     2323                       yGamma_dTs_phiT, yGamma_dQs_phiQ, &
     2324                       ydTs_ins, ydqs_ins, &
    21252325                       y_flux_t1, y_flux_q1, y_flux_u1, y_flux_v1, &
     2326!!!!                       HTRn_b, dd_HTRn, HTphiT_b, dd_HTphiT, &
     2327                       phiQ0_b, phiT0_b, &
    21262328                       y_flux_t1_x, y_flux_t1_w, &
    21272329                       y_flux_q1_x, y_flux_q1_w, &
     
    21292331                       y_flux_v1_x, y_flux_v1_w, &
    21302332                       yfluxlat_x, yfluxlat_w, &
    2131                        y_delta_tsurf &
     2333                       y_delta_qsats, &
     2334                       y_delta_tsurf_new, y_delta_qsurf &
    21322335                       )
     2336!
     2337         CALL wx_pbl_check(knon, dtime, ypplay, ypaprs, ywake_s, ybeta, iflag_split, &
     2338                       yTs, y_delta_tsurf,  &
     2339                       yqsurf, yTsurf_new,  &
     2340                       y_delta_tsurf_new, y_delta_qsats,  &
     2341                       AcoefH_x, AcoefH_w, &
     2342                       BcoefH_x, BcoefH_w, &
     2343                       AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0,  &
     2344                       AcoefH, AcoefQ, BcoefH, BcoefQ,  &
     2345                       y_flux_t1, y_flux_q1,  &
     2346                       y_flux_t1_x, y_flux_t1_w, &
     2347                       y_flux_q1_x, y_flux_q1_w)
     2348!
     2349         IF (nsrf .ne. is_oce) THEN
     2350           CALL wx_pbl_dts_check(knon, dtime, ypplay, ypaprs, ywake_s, ybeta, iflag_split, &
     2351                         yTs, y_delta_tsurf,  &
     2352                         yqsurf, yTsurf_new,  &
     2353                         y_delta_qsats, y_delta_tsurf_new, y_delta_qsurf,  &
     2354                         AcoefH_x, AcoefH_w, &
     2355                         BcoefH_x, BcoefH_w, &
     2356                         AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0,  &
     2357                         AcoefH, AcoefQ, BcoefH, BcoefQ,  &
     2358                         HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
     2359                         phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
     2360                         yg_T, yg_Q, &
     2361                         yGamma_dTs_phiT, yGamma_dQs_phiQ, &
     2362                         ydTs_ins, ydqs_ins, &
     2363                         y_flux_t1, y_flux_q1,  &
     2364                         y_flux_t1_x, y_flux_t1_w, &
     2365                         y_flux_q1_x, y_flux_q1_w )
     2366         ENDIF   ! (nsrf .ne. is_oce)
     2367!
     2368       ELSE  ! (iflag_split .ge. 1)
     2369         ybeta(1:knon) = 1.
     2370         yevap_pot(1:knon) = yevap(1:knon)
    21332371       ENDIF  ! (iflag_split .ge. 1)
     2372!
     2373       IF (prt_level >= 10) THEN
     2374         print *,'pbl_surface, ybeta , yevap, yevap_pot ', &
     2375                               ybeta , yevap, yevap_pot
     2376       ENDIF  ! (prt_level >= 10)
     2377!
    21342378!>jyg
    21352379!
     
    21802424         print*,'Chx,Chw,Ch', ycdragh_x(j), ycdragh_w(j), ycdragh(j)
    21812425         print*,'Khx,Khw,Kh', Kech_h_x(j), Kech_h_w(j), Kech_h(j)
    2182 !         print*,'tsurf_x,tsurf_w,tsurf,t1', ytsurf_th_x(j), ytsurf_th_w(j), ytsurf_th(j), yt(j,1)
    2183          print*,'tsurf_x,t1x,tsurf_w,t1w,tsurf,t1,t1_ancien', &
    2184  &               ytsurf_th_x(j), yt_x(j,1), ytsurf_th_w(j), yt_w(j,1), ytsurf_th(j), yt(j,1),t(j,1)
    2185          print*,'qsatsurf,qsatsurf_x,qsatsurf_w', yqsatsurf(j), yqsatsurf_x(j), yqsatsurf_w(j)
     2426         print*,'t1x, t1w, t1, t1_ancien', &
     2427 &               yt_x(j,1), yt_w(j,1),  yt(j,1), t(j,1)
    21862428         print*,'delta_coef,delta_flux,delta_tsurf,tau', delta_coef(j), y_delta_flux_t1(j), y_delta_tsurf(j), tau_eq(j)
    21872429        ENDDO
     
    21902432         print*,'fluxT_x, fluxT_w, y_flux_t1, fluxQ_x, fluxQ_w, yfluxlat, wakes' &
    21912433 &             , y_flux_t1_x(j), y_flux_t1_w(j), y_flux_t1(j), y_flux_q1_x(j)*RLVTT, y_flux_q1_w(j)*RLVTT, yfluxlat(j), ywake_s(j)
    2192          print*,'beta,ytsurf_new,yqsatsurf', ybeta(j), ytsurf_new(j), yqsatsurf(j)
    2193          print*,'inertia,facteur,cstar', inertia, facteur,wake_cstar(j)
     2434         print*,'beta, ytsurf_new ', ybeta(j), ytsurf_new(j)
     2435         print*,'inertia, facteur, cstar', inertia, facteur,wake_cstar(j)
    21942436        ENDDO
    21952437       ENDIF  ! (prt_level >=10)
     
    22942536       ENDIF  ! (iflag_split .eq.0)
    22952537!!!
    2296 
    2297         DO j = 1, knon
    2298           y_dflux_t(j) = y_dflux_t(j) * ypct(j)
    2299           y_dflux_q(j) = y_dflux_q(j) * ypct(j)
    2300         ENDDO
    2301 
     2538!!
     2539!!        DO j = 1, knon
     2540!!          y_dflux_t(j) = y_dflux_t(j) * ypct(j)
     2541!!          y_dflux_q(j) = y_dflux_q(j) * ypct(j)
     2542!!        ENDDO
     2543!!
    23022544!****************************************************************************************
    23032545! 13) Transform variables for output format :
     
    24142656          i = ni(j)
    24152657          evap(i,nsrf) = - flux_q(i,1,nsrf)                  !jyg
     2658          beta(i,nsrf) = ybeta(j)                             !jyg
    24162659          d_ts(i,nsrf) = y_d_ts(j)
    24172660!albedo SB >>>
     
    24292672          cdragh(i) = cdragh(i) + ycdragh(j)*ypct(j)
    24302673          cdragm(i) = cdragm(i) + ycdragm(j)*ypct(j)
    2431           dflux_t(i) = dflux_t(i) + y_dflux_t(j)
    2432           dflux_q(i) = dflux_q(i) + y_dflux_q(j)
     2674          dflux_t(i) = dflux_t(i) + y_dflux_t(j)*ypct(j)
     2675          dflux_q(i) = dflux_q(i) + y_dflux_q(j)*ypct(j)
    24332676       ENDDO
    24342677
     
    24462689!!! nrlmd le 13/06/2011
    24472690!!jyg20170131          delta_tsurf(i,nsrf)=y_delta_tsurf(j)*ypct(j)
    2448           delta_tsurf(i,nsrf)=y_delta_tsurf(j)
     2691!!jyg20210118          delta_tsurf(i,nsrf)=y_delta_tsurf(j)
     2692          delta_tsurf(i,nsrf)=y_delta_tsurf_new(j)
     2693!
     2694          delta_qsurf(i) = delta_qsurf(i) + y_delta_qsurf(j)*ypct(j)
    24492695!
    24502696          cdragh_x(i) = cdragh_x(i) + ycdragh_x(j)*ypct(j)
     
    26102856          sss(ni(:knon)) = ysss(:knon)
    26112857       end if
     2858
    26122859
    26132860!****************************************************************************************
     
    26472894               * (ypaprs(j,1)-ypplay(j,1))
    26482895          tairsol(j) = yts(j) + y_d_ts(j)
    2649           tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf(j)
     2896!!          tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf(j)
     2897          tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf_new(j)
    26502898          qairsol(j) = yqsurf(j)
    26512899        ENDDO
     
    26862934!!! jyg le 07/02/2012
    26872935       IF (iflag_split .eq.0) THEN
     2936        IF (iflag_new_t2mq2m==1) THEN
     2937         CALL stdlevvarn(klon, knon, nsrf, zxli, &
     2938            uzon, vmer, tair1, qair1, zgeo1, &
     2939            tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
     2940            yt2m, yq2m, yt10m, yq10m, yu10m, yustar, &
     2941            yn2mout(:, nsrf, :))
     2942        ELSE
    26882943        CALL stdlevvar(klon, knon, nsrf, zxli, &
    26892944            uzon, vmer, tair1, qair1, zgeo1, &
    26902945            tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
    26912946            yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
     2947        ENDIF
    26922948       ELSE  !(iflag_split .eq.0)
     2949        IF (iflag_new_t2mq2m==1) THEN
     2950         CALL stdlevvarn(klon, knon, nsrf, zxli, &
     2951            uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
     2952            tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, &
     2953            yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x, &
     2954            yn2mout_x(:, nsrf, :))
     2955         CALL stdlevvarn(klon, knon, nsrf, zxli, &
     2956            uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
     2957            tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
     2958            yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w, &
     2959            yn2mout_w(:, nsrf, :))
     2960        ELSE
    26932961        CALL stdlevvar(klon, knon, nsrf, zxli, &
    26942962            uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
     
    26992967            tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
    27002968            yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w)
     2969        ENDIF
    27012970!!!
    27022971       ENDIF  ! (iflag_split .eq.0)
     
    27122981          u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)**2+vmer(j)**2)
    27132982          v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)**2+vmer(j)**2)
     2983!
     2984          DO k = 1, 6
     2985           n2mout(i,nsrf,k) = yn2mout(j,nsrf,k)
     2986          END DO 
     2987!
    27142988        ENDDO
    27152989       ELSE  !(iflag_split .eq.0)
     
    27222996          u10m_x(i,nsrf)=(yu10m_x(j) * uzon_x(j))/SQRT(uzon_x(j)**2+vmer_x(j)**2)
    27232997          v10m_x(i,nsrf)=(yu10m_x(j) * vmer_x(j))/SQRT(uzon_x(j)**2+vmer_x(j)**2)
     2998!
     2999          DO k = 1, 6
     3000           n2mout_x(i,nsrf,k) = yn2mout_x(j,nsrf,k)
     3001          END DO 
     3002!
    27243003        ENDDO
    27253004        DO j=1, knon
     
    27353014          u10m(i,nsrf) = u10m_x(i,nsrf) + wake_s(i)*(u10m_w(i,nsrf)-u10m_x(i,nsrf))
    27363015          v10m(i,nsrf) = v10m_x(i,nsrf) + wake_s(i)*(v10m_w(i,nsrf)-v10m_x(i,nsrf))
     3016!
     3017          DO k = 1, 6
     3018           n2mout_w(i,nsrf,k) = yn2mout_w(j,nsrf,k)
     3019          END DO 
     3020!
    27373021        ENDDO
    27383022!!!
     
    29173201!****************************************************************************************
    29183202    ENDDO loop_nbsrf
     3203!
     3204!----------------------------------------------------------------------------------------
     3205!   Reset iflag_split
     3206!
     3207   iflag_split=iflag_split_ref
    29193208
    29203209!****************************************************************************************
     
    29863275    ENDDO
    29873276!!!
    2988    
     3277
    29893278!
    29903279! Incrementer la temperature du sol
    29913280!
    29923281    zxtsol(:) = 0.0  ; zxfluxlat(:) = 0.0
    2993     zt2m(:) = 0.0    ; zq2m(:) = 0.0
     3282    zt2m(:) = 0.0    ; zq2m(:) = 0.0 ; zn2mout(:,:) = 0
    29943283    zustar(:)=0.0 ; zu10m(:) = 0.0   ; zv10m(:) = 0.0
    29953284    s_pblh(:) = 0.0  ; s_plcl(:) = 0.0
     
    30443333          zt2m(i)  = zt2m(i)  + t2m(i,nsrf)  * pctsrf(i,nsrf)
    30453334          zq2m(i)  = zq2m(i)  + q2m(i,nsrf)  * pctsrf(i,nsrf)
     3335!
     3336          DO k = 1, 6
     3337           zn2mout(i,k)  = zn2mout(i,k)  + n2mout(i,nsrf,k)  * pctsrf(i,nsrf)
     3338          ENDDO 
     3339!
    30463340          zustar(i) = zustar(i) + ustar(i,nsrf) * pctsrf(i,nsrf)
    30473341          wstar(i,is_ave)=wstar(i,is_ave)+wstar(i,nsrf)*pctsrf(i,nsrf)
     
    30753369          zt2m(i)  = zt2m(i)  + (t2m_x(i,nsrf)+wake_s(i)*(t2m_w(i,nsrf)-t2m_x(i,nsrf))) * pctsrf(i,nsrf)
    30763370          zq2m(i)  = zq2m(i)  + q2m_x(i,nsrf)  * pctsrf(i,nsrf)
     3371!
     3372          DO k = 1, 6
     3373           zn2mout(i,k)  = zn2mout(i,k)  + n2mout_x(i,nsrf,k)  * pctsrf(i,nsrf)
     3374          ENDDO
     3375!
    30773376          zustar(i) = zustar(i) + ustar_x(i,nsrf) * pctsrf(i,nsrf)
    30783377          wstar(i,is_ave)=wstar(i,is_ave)+wstar_x(i,nsrf)*pctsrf(i,nsrf)
     
    31533452    DO nsrf = 1, nbsrf
    31543453       DO i = 1, klon
    3155           zxqsurf(i) = zxqsurf(i) + qsurf(i,nsrf) * pctsrf(i,nsrf)
     3454          zxqsurf(i) = zxqsurf(i) + MAX(qsurf(i,nsrf),0.0) * pctsrf(i,nsrf)
    31563455          zxsnow(i)  = zxsnow(i)  + snow(i,nsrf)  * pctsrf(i,nsrf)
    31573456       ENDDO
     
    31983497    IF (ALLOCATED(qsurf)) DEALLOCATE(qsurf)
    31993498    IF (ALLOCATED(ftsoil)) DEALLOCATE(ftsoil)
     3499    IF (ALLOCATED(ydTs0)) DEALLOCATE(ydTs0)
     3500    IF (ALLOCATED(ydqs0)) DEALLOCATE(ydqs0)
    32003501
    32013502!jyg<
Note: See TracChangeset for help on using the changeset viewer.