Changeset 5869 for LMDZ6/branches/PBLSURF_GPUPORT
- Timestamp:
- Nov 17, 2025, 4:01:45 PM (5 days ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/PBLSURF_GPUPORT/libf/phylmd/pbl_surface_mod.F90
r5868 r5869 324 324 325 325 END SUBROUTINE pbl_surface_init_iso 326 #endif327 328 329 #ifndef toSupress330 331 !332 !****************************************************************************************333 !334 335 SUBROUTINE pbl_surface( &336 dtime, date0, itap, jour, &337 debut, lafin, &338 rlon, rlat, rugoro, rmu0, &339 !GG lwdown_m, cldt, &340 lwdown_m, pphi, cldt, &341 !GG342 rain_f, snow_f, bs_f, solsw_m, solswfdiff_m, sollw_m, &343 gustiness, &344 t, q, qbs, u, v, &345 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012346 !! t_x, q_x, t_w, q_w, &347 wake_dlt, wake_dlq, &348 wake_cstar, wake_s, &349 !!!350 pplay, paprs, pctsrf, &351 ts,SFRWL, alb_dir, alb_dif,ustar, u10m, v10m,wstar, &352 cdragh, cdragm, zu1, zv1, &353 !jyg< (26/09/2019)354 beta, &355 !>jyg356 alb_dir_m, alb_dif_m, zxsens, zxevap, zxsnowerosion, &357 icesub_lic, alb3_lic, runoff, snowhgt, qsnow, to_ice, sissnow, &358 zxtsol, zxfluxlat, zt2m, qsat2m, zn2mout, &359 d_t, d_q, d_qbs, d_u, d_v, d_t_diss, &360 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012361 d_t_w, d_q_w, &362 d_t_x, d_q_x, &363 !! d_wake_dlt,d_wake_dlq, &364 zxsens_x, zxfluxlat_x,zxsens_w,zxfluxlat_w, &365 !!!366 !!! nrlmd le 13/06/2011367 delta_tsurf,wake_dens,cdragh_x,cdragh_w, &368 cdragm_x,cdragm_w,kh,kh_x,kh_w, &369 !!!370 zcoefh, zcoefm, slab_wfbils, &371 qsol, zq2m, s_pblh, s_plcl, &372 !!!373 !!! jyg le 08/02/2012374 s_pblh_x, s_plcl_x, s_pblh_w, s_plcl_w, &375 !!!376 s_capCL, s_oliqCL, s_cteiCL, s_pblT, &377 s_therm, s_trmb1, s_trmb2, s_trmb3, &378 zustar,zu10m, zv10m, fder_print, &379 zxqsurf, delta_qsurf, &380 rh2m, zxfluxu, zxfluxv, &381 z0m, z0h, agesno, sollw, solsw, &382 d_ts, evap, fluxlat, t2m, &383 wfbils, wfevap, &384 flux_t, flux_u, flux_v, &385 dflux_t, dflux_q, zxsnow, &386 !jyg<387 !! zxfluxt, zxfluxq, q2m, flux_q, tke, &388 zxfluxt, zxfluxq, zxfluxqbs, q2m, flux_q, flux_qbs, tke_x, eps_x, &389 !>jyg390 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012391 !! tke_x, tke_w &392 wake_dltke, &393 !GG treedrg &394 treedrg,hice ,tice, bilg_cumul, &395 fcds, fcdi, dh_basal_growth, dh_basal_melt, &396 dh_top_melt, dh_snow2sic, &397 dtice_melt, dtice_snow2sic , &398 !GG399 !FC400 !AM heterogeneous continental sub-surfaces401 tsurf_tersrf, tsoil_tersrf, qsurf_tersrf, tsurf_new_tersrf, &402 cdragm_tersrf, cdragh_tersrf, &403 swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf &404 !!!405 #ifdef ISO406 & ,xtrain_f, xtsnow_f,xt, &407 & wake_dlxt,zxxtevap,xtevap, &408 & d_xt,d_xt_w,d_xt_x, &409 & xtsol,dflux_xt,zxxtsnow,zxfluxxt,flux_xt, &410 & h1_diag,runoff_diag,xtrunoff_diag &411 #endif412 & )413 !****************************************************************************************414 ! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818415 ! Objet: interface de "couche limite" (diffusion verticale)416 !417 !AA REM:418 !AA-----419 !AA Tout ce qui a trait au traceurs est dans phytrac maintenant420 !AA pour l'instant le calcul de la couche limite pour les traceurs421 !AA se fait avec cltrac et ne tient pas compte de la differentiation422 !AA des sous-fraction de sol.423 !AA REM bis :424 !AA----------425 !AA Pour pouvoir extraire les coefficient d'echanges et le vent426 !AA dans la premiere couche, 3 champs supplementaires ont ete crees427 !AA zcoefh, zu1 et zv1. Pour l'instant nous avons moyenne les valeurs428 !AA de ces trois champs sur les 4 subsurfaces du modele. Dans l'avenir429 !AA si les informations des subsurfaces doivent etre prises en compte430 !AA il faudra sortir ces memes champs en leur ajoutant une dimension,431 !AA c'est a dire nbsrf (nbre de subsurface).432 !433 ! Arguments:434 !435 ! dtime----input-R- interval du temps (secondes)436 ! itap-----input-I- numero du pas de temps437 ! date0----input-R- jour initial438 ! t--------input-R- temperature (K)439 ! q--------input-R- vapeur d'eau (kg/kg)440 ! u--------input-R- vitesse u441 ! v--------input-R- vitesse v442 ! wake_dlt-input-R- temperatre difference between (w) and (x) (K)443 ! wake_dlq-input-R- humidity difference between (w) and (x) (kg/kg)444 !wake_cstar-input-R- wake gust front speed (m/s)445 ! wake_s---input-R- wake fractionnal area446 ! ts-------input-R- temperature du sol (en Kelvin)447 ! paprs----input-R- pression a intercouche (Pa)448 ! pplay----input-R- pression au milieu de couche (Pa)449 ! rlat-----input-R- latitude en degree450 ! z0m, z0h ----input-R- longeur de rugosite (en m)451 ! Martin452 ! cldt-----input-R- total cloud fraction453 ! Martin454 !GG455 ! pphi-----input-R- geopotentiel de chaque couche (g z) (reference sol)456 !GG457 !458 ! d_t------output-R- le changement pour "t"459 ! d_q------output-R- le changement pour "q"460 ! d_u------output-R- le changement pour "u"461 ! d_v------output-R- le changement pour "v"462 ! d_ts-----output-R- le changement pour "ts"463 ! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)464 ! (orientation positive vers le bas)465 ! tke_x---input/output-R- tke in the (x) region (kg/m**2/s)466 ! wake_dltke-input/output-R- tke difference between (w) and (x) (kg/m**2/s)467 ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)468 ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal469 ! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal470 ! dflux_t--output-R- derive du flux sensible471 ! dflux_q--output-R- derive du flux latent472 ! zu1------output-R- le vent dans la premiere couche473 ! zv1------output-R- le vent dans la premiere couche474 ! trmb1----output-R- deep_cape475 ! trmb2----output-R- inhibition476 ! trmb3----output-R- Point Omega477 ! cteiCL---output-R- Critere d'instab d'entrainmt des nuages de CL478 ! plcl-----output-R- Niveau de condensation479 ! pblh-----output-R- HCL480 ! pblT-----output-R- T au nveau HCL481 ! treedrg--output-R- tree drag (m)482 ! qsurf_tersrf--output-R- surface specific humidity of continental sub-surfaces483 ! cdragm_tersrf--output-R- momentum drag coefficient of continental sub-surfaces484 ! cdragh_tersrf--output-R- heat drag coefficient of continental sub-surfaces485 ! tsurf_new_tersrf--output-R- surface temperature of continental sub-surfaces486 ! swnet_tersrf--output-R- net shortwave radiation of continental sub-surfaces487 ! lwnet_tersrf--output-R- net longwave radiation of continental sub-surfaces488 ! fluxsens_tersrf--output-R- sensible heat flux of continental sub-surfaces489 ! fluxlat_tersrf--output-R- latent heat flux of continental sub-surfaces490 491 USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, carbon_cycle_tr, level_coupling_esm492 USE carbon_cycle_mod, ONLY : co2_send, nbcf_out, fields_out, yfields_out, cfname_out493 use hbtm_mod, only: hbtm494 USE indice_sol_mod495 USE time_phylmdz_mod, ONLY : day_ini,annee_ref,itau_phy496 USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat, grid1dto2d_glo497 USE print_control_mod, ONLY : prt_level,lunout498 #ifdef ISO499 USE isotopes_mod, ONLY: Rdefault,iso_eau500 #ifdef ISOVERIF501 USE isotopes_verif_mod502 #endif503 #ifdef ISOTRAC504 USE isotrac_mod, only: index_iso505 #endif506 #endif507 USE dimpft_mod_h508 USE flux_arp_mod_h509 USE compbl_mod_h510 USE yoethf_mod_h511 USE clesphys_mod_h512 USE ioipsl_getin_p_mod, ONLY : getin_p513 use phys_state_var_mod, only: ds_ns, dt_ns, delta_sst, delta_sal, dter, &514 dser, dt_ds, zsig, zmea, &515 frac_tersrf, z0m_tersrf, ratio_z0m_z0h_tersrf, albedo_tersrf !AM516 use phys_output_var_mod, only: tkt, tks, taur, sss517 use lmdz_blowing_snow_ini, only : zeta_bs518 use wxios_mod, ONLY: missing_val_xios => missing_val, using_xios519 USE netcdf, only: missing_val_netcdf => nf90_fill_real520 USE dimsoil_mod_h, ONLY: nsoilmx521 USE surf_param_mod, ONLY: eff_surf_param !AM522 523 USE yomcst_mod_h524 IMPLICIT NONE525 526 INCLUDE "FCTTRE.h"527 !FC528 529 !****************************************************************************************530 REAL, INTENT(IN) :: dtime ! time interval (s)531 REAL, INTENT(IN) :: date0 ! initial day532 INTEGER, INTENT(IN) :: itap ! time step533 INTEGER, INTENT(IN) :: jour ! current day of the year534 LOGICAL, INTENT(IN) :: debut ! true if first run step535 LOGICAL, INTENT(IN) :: lafin ! true if last run step536 REAL, DIMENSION(klon), INTENT(IN) :: rlon ! longitudes in degrees537 REAL, DIMENSION(klon), INTENT(IN) :: rlat ! latitudes in degrees538 REAL, DIMENSION(klon), INTENT(IN) :: rugoro ! rugosity length539 REAL, DIMENSION(klon), INTENT(IN) :: rmu0 ! cosine of solar zenith angle540 REAL, DIMENSION(klon), INTENT(IN) :: rain_f ! rain fall541 REAL, DIMENSION(klon), INTENT(IN) :: snow_f ! snow fall542 REAL, DIMENSION(klon), INTENT(IN) :: bs_f ! blowing snow fall543 REAL, DIMENSION(klon), INTENT(IN) :: solsw_m ! net shortwave radiation at mean surface544 REAL, DIMENSION(klon), INTENT(IN) :: solswfdiff_m ! diffuse fraction fordownward shortwave radiation at mean surface545 REAL, DIMENSION(klon), INTENT(IN) :: sollw_m ! net longwave radiation at mean surface546 REAL, DIMENSION(klon,klev), INTENT(IN) :: t ! temperature (K)547 REAL, DIMENSION(klon,klev), INTENT(IN) :: q ! water vapour (kg/kg)548 REAL, DIMENSION(klon,klev), INTENT(IN) :: qbs ! blowing snow specific content (kg/kg)549 REAL, DIMENSION(klon,klev), INTENT(IN) :: u ! u speed550 REAL, DIMENSION(klon,klev), INTENT(IN) :: v ! v speed551 REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! mid-layer pression (Pa)552 REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! pression between layers (Pa)553 REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf ! sub-surface fraction554 ! Martin555 REAL, DIMENSION(klon), INTENT(IN) :: lwdown_m ! downward longwave radiation at mean s556 REAL, DIMENSION(klon), INTENT(IN) :: gustiness ! gustiness557 558 !GG559 REAL, DIMENSION(klon,klev), INTENT(IN) :: pphi ! geopotential (m2/s2)560 !GG561 REAL, DIMENSION(klon), INTENT(IN) :: cldt ! total cloud562 563 #ifdef ISO564 REAL, DIMENSION(ntraciso,klon,klev), INTENT(IN) :: xt ! water vapour (kg/kg)565 REAL, DIMENSION(ntraciso,klon), INTENT(IN) :: xtrain_f ! rain fall566 REAL, DIMENSION(ntraciso,klon), INTENT(IN) :: xtsnow_f ! snow fall567 #endif568 569 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012570 !! REAL, DIMENSION(klon,klev), INTENT(IN) :: t_x ! Temp\'erature hors poche froide571 !! REAL, DIMENSION(klon,klev), INTENT(IN) :: t_w ! Temp\'erature dans la poches froide572 !! REAL, DIMENSION(klon,klev), INTENT(IN) :: q_x !573 !! REAL, DIMENSION(klon,klev), INTENT(IN) :: q_w ! Pareil pour l'humidit\'e574 REAL, DIMENSION(klon,klev), INTENT(IN) :: wake_dlt !temperature difference between (w) and (x) (K)575 REAL, DIMENSION(klon,klev), INTENT(IN) :: wake_dlq !humidity difference between (w) and (x) (K)576 REAL, DIMENSION(klon), INTENT(IN) :: wake_s ! Fraction de poches froides577 REAL, DIMENSION(klon), INTENT(IN) :: wake_cstar! Vitesse d'expansion des poches froides578 REAL, DIMENSION(klon), INTENT(IN) :: wake_dens579 !!!580 #ifdef ISO581 REAL, DIMENSION(ntraciso,klon,klev), INTENT(IN) :: wake_dlxt582 #endif583 ! Input/Output variables584 !****************************************************************************************585 !jyg<586 REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: beta ! Aridity factor587 !>jyg588 REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: ts ! temperature at surface (K)589 REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: delta_tsurf !surface temperature difference between590 !wake and off-wake regions591 !albedo SB >>>592 REAL, DIMENSIOn(6),intent(in) :: SFRWL593 REAL, DIMENSION(klon, nsw, nbsrf), INTENT(INOUT) :: alb_dir,alb_dif594 !albedo SB <<<595 !jyg Pourquoi ustar et wstar sont-elles INOUT ?596 REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: ustar ! u* (m/s)597 REAL, DIMENSION(klon, nbsrf+1), INTENT(INOUT) :: wstar ! w* (m/s)598 REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: u10m ! u speed at 10m599 REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: v10m ! v speed at 10m600 !jyg<601 !! REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke602 REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke_x603 !>jyg604 605 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012606 REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: wake_dltke ! TKE_w - TKE_x607 !!!608 609 ! Output variables610 !****************************************************************************************611 REAL, DIMENSION(klon,klev+1,nbsrf+1), INTENT(OUT) :: eps_x ! TKE dissipation rate612 613 REAL, DIMENSION(klon), INTENT(OUT) :: cdragh ! drag coefficient for T and Q614 REAL, DIMENSION(klon), INTENT(OUT) :: cdragm ! drag coefficient for wind615 REAL, DIMENSION(klon), INTENT(OUT) :: zu1 ! u wind speed in first layer616 REAL, DIMENSION(klon), INTENT(OUT) :: zv1 ! v wind speed in first layer617 !albedo SB >>>618 REAL, DIMENSION(klon, nsw), INTENT(OUT) :: alb_dir_m,alb_dif_m619 !albedo SB <<<620 ! Martin621 REAL, DIMENSION(klon), INTENT(OUT) :: alb3_lic622 ! Martin623 REAL, DIMENSION(klon), INTENT(OUT) :: zxsens ! sensible heat flux at surface with inversed sign624 ! (=> positive sign upwards)625 REAL, DIMENSION(klon), INTENT(OUT) :: zxevap ! water vapour flux at surface, positiv upwards626 REAL, DIMENSION(klon), INTENT(OUT) :: zxsnowerosion ! blowing snow flux at surface627 REAL, DIMENSION(klon), INTENT(OUT) :: icesub_lic ! ice (no snow!) sublimation over ice sheet628 REAL, DIMENSION(klon), INTENT(OUT) :: zxtsol ! temperature at surface, mean for each grid point629 !!! jyg le ???630 REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t_w ! !631 REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q_w ! ! Tendances dans les poches632 REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t_x ! !633 REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q_x ! ! Tendances hors des poches634 !!! jyg635 REAL, DIMENSION(klon), INTENT(OUT) :: zxfluxlat ! latent flux, mean for each grid point636 REAL, DIMENSION(klon), INTENT(OUT) :: zt2m ! temperature at 2m, mean for each grid point637 INTEGER, DIMENSION(klon, 6), INTENT(OUT) :: zn2mout ! number of times the 2m temperature is out of the [tsol,temp]638 REAL, DIMENSION(klon), INTENT(OUT) :: qsat2m639 REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_t ! change in temperature640 REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_t_diss ! change in temperature641 REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_q ! change in water vapour642 REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_u ! change in u speed643 REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_v ! change in v speed644 REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_qbs ! change in blowing snow specific content645 646 647 REAL, INTENT(OUT):: zcoefh(klon, klev+1, nbsrf + 1) ! (klon, klev, nbsrf + 1) => only use (klon, klev, nbsrf+1)648 ! coef for turbulent diffusion of T and Q, mean for each grid point649 650 REAL, INTENT(OUT):: zcoefm(klon, klev+1, nbsrf + 1) ! (klon, klev, nbsrf + 1) => only use (klon, klev, nbsrf+1)651 ! coef for turbulent diffusion of U and V (?), mean for each grid point652 #ifdef ISO653 REAL, DIMENSION(ntraciso,klon), INTENT(OUT) :: zxxtevap ! water vapour flux at surface, positiv upwards654 REAL, DIMENSION(ntraciso,klon, klev), INTENT(OUT) :: d_xt ! change in water vapour655 REAL, DIMENSION(klon), INTENT(OUT) :: runoff_diag656 REAL, DIMENSION(niso,klon), INTENT(OUT) :: xtrunoff_diag657 REAL, DIMENSION(ntraciso,klon,klev), INTENT(OUT) :: d_xt_w658 REAL, DIMENSION(ntraciso,klon,klev), INTENT(OUT) :: d_xt_x659 #endif660 661 662 663 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012664 REAL, DIMENSION(klon), INTENT(OUT) :: zxsens_x ! Flux sensible hors poche665 REAL, DIMENSION(klon), INTENT(OUT) :: zxsens_w ! Flux sensible dans la poche666 REAL, DIMENSION(klon), INTENT(OUT) :: zxfluxlat_x! Flux latent hors poche667 REAL, DIMENSION(klon), INTENT(OUT) :: zxfluxlat_w! Flux latent dans la poche668 !! REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_wake_dlt669 !! REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_wake_dlq670 671 ! Output only for diagnostics672 REAL, DIMENSION(klon), INTENT(OUT) :: cdragh_x673 REAL, DIMENSION(klon), INTENT(OUT) :: cdragh_w674 REAL, DIMENSION(klon), INTENT(OUT) :: cdragm_x675 REAL, DIMENSION(klon), INTENT(OUT) :: cdragm_w676 REAL, DIMENSION(klon), INTENT(OUT) :: kh677 REAL, DIMENSION(klon), INTENT(OUT) :: kh_x678 REAL, DIMENSION(klon), INTENT(OUT) :: kh_w679 !!!680 REAL, DIMENSION(klon), INTENT(OUT) :: slab_wfbils! heat balance at surface only for slab at ocean points681 REAL, DIMENSION(klon), INTENT(OUT) :: qsol ! water height in the soil (mm)682 REAL, DIMENSION(klon), INTENT(OUT) :: zq2m ! water vapour at 2m, mean for each grid point683 REAL, DIMENSION(klon), INTENT(OUT) :: s_pblh ! height of the planetary boundary layer(HPBL)684 !!! jyg le 08/02/2012685 REAL, DIMENSION(klon), INTENT(OUT) :: s_pblh_x ! height of the PBL in the off-wake region686 REAL, DIMENSION(klon), INTENT(OUT) :: s_pblh_w ! height of the PBL in the wake region687 !!!688 REAL, DIMENSION(klon), INTENT(OUT) :: s_plcl ! condensation level689 !!! jyg le 08/02/2012690 REAL, DIMENSION(klon), INTENT(OUT) :: s_plcl_x ! condensation level in the off-wake region691 REAL, DIMENSION(klon), INTENT(OUT) :: s_plcl_w ! condensation level in the wake region692 !!!693 REAL, DIMENSION(klon), INTENT(OUT) :: s_capCL ! CAPE of PBL694 REAL, DIMENSION(klon), INTENT(OUT) :: s_oliqCL ! liquid water intergral of PBL695 REAL, DIMENSION(klon), INTENT(OUT) :: s_cteiCL ! cloud top instab. crit. of PBL696 REAL, DIMENSION(klon), INTENT(OUT) :: s_pblT ! temperature at PBLH697 REAL, DIMENSION(klon), INTENT(OUT) :: s_therm ! thermal virtual temperature excess698 REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb1 ! deep cape, mean for each grid point699 REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb2 ! inhibition, mean for each grid point700 REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb3 ! point Omega, mean for each grid point701 REAL, DIMENSION(klon), INTENT(OUT) :: zustar ! u*702 REAL, DIMENSION(klon), INTENT(OUT) :: zu10m ! u speed at 10m, mean for each grid point703 REAL, DIMENSION(klon), INTENT(OUT) :: zv10m ! v speed at 10m, mean for each grid point704 REAL, DIMENSION(klon), INTENT(OUT) :: fder_print ! fder for printing (=fder(i) + dflux_t(i) + dflux_q(i))705 REAL, DIMENSION(klon), INTENT(OUT) :: zxqsurf ! humidity at surface, mean for each grid point706 REAL, DIMENSION(klon), INTENT(OUT) :: delta_qsurf! humidity difference at surface, mean for each grid point707 REAL, DIMENSION(klon), INTENT(OUT) :: rh2m ! relative humidity at 2m708 REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxu ! u wind tension, mean for each grid point709 REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxv ! v wind tension, mean for each grid point710 REAL, DIMENSION(klon, nbsrf+1), INTENT(INOUT) :: z0m,z0h ! rugosity length (m)711 REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: agesno ! age of snow at surface712 REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: solsw ! net shortwave radiation at surface713 REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: sollw ! net longwave radiation at surface714 REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: d_ts ! change in temperature at surface715 REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: evap ! evaporation at surface716 REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: fluxlat ! latent flux717 REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: t2m ! temperature at 2 meter height718 REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfbils ! heat balance at surface719 REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfevap ! water balance (evap) at surface weighted by srf720 REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_t ! sensible heat flux (CpT) J/m**2/s (W/m**2)721 ! positve orientation downwards722 REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_u ! u wind tension (kg m/s)/(m**2 s) or Pascal723 REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_v ! v wind tension (kg m/s)/(m**2 s) or Pascal724 !FC725 REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: treedrg ! tree drag (m)726 !AM heterogeneous continental sub-surfaces727 REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: tsurf_tersrf ! surface temperature of continental sub-surfaces (K)728 REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: qsurf_tersrf ! surface specific humidity of continental sub-surfaces (kg/kg)729 REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: tsurf_new_tersrf ! surface temperature of continental sub-surfaces (K)730 REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: cdragm_tersrf ! momentum drag coefficient of continental sub-surfaces (-)731 REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: cdragh_tersrf ! heat drag coefficient of continental sub-surfaces (-)732 REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: swnet_tersrf ! net shortwave radiation of continental sub-surfaces (W/m2)733 REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: lwnet_tersrf ! net longwave radiation of continental sub-surfaces (W/m2)734 REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: fluxsens_tersrf ! sensible heat flux of continental sub-surfaces (W/m2)735 REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: fluxlat_tersrf ! latent heat flux of continental sub-surfaces (W/m2)736 REAL, DIMENSION(klon, nsoilmx, nbtersrf), INTENT(INOUT) :: tsoil_tersrf ! soil temperature of continental sub-surfaces (K)737 #ifdef ISO738 REAL, DIMENSION(niso,klon), INTENT(OUT) :: xtsol ! water height in the soil (mm)739 REAL, DIMENSION(ntraciso,klon, nbsrf) :: xtevap ! evaporation at surface740 REAL, DIMENSION(klon), INTENT(OUT) :: h1_diag ! just diagnostic, not useful741 #endif742 743 744 ! Output not needed745 REAL, DIMENSION(klon), INTENT(OUT) :: dflux_t ! change of sensible heat flux746 REAL, DIMENSION(klon), INTENT(OUT) :: dflux_q ! change of water vapour flux747 REAL, DIMENSION(klon), INTENT(OUT) :: zxsnow ! snow at surface, mean for each grid point748 REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxt ! sensible heat flux, mean for each grid point749 REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxq ! water vapour flux, mean for each grid point750 REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxqbs ! blowing snow flux, mean for each grid point751 REAL, DIMENSION(klon, nbsrf),INTENT(OUT) :: q2m ! water vapour at 2 meter height752 REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_q ! water vapour flux(latent flux) (kg/m**2/s)753 REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_qbs ! blowind snow vertical flux (kg/m**2754 755 #ifdef ISO756 REAL, DIMENSION(ntraciso,klon), INTENT(OUT) :: dflux_xt ! change of water vapour flux757 REAL, DIMENSION(niso,klon), INTENT(OUT) :: zxxtsnow ! snow at surface, mean for each grid point758 REAL, DIMENSION(ntraciso,klon, klev), INTENT(OUT) :: zxfluxxt ! water vapour flux, mean for each grid point759 REAL, DIMENSION(ntraciso,klon, klev, nbsrf), INTENT(OUT) :: flux_xt ! water vapour flux(latent flux) (kg/m**2/s)760 #endif761 762 ! Martin763 ! inlandsis764 REAL, DIMENSION(klon), INTENT(OUT) :: qsnow ! snow water content765 REAL, DIMENSION(klon), INTENT(OUT) :: snowhgt ! snow height766 REAL, DIMENSION(klon), INTENT(OUT) :: to_ice ! snow passed to ice767 REAL, DIMENSION(klon), INTENT(OUT) :: sissnow ! snow in snow model768 REAL, DIMENSION(klon), INTENT(OUT) :: runoff ! runoff on land ice769 ! Martin770 !GG771 REAL, DIMENSION(klon), INTENT(INOUT) :: hice ! hice772 REAL, DIMENSION(klon), INTENT(INOUT) :: tice ! tice773 REAL, DIMENSION(klon), INTENT(INOUT) :: bilg_cumul ! flux cumulated774 REAL, DIMENSION(klon), INTENT(INOUT) :: fcds775 REAL, DIMENSION(klon), INTENT(INOUT) :: fcdi776 REAL, DIMENSION(klon), INTENT(INOUT) :: dh_basal_growth777 REAL, DIMENSION(klon), INTENT(INOUT) :: dh_basal_melt778 REAL, DIMENSION(klon), INTENT(INOUT) :: dh_top_melt779 REAL, DIMENSION(klon), INTENT(INOUT) :: dh_snow2sic780 REAL, DIMENSION(klon), INTENT(INOUT) :: dtice_melt781 REAL, DIMENSION(klon), INTENT(INOUT) :: dtice_snow2sic782 !GG783 784 ! Other local variables785 !****************************************************************************************786 ! >> PC787 INTEGER :: ierr788 INTEGER :: n789 ! << PC790 INTEGER :: iflag_split, iflag_split_ref791 INTEGER :: i, k, nsrf792 INTEGER :: knon, j793 INTEGER :: idayref794 INTEGER , DIMENSION(klon) :: ni795 REAL :: yt1_new796 REAL :: zx_alf1, zx_alf2 !valeur ambiante par extrapola797 REAL :: amn, amx798 REAL :: f1 ! fraction de longeurs visibles parmi tout SW intervalle799 REAL, DIMENSION(klon) :: r_co2_ppm ! taux CO2 atmosphere800 REAL, DIMENSION(klon) :: yts, yz0m, yz0h, ypct801 REAL, DIMENSION(klon) :: yz0h_old802 !albedo SB >>>803 REAL, DIMENSION(klon) :: yalb,yalb_vis804 !albedo SB <<<805 REAL, DIMENSION(klon) :: yt1, yq1, yu1, yv1, yqbs1806 REAL, DIMENSION(klon) :: yqa807 REAL, DIMENSION(klon) :: ysnow, yqsurf, yagesno, yqsol808 REAL, DIMENSION(klon) :: yrain_f, ysnow_f, ybs_f809 #ifdef ISO810 REAL, DIMENSION(ntraciso,klon) :: yxt1811 REAL, DIMENSION(niso,klon) :: yxtsnow, yxtsol812 REAL, DIMENSION(ntraciso,klon) :: yxtrain_f, yxtsnow_f813 REAL, DIMENSION(klon) :: yrunoff_diag814 REAL, DIMENSION(niso,klon) :: yxtrunoff_diag815 REAL, DIMENSION(niso,klon) :: yRland_ice816 #endif817 REAL, DIMENSION(klon) :: ysolsw, ysollw818 REAL, DIMENSION(klon) :: yfder819 REAL, DIMENSION(klon) :: yrugoro820 REAL, DIMENSION(klon) :: yfluxlat821 REAL, DIMENSION(klon) :: yfluxbs822 REAL, DIMENSION(klon) :: y_d_ts823 REAL, DIMENSION(klon) :: y_flux_t1, y_flux_q1824 REAL, DIMENSION(klon) :: y_dflux_t, y_dflux_q825 #ifdef ISO826 REAL, DIMENSION(ntraciso,klon) :: y_flux_xt1827 REAL, DIMENSION(ntraciso,klon) :: y_dflux_xt828 #endif829 REAL, DIMENSION(klon) :: y_flux_u1, y_flux_v1830 REAL, DIMENSION(klon) :: y_flux_bs, y_flux0831 REAL, DIMENSION(klon) :: yt2m, yq2m, yu10m832 INTEGER, DIMENSION(klon, nbsrf, 6) :: yn2mout, yn2mout_x, yn2mout_w833 INTEGER, DIMENSION(klon, nbsrf, 6) :: n2mout, n2mout_x, n2mout_w834 REAL, DIMENSION(klon) :: yustar835 REAL, DIMENSION(klon) :: ywstar836 REAL, DIMENSION(klon) :: ywindsp837 REAL, DIMENSION(klon) :: yt10m, yq10m838 REAL, DIMENSION(klon) :: ypblh839 REAL, DIMENSION(klon) :: ylcl840 REAL, DIMENSION(klon) :: ycapCL841 REAL, DIMENSION(klon) :: yoliqCL842 REAL, DIMENSION(klon) :: ycteiCL843 REAL, DIMENSION(klon) :: ypblT844 REAL, DIMENSION(klon) :: ytherm845 REAL, DIMENSION(klon) :: ytrmb1846 REAL, DIMENSION(klon) :: ytrmb2847 REAL, DIMENSION(klon) :: ytrmb3848 REAL, DIMENSION(klon) :: uzon, vmer849 REAL, DIMENSION(klon) :: tair1, qair1, tairsol850 REAL, DIMENSION(klon) :: psfce, patm851 REAL, DIMENSION(klon) :: qairsol, zgeo1, speed, zri1, pref !speed, zri1, pref, added by Fuxing WANG, 04/03/2015852 REAL, DIMENSION(klon) :: yz0h_oupas853 REAL, DIMENSION(klon) :: yfluxsens854 REAL, DIMENSION(klon) :: AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0855 REAL, DIMENSION(klon) :: AcoefH, AcoefQ, BcoefH, BcoefQ856 #ifdef ISO857 REAL, DIMENSION(ntraciso,klon) :: AcoefXT, BcoefXT858 #endif859 REAL, DIMENSION(klon) :: AcoefU, AcoefV, BcoefU, BcoefV860 REAL, DIMENSION(klon) :: AcoefQBS, BcoefQBS861 REAL, DIMENSION(klon) :: ypsref862 REAL, DIMENSION(klon) :: yevap, yevap_pot, ytsurf_new, yalb3_new, yicesub_lic863 !albedo SB >>>864 REAL, DIMENSION(klon,nsw) :: yalb_dir_new, yalb_dif_new865 !albedo SB <<<866 REAL, DIMENSION(klon) :: ztsol867 REAL, DIMENSION(klon) :: meansqT ! mean square deviation of subsurface temperatures868 REAL, DIMENSION(klon) :: alb_m ! mean albedo for whole SW interval869 REAL, DIMENSION(klon,klev) :: y_d_t, y_d_q, y_d_t_diss, y_d_qbs870 REAL, DIMENSION(klon,klev) :: y_d_u, y_d_v871 REAL, DIMENSION(klon,klev) :: y_flux_t, y_flux_q, y_flux_qbs872 REAL, DIMENSION(klon,klev) :: y_flux_u, y_flux_v873 REAL, DIMENSION(klon,klev) :: ycoefh,ycoefm,ycoefq,ycoefqbs874 REAL, DIMENSION(klon) :: ycdragh, ycdragq, ycdragm875 REAL, DIMENSION(klon,klev) :: yu, yv876 REAL, DIMENSION(klon,klev) :: yt, yq, yqbs877 #ifdef ISO878 REAL, DIMENSION(ntraciso,klon) :: yxtevap879 REAL, DIMENSION(ntraciso,klon,klev) :: y_d_xt880 REAL, DIMENSION(ntraciso,klon,klev) :: y_flux_xt881 REAL, DIMENSION(ntraciso,klon,klev) :: yxt882 #endif883 REAL, DIMENSION(klon,klev) :: ypplay, ydelp884 REAL, DIMENSION(klon,klev) :: delp885 REAL, DIMENSION(klon,klev+1) :: ypaprs886 REAL, DIMENSION(klon,klev+1) :: ytke, yeps887 REAL, DIMENSION(klon,nsoilmx) :: ytsoil888 !FC889 REAL, DIMENSION(klon,nvm_lmdz) :: yveget890 REAL, DIMENSION(klon,nvm_lmdz) :: ylai891 REAL, DIMENSION(klon,nvm_lmdz) :: yheight892 REAL, DIMENSION(klon,klev) :: y_d_u_frein893 REAL, DIMENSION(klon,klev) :: y_d_v_frein894 REAL, DIMENSION(klon,klev) :: y_treedrg895 !FC896 897 CHARACTER(len=80) :: abort_message898 CHARACTER(len=20) :: modname = 'pbl_surface'899 LOGICAL, PARAMETER :: zxli=.FALSE. ! utiliser un jeu de fonctions simples900 LOGICAL, PARAMETER :: check=.FALSE.901 902 !!! nrlmd le 02/05/2011903 !!! jyg le 07/02/2012904 REAL, DIMENSION(klon) :: ywake_s, ywake_cstar, ywake_dens905 !!!906 REAL, DIMENSION(klon,klev+1) :: ytke_x, ytke_w, yeps_x, yeps_w907 REAL, DIMENSION(klon,klev+1) :: ywake_dltke908 REAL, DIMENSION(klon,klev) :: yu_x, yv_x, yu_w, yv_w909 REAL, DIMENSION(klon,klev) :: yt_x, yq_x, yt_w, yq_w910 REAL, DIMENSION(klon,klev) :: ycoefh_x, ycoefm_x, ycoefh_w, ycoefm_w911 REAL, DIMENSION(klon,klev) :: ycoefq_x, ycoefq_w912 REAL, DIMENSION(klon) :: ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w913 REAL, DIMENSION(klon) :: ycdragm_x, ycdragm_w914 REAL, DIMENSION(klon) :: AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x915 REAL, DIMENSION(klon) :: AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w916 REAL, DIMENSION(klon) :: AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x917 REAL, DIMENSION(klon) :: AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w918 REAL, DIMENSION(klon) :: y_flux_t1_x, y_flux_q1_x, y_flux_t1_w, y_flux_q1_w919 REAL, DIMENSION(klon) :: y_flux_u1_x, y_flux_v1_x, y_flux_u1_w, y_flux_v1_w920 REAL, DIMENSION(klon,klev) :: y_flux_t_x, y_flux_q_x, y_flux_t_w, y_flux_q_w921 REAL, DIMENSION(klon,klev) :: y_flux_u_x, y_flux_v_x, y_flux_u_w, y_flux_v_w922 REAL, DIMENSION(klon) :: yfluxlat_x, yfluxlat_w923 REAL, DIMENSION(klon,klev) :: y_d_t_x, y_d_q_x, y_d_t_w, y_d_q_w924 REAL, DIMENSION(klon,klev) :: y_d_t_diss_x, y_d_t_diss_w925 REAL, DIMENSION(klon,klev) :: d_t_diss_x, d_t_diss_w926 REAL, DIMENSION(klon,klev) :: y_d_u_x, y_d_v_x, y_d_u_w, y_d_v_w927 REAL, DIMENSION(klon, klev, nbsrf) :: flux_t_x, flux_q_x, flux_t_w, flux_q_w928 REAL, DIMENSION(klon, klev, nbsrf) :: flux_u_x, flux_v_x, flux_u_w, flux_v_w929 REAL, DIMENSION(klon, nbsrf) :: fluxlat_x, fluxlat_w930 REAL, DIMENSION(klon, klev) :: zxfluxt_x, zxfluxq_x, zxfluxt_w, zxfluxq_w931 REAL, DIMENSION(klon, klev) :: zxfluxu_x, zxfluxv_x, zxfluxu_w, zxfluxv_w932 REAL :: zx_qs_surf, zcor_surf, zdelta_surf933 !jyg<934 REAL, DIMENSION(klon) :: ybeta935 REAL, DIMENSION(klon) :: ybeta_prev936 !>jyg937 REAL, DIMENSION(klon, klev) :: d_u_x938 REAL, DIMENSION(klon, klev) :: d_u_w939 REAL, DIMENSION(klon, klev) :: d_v_x940 REAL, DIMENSION(klon, klev) :: d_v_w941 942 REAL, DIMENSION(klon,klev) :: CcoefH, CcoefQ, DcoefH, DcoefQ943 REAL, DIMENSION(klon,klev) :: CcoefU, CcoefV, DcoefU, DcoefV944 REAL, DIMENSION(klon,klev) :: CcoefQBS, DcoefQBS945 REAL, DIMENSION(klon,klev) :: CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x946 REAL, DIMENSION(klon,klev) :: CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w947 REAL, DIMENSION(klon,klev) :: CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x948 REAL, DIMENSION(klon,klev) :: CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w949 REAL, DIMENSION(klon,klev) :: Kcoef_hq, Kcoef_m, gama_h, gama_q950 REAL, DIMENSION(klon,klev) :: gama_qbs, Kcoef_qbs951 REAL, DIMENSION(klon,klev) :: Kcoef_hq_x, Kcoef_m_x, gama_h_x, gama_q_x952 REAL, DIMENSION(klon,klev) :: Kcoef_hq_w, Kcoef_m_w, gama_h_w, gama_q_w953 REAL, DIMENSION(klon) :: alf_1, alf_2, alf_1_x, alf_2_x, alf_1_w, alf_2_w954 #ifdef ISO955 REAL, DIMENSION(ntraciso,klon,klev) :: yxt_x, yxt_w956 REAL, DIMENSION(ntraciso,klon) :: y_flux_xt1_x , y_flux_xt1_w957 REAL, DIMENSION(ntraciso,klon,klev) :: y_flux_xt_x,y_d_xt_x,zxfluxxt_x958 REAL, DIMENSION(ntraciso,klon,klev) :: y_flux_xt_w,y_d_xt_w,zxfluxxt_w959 REAL, DIMENSION(ntraciso,klon,klev,nbsrf) :: flux_xt_x, flux_xt_w960 REAL, DIMENSION(ntraciso,klon) :: AcoefXT_x, BcoefXT_x961 REAL, DIMENSION(ntraciso,klon) :: AcoefXT_w, BcoefXT_w962 REAL, DIMENSION(ntraciso,klon,klev) :: CcoefXT, DcoefXT963 REAL, DIMENSION(ntraciso,klon,klev) :: CcoefXT_x, DcoefXT_x964 REAL, DIMENSION(ntraciso,klon,klev) :: CcoefXT_w, DcoefXT_w965 REAL, DIMENSION(ntraciso,klon,klev) :: gama_xt,gama_xt_x,gama_xt_w966 #endif967 !!!968 !!!jyg le 08/02/2012969 REAL, DIMENSION(klon, nbsrf) :: windsp970 !971 REAL, DIMENSION(klon, nbsrf) :: t2m_x972 REAL, DIMENSION(klon, nbsrf) :: q2m_x973 REAL, DIMENSION(klon) :: rh2m_x974 REAL, DIMENSION(klon) :: qsat2m_x975 REAL, DIMENSION(klon, nbsrf) :: u10m_x976 REAL, DIMENSION(klon, nbsrf) :: v10m_x977 REAL, DIMENSION(klon, nbsrf) :: ustar_x978 REAL, DIMENSION(klon, nbsrf) :: wstar_x979 !980 REAL, DIMENSION(klon, nbsrf) :: pblh_x981 REAL, DIMENSION(klon, nbsrf) :: plcl_x982 REAL, DIMENSION(klon, nbsrf) :: capCL_x983 REAL, DIMENSION(klon, nbsrf) :: oliqCL_x984 REAL, DIMENSION(klon, nbsrf) :: cteiCL_x985 REAL, DIMENSION(klon, nbsrf) :: pblt_x986 REAL, DIMENSION(klon, nbsrf) :: therm_x987 REAL, DIMENSION(klon, nbsrf) :: trmb1_x988 REAL, DIMENSION(klon, nbsrf) :: trmb2_x989 REAL, DIMENSION(klon, nbsrf) :: trmb3_x990 !991 REAL, DIMENSION(klon, nbsrf) :: t2m_w992 REAL, DIMENSION(klon, nbsrf) :: q2m_w993 REAL, DIMENSION(klon) :: rh2m_w994 REAL, DIMENSION(klon) :: qsat2m_w995 REAL, DIMENSION(klon, nbsrf) :: u10m_w996 REAL, DIMENSION(klon, nbsrf) :: v10m_w997 REAL, DIMENSION(klon, nbsrf) :: ustar_w998 REAL, DIMENSION(klon, nbsrf) :: wstar_w999 !1000 REAL, DIMENSION(klon, nbsrf) :: pblh_w1001 REAL, DIMENSION(klon, nbsrf) :: plcl_w1002 REAL, DIMENSION(klon, nbsrf) :: capCL_w1003 REAL, DIMENSION(klon, nbsrf) :: oliqCL_w1004 REAL, DIMENSION(klon, nbsrf) :: cteiCL_w1005 REAL, DIMENSION(klon, nbsrf) :: pblt_w1006 REAL, DIMENSION(klon, nbsrf) :: therm_w1007 REAL, DIMENSION(klon, nbsrf) :: trmb1_w1008 REAL, DIMENSION(klon, nbsrf) :: trmb2_w1009 REAL, DIMENSION(klon, nbsrf) :: trmb3_w1010 !1011 REAL, DIMENSION(klon) :: yt2m_x1012 REAL, DIMENSION(klon) :: yq2m_x1013 REAL, DIMENSION(klon) :: yt10m_x1014 REAL, DIMENSION(klon) :: yq10m_x1015 REAL, DIMENSION(klon) :: yu10m_x1016 REAL, DIMENSION(klon) :: yv10m_x1017 REAL, DIMENSION(klon) :: yustar_x1018 REAL, DIMENSION(klon) :: ywstar_x1019 !1020 REAL, DIMENSION(klon) :: ypblh_x1021 REAL, DIMENSION(klon) :: ylcl_x1022 REAL, DIMENSION(klon) :: ycapCL_x1023 REAL, DIMENSION(klon) :: yoliqCL_x1024 REAL, DIMENSION(klon) :: ycteiCL_x1025 REAL, DIMENSION(klon) :: ypblt_x1026 REAL, DIMENSION(klon) :: ytherm_x1027 REAL, DIMENSION(klon) :: ytrmb1_x1028 REAL, DIMENSION(klon) :: ytrmb2_x1029 REAL, DIMENSION(klon) :: ytrmb3_x1030 !1031 REAL, DIMENSION(klon) :: yt2m_w1032 REAL, DIMENSION(klon) :: yq2m_w1033 REAL, DIMENSION(klon) :: yt10m_w1034 REAL, DIMENSION(klon) :: yq10m_w1035 REAL, DIMENSION(klon) :: yu10m_w1036 REAL, DIMENSION(klon) :: yv10m_w1037 REAL, DIMENSION(klon) :: yustar_w1038 REAL, DIMENSION(klon) :: ywstar_w1039 !1040 REAL, DIMENSION(klon) :: ypblh_w1041 REAL, DIMENSION(klon) :: ylcl_w1042 REAL, DIMENSION(klon) :: ycapCL_w1043 REAL, DIMENSION(klon) :: yoliqCL_w1044 REAL, DIMENSION(klon) :: ycteiCL_w1045 REAL, DIMENSION(klon) :: ypblt_w1046 REAL, DIMENSION(klon) :: ytherm_w1047 REAL, DIMENSION(klon) :: ytrmb1_w1048 REAL, DIMENSION(klon) :: ytrmb2_w1049 REAL, DIMENSION(klon) :: ytrmb3_w1050 !1051 REAL, DIMENSION(klon) :: uzon_x, vmer_x, speed_x, zri1_x, pref_x !speed_x, zri1_x, pref_x, added by Fuxing WANG, 04/03/20151052 REAL, DIMENSION(klon) :: zgeo1_x, tair1_x, qair1_x, tairsol_x1053 !1054 REAL, DIMENSION(klon) :: uzon_w, vmer_w, speed_w, zri1_w, pref_w !speed_w, zri1_w, pref_w, added by Fuxing WANG, 04/03/20151055 REAL, DIMENSION(klon) :: zgeo1_w, tair1_w, qair1_w, tairsol_w1056 REAL, DIMENSION(klon) :: yus0, yvs01057 1058 !!! jyg le 25/03/20131059 !! Variables intermediaires pour le raccord des deux colonnes \`a la surface1060 !jyg<1061 !! REAL :: dd_Ch1062 !! REAL :: dd_Cm1063 !! REAL :: dd_Kh1064 !! REAL :: dd_Km1065 !! REAL :: dd_u1066 !! REAL :: dd_v1067 !! REAL :: dd_t1068 !! REAL :: dd_q1069 !! REAL :: dd_AH1070 !! REAL :: dd_AQ1071 !! REAL :: dd_AU1072 !! REAL :: dd_AV1073 !! REAL :: dd_BH1074 !! REAL :: dd_BQ1075 !! REAL :: dd_BU1076 !! REAL :: dd_BV1077 !!1078 !! REAL :: dd_KHp1079 !! REAL :: dd_KQp1080 !! REAL :: dd_KUp1081 !! REAL :: dd_KVp1082 !>jyg1083 1084 !!!1085 !!! nrlmd le 13/06/20111086 REAL, DIMENSION(klon) :: y_delta_flux_t1, y_delta_flux_q1, y_delta_flux_u1, y_delta_flux_v11087 REAL, DIMENSION(klon) :: y_delta_tsurf, y_delta_tsurf_new1088 REAL, DIMENSION(klon) :: delta_coef, tau_eq1089 REAL, DIMENSION(klon) :: HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn1090 REAL, DIMENSION(klon) :: phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn01091 REAL, DIMENSION(klon) :: y_delta_qsurf1092 REAL, DIMENSION(klon) :: y_delta_qsats1093 REAL, DIMENSION(klon) :: yg_T, yg_Q1094 REAL, DIMENSION(klon) :: yGamma_dTs_phiT, yGamma_dQs_phiQ1095 REAL, DIMENSION(klon) :: ydTs_ins, ydqs_ins1096 !1097 REAL, PARAMETER :: facteur=2./sqrt(3.14)1098 REAL, PARAMETER :: inertia=2000.1099 REAL, DIMENSION(klon) :: ydtsurf_th1100 REAL :: zdelta_surf_x,zdelta_surf_w,zx_qs_surf_x,zx_qs_surf_w1101 REAL :: zcor_surf_x,zcor_surf_w1102 REAL :: mod_wind_x, mod_wind_w1103 REAL :: rho11104 REAL, DIMENSION(klon) :: Kech_h ! Coefficient d'echange pour l'energie1105 REAL, DIMENSION(klon) :: Kech_h_x, Kech_h_w1106 REAL, DIMENSION(klon) :: Kech_m1107 REAL, DIMENSION(klon) :: Kech_m_x, Kech_m_w1108 REAL, DIMENSION(klon) :: yts_x, yts_w1109 REAL, DIMENSION(klon) :: yqsatsrf0_x, yqsatsrf0_w1110 REAL, DIMENSION(klon) :: yqsurf_x, yqsurf_w1111 !jyg<1112 !! REAL, DIMENSION(klon) :: Kech_Hp, Kech_H_xp, Kech_H_wp1113 !! REAL, DIMENSION(klon) :: Kech_Qp, Kech_Q_xp, Kech_Q_wp1114 !! REAL, DIMENSION(klon) :: Kech_Up, Kech_U_xp, Kech_U_wp1115 !! REAL, DIMENSION(klon) :: Kech_Vp, Kech_V_xp, Kech_V_wp1116 !>jyg1117 1118 REAL :: fact_cdrag1119 REAL :: z1lay1120 1121 REAL :: vent1122 !1123 ! For debugging with IOIPSL1124 INTEGER, DIMENSION(nbp_lon*nbp_lat) :: ndexbg1125 REAL :: zjulian1126 REAL, DIMENSION(klon) :: tabindx1127 REAL, DIMENSION(nbp_lon,nbp_lat) :: zx_lon, zx_lat1128 REAL, DIMENSION(nbp_lon,nbp_lat) :: debugtab1129 1130 1131 REAL, DIMENSION(klon,nbsrf) :: pblh ! height of the planetary boundary layer1132 REAL, DIMENSION(klon,nbsrf) :: plcl ! condensation level1133 REAL, DIMENSION(klon,nbsrf) :: capCL1134 REAL, DIMENSION(klon,nbsrf) :: oliqCL1135 REAL, DIMENSION(klon,nbsrf) :: cteiCL1136 REAL, DIMENSION(klon,nbsrf) :: pblT1137 REAL, DIMENSION(klon,nbsrf) :: therm1138 REAL, DIMENSION(klon,nbsrf) :: trmb1 ! deep cape1139 REAL, DIMENSION(klon,nbsrf) :: trmb2 ! inhibition1140 REAL, DIMENSION(klon,nbsrf) :: trmb3 ! point Omega1141 REAL, DIMENSION(klon,nbsrf) :: zx_rh2m, zx_qsat2m1142 REAL, DIMENSION(klon,nbsrf) :: zx_t11143 REAL, DIMENSION(klon, nbsrf) :: alb ! mean albedo for whole SW interval1144 REAL, DIMENSION(klon,nbsrf) :: snowerosion1145 REAL, DIMENSION(klon) :: ylwdown ! jg : temporary (ysollwdown)1146 REAL, DIMENSION(klon) :: ygustiness ! jg : temporary (ysollwdown)1147 1148 REAL :: zx_qs1, zcor1, zdelta11149 1150 ! Martin1151 REAL, DIMENSION(klon, nbsrf) :: sollwd ! net longwave radiation at surface1152 REAL, DIMENSION(klon) :: ytoice1153 REAL, DIMENSION(klon) :: ysnowhgt, yqsnow, ysissnow, yrunoff1154 REAL, DIMENSION(klon) :: yzmea1155 REAL, DIMENSION(klon) :: yzsig1156 REAL, DIMENSION(klon) :: ycldt1157 REAL, DIMENSION(klon) :: yrmu01158 ! Martin1159 REAL, DIMENSION(klon) :: yri01160 1161 REAL, DIMENSION(klon):: ydelta_sst, ydelta_sal, yds_ns, ydt_ns, ydter, &1162 ydser, ydt_ds, ytkt, ytks, ytaur, ysss1163 ! compression of delta_sst, delta_sal, ds_ns, dt_ns, dter, dser,1164 ! dt_ds, tkt, tks, taur, sss on ocean points1165 REAL :: missing_val1166 1167 ! GG1168 REAL, DIMENSION(klon,klev) :: ytheta1169 REAL, DIMENSION(klon,klev) :: ypphii1170 REAL, DIMENSION(klon,klev) :: ypphi1171 REAL, DIMENSION(klon,klev) :: ydthetadz1172 REAL, DIMENSION(klon) :: ydthetadz3001173 REAL, DIMENSION(klon) :: Ampl1174 ! GG1175 1176 ! AM !1177 REAL, DIMENSION(klon) :: albedo_eff1178 REAL, DIMENSION(klon, nbtersrf) :: yfrac_tersrf, yz0m_tersrf, yz0h_tersrf1179 #ifdef ISO1180 REAL, DIMENSION(klon) :: h11181 INTEGER :: ixt1182 !#ifdef ISOVERIF1183 ! integer iso_verif_positif_nostop1184 !#endif1185 #endif1186 1187 !****************************************************************************************1188 ! End of declarations1189 !****************************************************************************************1190 IF (using_xios) THEN1191 missing_val=missing_val_xios1192 ELSE1193 missing_val=missing_val_netcdf1194 ENDIF1195 1196 IF (prt_level >=10) print *,' -> pbl_surface, itap ',itap1197 !1198 !!jyg iflag_split = mod(iflag_pbl_split,2)1199 !!jyg iflag_split = mod(iflag_pbl_split,10)1200 !1201 ! Flags controlling the splitting of the turbulent boundary layer:1202 ! iflag_split_ref = 0 ==> no splitting1203 ! = 1 ==> splitting without coupling with surface temperature1204 ! = 2 ==> splitting with coupling with surface temperature over land1205 ! = 3 ==> splitting over ocean; no splitting over land1206 ! iflag_split: actual flag controlling the splitting.1207 ! iflag_split = iflag_split_ref outside the sub-surface loop1208 ! = iflag_split_ref if iflag_split_ref = 0, 1, or 21209 ! = 0 over land if iflga_split_ref = 31210 ! = 1 over ocean if iflga_split_ref = 31211 1212 iflag_split_ref = mod(iflag_pbl_split,10)1213 iflag_split = iflag_split_ref1214 1215 #ifdef ISO1216 #ifdef ISOVERIF1217 DO i=1,klon1218 DO ixt=1,niso1219 CALL iso_verif_noNaN(xtsol(ixt,i),'pbl_surface 608')1220 ENDDO1221 ENDDO1222 #endif1223 #ifdef ISOVERIF1224 DO i=1,klon1225 IF (iso_eau >= 0) THEN1226 CALL iso_verif_egalite_choix(Rland_ice(iso_eau,i),1.0, &1227 & 'pbl_surf_mod 585',errmax,errmaxrel)1228 CALL iso_verif_egalite_choix(xtsnow_f(iso_eau,i),snow_f(i), &1229 & 'pbl_surf_mod 594',errmax,errmaxrel)1230 IF (iso_verif_egalite_choix_nostop(xtsol(iso_eau,i),qsol(i), &1231 & 'pbl_surf_mod 596',errmax,errmaxrel) == 1) THEN1232 WRITE(*,*) 'i=',i1233 STOP1234 ENDIF1235 DO nsrf=1,nbsrf1236 CALL iso_verif_egalite_choix(xtsnow(iso_eau,i,nsrf),snow(i,nsrf), &1237 & 'pbl_surf_mod 598',errmax,errmaxrel)1238 ENDDO1239 ENDIF !IF (iso_eau >= 0) THEN1240 ENDDO !DO i=1,knon1241 DO k=1,klev1242 DO i=1,klon1243 IF (iso_eau >= 0) THEN1244 CALL iso_verif_egalite_choix(xt(iso_eau,i,k),q(i,k), &1245 & 'pbl_surf_mod 595',errmax,errmaxrel)1246 ENDIF !IF (iso_eau >= 0) THEN1247 ENDDO !DO i=1,knon1248 ENDDO !DO k=1,klev1249 #endif1250 #endif1251 1252 1253 1254 !****************************************************************************************1255 ! Force soil water content to qsol0 if qsol0>0 and VEGET=F (use bucket1256 ! instead of ORCHIDEE)1257 IF (qsol0>=0.) THEN1258 PRINT*,'WARNING : On impose qsol=',qsol01259 qsol(:)=qsol01260 #ifdef ISO1261 DO ixt=1,niso1262 xtsol(ixt,:)=qsol0*Rdefault(ixt)1263 ENDDO1264 #ifdef ISOTRAC1265 DO ixt=1+niso,ntraciso1266 xtsol(ixt,:)=qsol0*Rdefault(index_iso(ixt))1267 ENDDO1268 #endif1269 #endif1270 ENDIF1271 !****************************************************************************************1272 1273 !****************************************************************************************1274 ! 2) Initialization to zero1275 !****************************************************************************************1276 !1277 ! 2a) Initialization of all argument variables with INTENT(OUT)1278 !****************************************************************************************1279 cdragh(:)=0. ; cdragm(:)=0.1280 zu1(:)=0. ; zv1(:)=0.1281 yus0(:)=0. ; yvs0(:)=0.1282 !albedo SB >>>1283 alb_dir_m=0. ; alb_dif_m=0. ; alb3_lic(:)=0.1284 !albedo SB <<<1285 zxsens(:)=0. ; zxevap(:)=0. ; zxtsol(:)=0. ; zxsnowerosion(:)=0.1286 d_t_w(:,:)=0. ; d_q_w(:,:)=0. ; d_t_x(:,:)=0. ; d_q_x(:,:)=0.1287 zxfluxlat(:)=0.1288 zt2m(:)=0. ; zq2m(:)=0. ; qsat2m(:)=0. ; rh2m(:)=0.1289 zn2mout(:,:)=0 ;1290 d_t(:,:)=0. ; d_t_diss(:,:)=0. ; d_q(:,:)=0. ; d_qbs(:,:)=0. ; d_u(:,:)=0. ; d_v(:,:)=0.1291 zcoefh(:,:,:)=0. ; zcoefm(:,:,:)=0.1292 zxsens_x(:)=0. ; zxsens_w(:)=0. ; zxfluxlat_x(:)=0. ; zxfluxlat_w(:)=0.1293 cdragh_x(:)=0. ; cdragh_w(:)=0. ; cdragm_x(:)=0. ; cdragm_w(:)=0.1294 kh(:)=0. ; kh_x(:)=0. ; kh_w(:)=0.1295 slab_wfbils(:)=0.1296 s_pblh(:)=0. ; s_pblh_x(:)=0. ; s_pblh_w(:)=0.1297 s_plcl(:)=0. ; s_plcl_x(:)=0. ; s_plcl_w(:)=0.1298 s_capCL(:)=0. ; s_oliqCL(:)=0. ; s_cteiCL(:)=0. ; s_pblT(:)=0.1299 s_therm(:)=0.1300 s_trmb1(:)=0. ; s_trmb2(:)=0. ; s_trmb3(:)=0.1301 zustar(:)=0.1302 zu10m(:)=0. ; zv10m(:)=0.1303 fder_print(:)=0.1304 zxqsurf(:)=0.1305 delta_qsurf(:) = 0.1306 zxfluxu(:,:)=0. ; zxfluxv(:,:)=0.1307 solsw(:,:)=0. ; sollw(:,:)=0.1308 d_ts(:,:)=0.1309 evap(:,:)=0.1310 snowerosion(:,:)=0.1311 fluxlat(:,:)=0.1312 wfbils(:,:)=0. ; wfevap(:,:)=0. ;1313 flux_t(:,:,:)=0. ; flux_q(:,:,:)=0. ; flux_u(:,:,:)=0. ; flux_v(:,:,:)=0.1314 flux_qbs(:,:,:)=0.1315 dflux_t(:)=0. ; dflux_q(:)=0.1316 zxsnow(:)=0.1317 zxfluxt(:,:)=0. ; zxfluxq(:,:)=0.; zxfluxqbs(:,:)=0.1318 qsnow(:)=0. ; snowhgt(:)=0. ; to_ice(:)=0. ; sissnow(:)=0.1319 runoff(:)=0. ; icesub_lic(:)=0.1320 #ifdef ISO1321 zxxtevap(:,:)=0.1322 d_xt(:,:,:)=0.1323 d_xt_x(:,:,:)=0.1324 d_xt_w(:,:,:)=0.1325 flux_xt(:,:,:,:)=0.1326 ! xtsnow(:,:,:)=0.! attention, xtsnow est l'équivalent de snow et non de qsnow1327 xtevap(:,:,:)=0.1328 #endif1329 IF (iflag_pbl<20.or.iflag_pbl>=30) THEN1330 zcoefh(:,:,:) = 0.01331 zcoefh(:,1,:) = 999999. ! zcoefh(:,k=1) should never be used1332 zcoefm(:,:,:) = 0.01333 zcoefm(:,1,:) = 999999. !1334 ELSE1335 zcoefm(:,:,is_ave)=0.1336 zcoefh(:,:,is_ave)=0.1337 ENDIF1338 !!1339 ! The components "is_ave" of tke_x and wake_deltke are "OUT" variables1340 !jyg<1341 !! tke(:,:,is_ave)=0.1342 tke_x(:,:,is_ave)=0.1343 eps_x(:,:,is_ave)=0.1344 1345 wake_dltke(:,:,is_ave)=0.1346 !>jyg1347 !!! jyg le 23/02/20131348 t2m(:,:) = 999999. ! t2m and q2m are meaningfull only over sub-surfaces1349 q2m(:,:) = 999999. ! actually present in the grid cell.1350 !!!1351 rh2m(:) = 0. ; qsat2m(:) = 0.1352 !!!1353 !!! jyg le 10/02/20121354 rh2m_x(:) = 0. ; qsat2m_x(:) = 0. ; rh2m_w(:) = 0. ; qsat2m_w(:) = 0.1355 1356 ! 2b) Initialization of all local variables that will be compressed later1357 !****************************************************************************************1358 !! cdragh = 0.0 ; cdragm = 0.0 ; dflux_t = 0.0 ; dflux_q = 0.01359 ypct = 0.0 ; yts = 0.0 ; ysnow = 0.01360 !! zv1 = 0.0 ; yqsurf = 0.01361 !albedo SB >>>1362 yqsurf = 0.0 ; yalb = 0.0 ; yalb_vis = 0.01363 !albedo SB <<<1364 yrain_f = 0.0 ; ysnow_f = 0.0 ; ybs_f=0.0 ; yfder = 0.0 ; ysolsw = 0.01365 ysollw = 0.0 ; yz0m = 0.0 ; yz0h = 0.0 ; yz0h_oupas = 0.0 ; yu1 = 0.01366 yv1 = 0.0 ; ypaprs = 0.0 ; ypplay = 0.0 ; yqbs1 = 0.01367 ydelp = 0.0 ; yu = 0.0 ; yv = 0.0 ; yt = 0.01368 yq = 0.0 ; y_dflux_t = 0.0 ; y_dflux_q = 0.01369 yqbs(:,:)=0.01370 yrugoro = 0.0 ; ywindsp = 0.01371 !! d_ts = 0.0 ; yfluxlat=0.0 ; flux_t = 0.0 ; flux_q = 0.01372 yfluxlat=0.0 ; y_flux0(:)=0.01373 !! flux_u = 0.0 ; flux_v = 0.0 ; d_t = 0.0 ; d_q = 0.01374 !! d_t_diss= 0.0 ;d_u = 0.0 ; d_v = 0.01375 yqsol = 0.01376 1377 ytke=0.1378 yeps=0.1379 yri0(:)=0.1380 !FC1381 y_treedrg=0.1382 1383 ! Martin1384 ysnowhgt = 0.0; yqsnow = 0.0 ; yrunoff = 0.0 ; ytoice =0.01385 yalb3_new = 0.0 ; ysissnow = 0.01386 ycldt = 0.0 ; yrmu0 = 0.01387 ! Martin1388 y_d_qbs(:,:)=0.01389 1390 !!! nrlmd+jyg le 02/05/2011 et le 20/02/20121391 ytke_x=0. ; ytke_w=0. ; ywake_dltke=0.1392 yeps_x=0. ; yeps_w=0.1393 y_d_t_x=0. ; y_d_t_w=0. ; y_d_q_x=0. ; y_d_q_w=0.1394 !! d_t_w=0. ; d_q_w=0.1395 !! d_t_x=0. ; d_q_x=0.1396 !! d_wake_dlt=0. ; d_wake_dlq=0.1397 yfluxlat_x=0. ; yfluxlat_w=0.1398 ywake_s=0. ; ywake_cstar=0. ;ywake_dens=0.1399 !!!1400 !!! nrlmd le 13/06/20111401 tau_eq=0. ; delta_coef=0.1402 y_delta_flux_t1=0.1403 ydtsurf_th=0.1404 yts_x(:)=0. ; yts_w(:)=0.1405 y_delta_tsurf(:)=0. ; y_delta_qsurf(:)=0.1406 yqsurf_x(:)=0. ; yqsurf_w(:)=0.1407 yg_T(:) = 0. ; yg_Q(:) = 0.1408 yGamma_dTs_phiT(:) = 0. ; yGamma_dQs_phiQ(:) = 0.1409 ydTs_ins(:) = 0. ; ydqs_ins(:) = 0.1410 1411 !!!1412 ytsoil = 999999.1413 !FC1414 y_d_u_frein(:,:)=0.1415 y_d_v_frein(:,:)=0.1416 !FC1417 1418 #ifdef ISO1419 yxtrain_f = 0.0 ; yxtsnow_f = 0.01420 yxtsnow = 0.01421 yxt = 0.01422 yxtsol = 0.01423 flux_xt = 0.01424 yRland_ice = 0.01425 ! d_xt = 0.01426 y_dflux_xt = 0.01427 dflux_xt=0.01428 y_d_xt_x=0. ; y_d_xt_w=0.1429 #endif1430 1431 ! >> PC1432 !the yfields_out variable is defined in (klon,nbcf_out) even if it is used on1433 !the ORCHIDEE grid and as such should be defined in yfields_out(knon,nbcf_out) but1434 !the knon variable is not known at that level of pbl_surface_mod1435 1436 !the yfields_in variable is defined in (klon,nbcf_in) even if it is used on the1437 !ORCHIDEE grid and as such should be defined in yfields_in(knon,nbcf_in) but the1438 !knon variable is not known at that level of pbl_surface_mod1439 1440 yfields_out(:,:) = 0.1441 ! << PC1442 1443 !GG1444 ypphi = 0.01445 !GG1446 1447 1448 ! 2c) Initialization of all local variables computed within the subsurface loop and used later on1449 !****************************************************************************************1450 d_t_diss_x(:,:) = 0. ; d_t_diss_w(:,:) = 0.1451 d_u_x(:,:)=0. ; d_u_w(:,:)=0.1452 d_v_x(:,:)=0. ; d_v_w(:,:)=0.1453 flux_t_x(:,:,:)=0. ; flux_t_w(:,:,:)=0.1454 flux_q_x(:,:,:)=0. ; flux_q_w(:,:,:)=0.1455 !1456 !jyg<1457 flux_u_x(:,:,:)=0. ; flux_u_w(:,:,:)=0.1458 flux_v_x(:,:,:)=0. ; flux_v_w(:,:,:)=0.1459 fluxlat_x(:,:)=0. ; fluxlat_w(:,:)=0.1460 !>jyg1461 #ifdef ISO1462 flux_xt_x(:,:,:,:)=0. ; flux_xt_w(:,:,:,:)=0.1463 #endif1464 !1465 !jyg<1466 ! pblh,plcl,capCL,cteiCL ... are meaningfull only over sub-surfaces1467 ! actually present in the grid cell ==> value set to 999999.1468 !1469 !jyg<1470 ustar(:,:) = 999999.1471 wstar(:,:) = 999999.1472 windsp(:,:) = SQRT(u10m(:,:)**2 + v10m(:,:)**2 )1473 u10m(:,:) = 999999.1474 v10m(:,:) = 999999.1475 !>jyg1476 !1477 pblh(:,:) = 999999. ! Hauteur de couche limite1478 plcl(:,:) = 999999. ! Niveau de condensation de la CLA1479 capCL(:,:) = 999999. ! CAPE de couche limite1480 oliqCL(:,:) = 999999. ! eau_liqu integree de couche limite1481 cteiCL(:,:) = 999999. ! cloud top instab. crit. couche limite1482 pblt(:,:) = 999999. ! T a la Hauteur de couche limite1483 therm(:,:) = 999999.1484 trmb1(:,:) = 999999. ! deep_cape1485 trmb2(:,:) = 999999. ! inhibition1486 trmb3(:,:) = 999999. ! Point Omega1487 !1488 t2m_x(:,:) = 999999.1489 q2m_x(:,:) = 999999.1490 ustar_x(:,:) = 999999.1491 wstar_x(:,:) = 999999.1492 u10m_x(:,:) = 999999.1493 v10m_x(:,:) = 999999.1494 !1495 pblh_x(:,:) = 999999. ! Hauteur de couche limite1496 plcl_x(:,:) = 999999. ! Niveau de condensation de la CLA1497 capCL_x(:,:) = 999999. ! CAPE de couche limite1498 oliqCL_x(:,:) = 999999. ! eau_liqu integree de couche limite1499 cteiCL_x(:,:) = 999999. ! cloud top instab. crit. couche limite1500 pblt_x(:,:) = 999999. ! T a la Hauteur de couche limite1501 therm_x(:,:) = 999999.1502 trmb1_x(:,:) = 999999. ! deep_cape1503 trmb2_x(:,:) = 999999. ! inhibition1504 trmb3_x(:,:) = 999999. ! Point Omega1505 !1506 t2m_w(:,:) = 999999.1507 q2m_w(:,:) = 999999.1508 ustar_w(:,:) = 999999.1509 wstar_w(:,:) = 999999.1510 u10m_w(:,:) = 999999.1511 v10m_w(:,:) = 999999.1512 1513 pblh_w(:,:) = 999999. ! Hauteur de couche limite1514 plcl_w(:,:) = 999999. ! Niveau de condensation de la CLA1515 capCL_w(:,:) = 999999. ! CAPE de couche limite1516 oliqCL_w(:,:) = 999999. ! eau_liqu integree de couche limite1517 cteiCL_w(:,:) = 999999. ! cloud top instab. crit. couche limite1518 pblt_w(:,:) = 999999. ! T a la Hauteur de couche limite1519 therm_w(:,:) = 999999.1520 trmb1_w(:,:) = 999999. ! deep_cape1521 trmb2_w(:,:) = 999999. ! inhibition1522 trmb3_w(:,:) = 999999. ! Point Omega1523 !!!1524 !1525 !!!1526 !****************************************************************************************1527 ! 3) - Calculate pressure thickness of each layer1528 ! - Calculate the wind at first layer1529 ! - Mean calculations of albedo1530 ! - Calculate net radiance at sub-surface1531 !****************************************************************************************1532 DO k = 1, klev1533 DO i = 1, klon1534 delp(i,k) = paprs(i,k)-paprs(i,k+1)1535 ENDDO1536 ENDDO1537 1538 !****************************************************************************************1539 ! Test for rugos........ from physiq.. A la fin plutot???1540 !1541 !****************************************************************************************1542 1543 DO nsrf = 1, nbsrf1544 DO i = 1, klon1545 z0m(i,nsrf) = MAX(z0m(i,nsrf),z0min)1546 z0h(i,nsrf) = MAX(z0h(i,nsrf),z0min)1547 ENDDO1548 ENDDO1549 1550 ! AM heterogeneous continental subsurfaces1551 ! compute time-independent effective surface parameters1552 IF (iflag_hetero_surf .GT. 0) THEN1553 CALL eff_surf_param(klon, nbtersrf, albedo_tersrf, frac_tersrf, 'ARI', albedo_eff)1554 ENDIF1555 1556 ! Mean calculations of albedo1557 !1558 ! * alb : mean albedo for whole SW interval1559 !1560 ! Mean albedo for grid point1561 ! * alb_m : mean albedo at whole SW interval1562 1563 alb_dir_m(:,:) = 0.01564 alb_dif_m(:,:) = 0.01565 DO k = 1, nsw1566 DO nsrf = 1, nbsrf1567 DO i = 1, klon1568 ! AM heterogeneous continental sub-surfaces1569 IF (nsrf .EQ. is_ter .AND. iflag_hetero_surf .GT. 0) THEN1570 alb_dir(i,k,nsrf) = albedo_eff(i)1571 alb_dif(i,k,nsrf) = albedo_eff(i)1572 ENDIF1573 !1574 alb_dir_m(i,k) = alb_dir_m(i,k) + alb_dir(i,k,nsrf) * pctsrf(i,nsrf)1575 alb_dif_m(i,k) = alb_dif_m(i,k) + alb_dif(i,k,nsrf) * pctsrf(i,nsrf)1576 ENDDO1577 ENDDO1578 ENDDO1579 1580 ! We here suppose the fraction f1 of incoming radiance of visible radiance1581 ! as a fraction of all shortwave radiance1582 f1 = 0.51583 ! f1 = 1 ! put f1=1 to recreate old calculations1584 1585 !f1 is already included with SFRWL values in each surf files1586 alb=0.01587 DO k=1,nsw1588 DO nsrf = 1, nbsrf1589 DO i = 1, klon1590 alb(i,nsrf) = alb(i,nsrf) + alb_dir(i,k,nsrf)*SFRWL(k)1591 ENDDO1592 ENDDO1593 ENDDO1594 1595 alb_m=0.01596 DO k = 1,nsw1597 DO i = 1, klon1598 alb_m(i) = alb_m(i) + alb_dir_m(i,k)*SFRWL(k)1599 ENDDO1600 ENDDO1601 !albedo SB <<<1602 1603 1604 1605 ! Calculation of mean temperature at surface grid points1606 ztsol(:) = 0.01607 DO nsrf = 1, nbsrf1608 DO i = 1, klon1609 ztsol(i) = ztsol(i) + ts(i,nsrf)*pctsrf(i,nsrf)1610 ENDDO1611 ENDDO1612 1613 ! Linear distrubution on sub-surface of long- and shortwave net radiance1614 DO nsrf = 1, nbsrf1615 DO i = 1, klon1616 sollw(i,nsrf) = sollw_m(i) + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ts(i,nsrf))1617 !--OB this line is not satisfactory because alb is the direct albedo not total albedo1618 solsw(i,nsrf) = solsw_m(i) * (1.-alb(i,nsrf)) / (1.-alb_m(i))1619 ENDDO1620 ENDDO1621 !1622 !<al1: second order corrections1623 !- net = dwn -up; up=sig( T4 + 4sum%T3T' + 6sum%T2T'2 +...)1624 IF (iflag_order2_sollw == 1) THEN1625 meansqT(:) = 0. ! as working buffer1626 DO nsrf = 1, nbsrf1627 DO i = 1, klon1628 meansqT(i) = meansqT(i)+(ts(i,nsrf)-ztsol(i))**2 *pctsrf(i,nsrf)1629 ENDDO1630 ENDDO1631 DO nsrf = 1, nbsrf1632 DO i = 1, klon1633 sollw(i,nsrf) = sollw(i,nsrf) &1634 + 6.0*RSIGMA*ztsol(i)**2 *(meansqT(i)-(ztsol(i)-ts(i,nsrf))**2)1635 ENDDO1636 ENDDO1637 ENDIF ! iflag_order2_sollw == 11638 !>al11639 1640 !--OB add diffuse fraction of SW down1641 DO n=1,nbcf_out1642 IF (cfname_out(n) == "swdownfdiff" ) fields_out(:,n) = solswfdiff_m(:)1643 ENDDO1644 ! >> PC1645 IF (carbon_cycle_cpl .AND. carbon_cycle_tr .AND. nbcf_out.GT.0 ) THEN1646 r_co2_ppm(:) = co2_send(:)1647 DO n=1,nbcf_out1648 IF (cfname_out(n) == "atmco2" ) fields_out(:,n) = co2_send(:)1649 ENDDO1650 ENDIF1651 IF ( .NOT. carbon_cycle_tr .AND. nbcf_out.GT.0 ) THEN1652 r_co2_ppm(:) = co2_ppm ! Constant field1653 DO n=1,nbcf_out1654 IF (cfname_out(n) == "atmco2" ) fields_out(:,n) = co2_ppm1655 ENDDO1656 ENDIF1657 ! << PC1658 1659 !****************************************************************************************1660 ! 4) Loop over different surfaces1661 !1662 ! Only points containing a fraction of the sub surface will be treated.1663 !1664 !****************************************************************************************1665 !<<<<<<<<<<<<<1666 loop_nbsrf: DO nsrf = 1, nbsrf !<<<<<<<<<<<<<1667 !<<<<<<<<<<<<<1668 IF (prt_level >=10) print *,' Loop nsrf ',nsrf1669 !1670 IF (iflag_split_ref == 3) THEN1671 IF (nsrf == is_oce) THEN1672 iflag_split = 11673 ELSE1674 iflag_split=01675 ENDIF !! (nsrf == is_oce)1676 ELSE1677 iflag_split = iflag_split_ref1678 ENDIF !! (iflag_split_ref == 3)1679 1680 ! Search for index(ni) and size(knon) of domaine to treat1681 ni(:) = 01682 knon = 01683 DO i = 1, klon1684 IF (pctsrf(i,nsrf) > 0.) THEN1685 knon = knon + 11686 ni(knon) = i1687 ENDIF1688 ENDDO1689 1690 !!! jyg le 19/08/20121691 ! IF (knon <= 0) THEN1692 ! IF (prt_level >= 10) print *,' no grid point for nsrf= ',nsrf1693 ! cycle loop_nbsrf1694 ! ENDIF1695 !!!1696 1697 1698 !****************************************************************************************1699 ! 5) Compress variables1700 !1701 !****************************************************************************************1702 1703 !1704 !jyg< (20190926)1705 ! Provisional : set ybeta to standard values1706 IF (nsrf .NE. is_ter) THEN1707 ybeta(1:knon) = 1.1708 ELSE1709 IF (iflag_split .EQ. 0) THEN1710 ybeta(1:knon) = 1.1711 ELSE1712 DO j = 1, knon1713 i = ni(j)1714 ybeta(j) = beta(i,nsrf)1715 ENDDO1716 ENDIF ! (iflag_split .LE.1)1717 ENDIF ! (nsrf .NE. is_ter)1718 !>jyg1719 !1720 DO j = 1, knon1721 i = ni(j)1722 ypct(j) = pctsrf(i,nsrf)1723 yts(j) = ts(i,nsrf)1724 ysnow(j) = snow(i,nsrf)1725 yqsurf(j) = qsurf(i,nsrf)1726 yalb(j) = alb(i,nsrf)1727 !albedo SB >>>1728 yalb_vis(j) = alb_dir(i,1,nsrf)1729 IF (nsw==6) THEN1730 yalb_vis(j)=(alb_dir(i,1,nsrf)*SFRWL(1)+alb_dir(i,2,nsrf)*SFRWL(2) &1731 +alb_dir(i,3,nsrf)*SFRWL(3))/(SFRWL(1)+SFRWL(2)+SFRWL(3))1732 ENDIF1733 !albedo SB <<<1734 yrain_f(j) = rain_f(i)1735 ysnow_f(j) = snow_f(i)1736 ybs_f(j) = bs_f(i)1737 yagesno(j) = agesno(i,nsrf)1738 yfder(j) = fder(i)1739 ylwdown(j) = lwdown_m(i)1740 ygustiness(j) = gustiness(i)1741 ysolsw(j) = solsw(i,nsrf)1742 ysollw(j) = sollw(i,nsrf)1743 yz0m(j) = z0m(i,nsrf)1744 yz0h(j) = z0h(i,nsrf)1745 yrugoro(j) = rugoro(i)1746 yu1(j) = u(i,1)1747 yv1(j) = v(i,1)1748 yqbs1(j) = qbs(i,1)1749 ypaprs(j,klev+1) = paprs(i,klev+1)1750 !jyg<1751 !! ywindsp(j) = SQRT(u10m(i,nsrf)**2 + v10m(i,nsrf)**2 )1752 ywindsp(j) = windsp(i,nsrf)1753 !>jyg1754 ! Martin and Etienne1755 yzmea(j) = zmea(i)1756 yzsig(j) = zsig(i)1757 ycldt(j) = cldt(i)1758 yrmu0(j) = rmu0(i)1759 ! Martin1760 !!! nrlmd le 13/06/20111761 y_delta_tsurf(j)=delta_tsurf(i,nsrf)1762 yfluxbs(j)=0.01763 y_flux_bs(j) = 0.01764 !!!1765 #ifdef ISO1766 DO ixt=1,ntraciso1767 yxtrain_f(ixt,j) = xtrain_f(ixt,i)1768 yxtsnow_f(ixt,j) = xtsnow_f(ixt,i)1769 ENDDO1770 DO ixt=1,niso1771 yxtsnow(ixt,j) = xtsnow(ixt,i,nsrf)1772 ENDDO1773 !IF (nsrf == is_lic) THEN1774 DO ixt=1,niso1775 yRland_ice(ixt,j)= Rland_ice(ixt,i)1776 ENDDO1777 !endif !IF (nsrf == is_lic) THEN1778 #ifdef ISOVERIF1779 IF (iso_eau >= 0) THEN1780 call iso_verif_egalite_choix(ysnow_f(j), &1781 & yxtsnow_f(iso_eau,j),'pbl_surf_mod 862', &1782 & errmax,errmaxrel)1783 call iso_verif_egalite_choix(ysnow(j), &1784 & yxtsnow(iso_eau,j),'pbl_surf_mod 872', &1785 & errmax,errmaxrel)1786 ENDIF1787 #endif1788 #ifdef ISOVERIF1789 DO ixt=1,ntraciso1790 call iso_verif_noNaN(yxtsnow_f(ixt,j),'pbl_surf_mod 921')1791 ENDDO1792 #endif1793 #endif1794 ENDDO1795 ! >> PC1796 !--compressing fields_out onto ORCHIDEE grid1797 !--these fields are shared and used directly surf_land_orchidee_mod1798 DO n = 1, nbcf_out1799 DO j = 1, knon1800 i = ni(j)1801 yfields_out(j,n) = fields_out(i,n)1802 ENDDO1803 ENDDO1804 ! << PC1805 DO k = 1, klev1806 DO j = 1, knon1807 i = ni(j)1808 ypaprs(j,k) = paprs(i,k)1809 ypplay(j,k) = pplay(i,k)1810 ydelp(j,k) = delp(i,k)1811 ENDDO1812 ENDDO1813 !1814 !!! jyg le 07/02/2012 et le 10/04/20131815 DO k = 1, klev+11816 DO j = 1, knon1817 i = ni(j)1818 !jyg<1819 !! ytke(j,k) = tke(i,k,nsrf)1820 ytke(j,k) = tke_x(i,k,nsrf)1821 ENDDO1822 ENDDO1823 !>jyg1824 DO k = 1, klev1825 DO j = 1, knon1826 i = ni(j)1827 y_treedrg(j,k) = treedrg(i,k,nsrf)1828 yu(j,k) = u(i,k)1829 yv(j,k) = v(i,k)1830 yt(j,k) = t(i,k)1831 yq(j,k) = q(i,k)1832 yqbs(j,k)=qbs(i,k)1833 !! GG1834 ypphi(j,k) = pphi(i,k)1835 !!1836 1837 #ifdef ISO1838 DO ixt=1,ntraciso1839 yxt(ixt,j,k) = xt(ixt,i,k)1840 ENDDO !DO ixt=1,ntraciso1841 #endif1842 ENDDO1843 ENDDO1844 !1845 IF (iflag_split.GE.1) THEN1846 !!! nrlmd le 02/05/20111847 DO k = 1, klev1848 DO j = 1, knon1849 i = ni(j)1850 yu_x(j,k) = u(i,k)1851 yv_x(j,k) = v(i,k)1852 yt_x(j,k) = t(i,k)-wake_s(i)*wake_dlt(i,k)1853 yq_x(j,k) = q(i,k)-wake_s(i)*wake_dlq(i,k)1854 yu_w(j,k) = u(i,k)1855 yv_w(j,k) = v(i,k)1856 yt_w(j,k) = t(i,k)+(1.-wake_s(i))*wake_dlt(i,k)1857 yq_w(j,k) = q(i,k)+(1.-wake_s(i))*wake_dlq(i,k)1858 !!!1859 #ifdef ISO1860 DO ixt=1,ntraciso1861 yxt_x(ixt,j,k) = xt(ixt,i,k)-wake_s(i)*wake_dlxt(ixt,i,k)1862 yxt_w(ixt,j,k) = xt(ixt,i,k)+(1.-wake_s(i))*wake_dlxt(ixt,i,k)1863 ENDDO1864 #endif1865 ENDDO1866 ENDDO1867 1868 IF (prt_level .ge. 10) THEN1869 print *,'pbl_surface, wake_s(1), wake_dlt(1,:) ', wake_s(1), wake_dlt(1,:)1870 print *,'pbl_surface, wake_s(1), wake_dlq(1,:) ', wake_s(1), wake_dlq(1,:)1871 ENDIF1872 1873 !!! nrlmd le 02/05/20111874 DO k = 1, klev+11875 DO j = 1, knon1876 i = ni(j)1877 !jyg<1878 !! ytke_x(j,k) = tke(i,k,nsrf)-wake_s(i)*wake_dltke(i,k,nsrf)1879 !! ytke_w(j,k) = tke(i,k,nsrf)+(1.-wake_s(i))*wake_dltke(i,k,nsrf)1880 !! ywake_dltke(j,k) = wake_dltke(i,k,nsrf)1881 !! ytke(j,k) = tke(i,k,nsrf)1882 !1883 ytke_x(j,k) = tke_x(i,k,nsrf)1884 ytke(j,k) = tke_x(i,k,nsrf)+wake_s(i)*wake_dltke(i,k,nsrf)1885 ytke_w(j,k) = tke_x(i,k,nsrf)+wake_dltke(i,k,nsrf)1886 ywake_dltke(j,k) = wake_dltke(i,k,nsrf)1887 1888 !>jyg1889 ENDDO1890 ENDDO1891 !!!1892 !!! jyg le 07/02/20121893 DO j = 1, knon1894 i = ni(j)1895 ywake_s(j)=wake_s(i)1896 ywake_cstar(j)=wake_cstar(i)1897 ywake_dens(j)=wake_dens(i)1898 ENDDO1899 !!!1900 !!! nrlmd le 13/06/20111901 DO j=1,knon1902 yts_x(j)=yts(j)-ywake_s(j)*y_delta_tsurf(j)1903 yts_w(j)=yts(j)+(1.-ywake_s(j))*y_delta_tsurf(j)1904 ENDDO1905 !!!1906 ENDIF ! (iflag_split .ge.1)1907 !!!1908 DO k = 1, nsoilmx1909 DO j = 1, knon1910 i = ni(j)1911 ytsoil(j,k) = ftsoil(i,k,nsrf)1912 ENDDO1913 ENDDO1914 1915 ! qsol(water height in soil) only for bucket continental model1916 IF ( nsrf .EQ. is_ter .AND. .NOT. ok_veget ) THEN1917 DO j = 1, knon1918 i = ni(j)1919 yqsol(j) = qsol(i)1920 #ifdef ISO1921 DO ixt=1,niso1922 yxtsol(ixt,j) = xtsol(ixt,i)1923 ENDDO1924 #endif1925 ENDDO1926 ENDIF1927 1928 if (nsrf == is_oce .and. activate_ocean_skin >= 1) then1929 if (activate_ocean_skin == 2 .and. type_ocean == "couple") then1930 ydelta_sal(:knon) = delta_sal(ni(:knon))1931 ydelta_sst(:knon) = delta_sst(ni(:knon))1932 ydter(:knon) = dter(ni(:knon))1933 ydser(:knon) = dser(ni(:knon))1934 ydt_ds(:knon) = dt_ds(ni(:knon))1935 end if1936 1937 yds_ns(:knon) = ds_ns(ni(:knon))1938 ydt_ns(:knon) = dt_ns(ni(:knon))1939 end if1940 1941 !****************************************************************************************1942 ! 6a) Calculate coefficients for turbulent diffusion at surface, cdragh et cdragm.1943 !1944 !****************************************************************************************1945 1946 1947 !!! jyg le 07/02/20121948 IF (iflag_split .eq.0) THEN1949 !!!1950 !!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/20121951 ! Faire disparaitre les lignes commentees fin 2015 (le temps des tests)1952 ! CALL clcdrag( knon, nsrf, ypaprs, ypplay, &1953 ! yu(:,1), yv(:,1), yt(:,1), yq(:,1), &1954 ! yts, yqsurf, yrugos, &1955 ! ycdragm, ycdragh )1956 ! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag1957 DO i = 1, knon1958 ! print*,'PBL ',i,RD1959 ! print*,'PBL ',yt(i,1),ypaprs(i,1),ypplay(i,1)1960 zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &1961 * (ypaprs(i,1)-ypplay(i,1))1962 speed(i) = SQRT(yu(i,1)**2+yv(i,1)**2)1963 ENDDO1964 1965 !!! AM heterogeneous continental subsurfaces1966 IF (nsrf .EQ. is_ter) THEN1967 ! compute time-dependent effective surface parameters (function of zgeo1) !! AM1968 IF (iflag_hetero_surf .GT. 0) THEN1969 DO n=1,nbtersrf1970 DO j=1,knon1971 i = ni(j)1972 yfrac_tersrf(j,n) = frac_tersrf(i,n)1973 yz0m_tersrf(j,n) = z0m_tersrf(i,n)1974 IF (ratio_z0m_z0h_tersrf(i,n) .NE. 0.) THEN1975 yz0h_tersrf(j,n) = z0m_tersrf(i,n) / ratio_z0m_z0h_tersrf(i,n)1976 ELSE1977 yz0h_tersrf(j,n) = 0.1978 ENDIF1979 ENDDO1980 ENDDO1981 !1982 CALL eff_surf_param(knon, nbtersrf, yz0m_tersrf, yfrac_tersrf, 'CDN', yz0m, zgeo1/RG)1983 CALL eff_surf_param(knon, nbtersrf, yz0h_tersrf, yfrac_tersrf, 'CDN', yz0h, zgeo1/RG)1984 !1985 ENDIF1986 ENDIF1987 1988 !1989 CALL cdrag(knon, nsrf, &1990 speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1), s_pblh, &1991 yts, yqsurf, yz0m, yz0h, yri0, 0, &1992 ycdragm, ycdragh, zri1, pref, rain_f, zxtsol, ypplay(:,1))1993 1994 ! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 050820131995 IF (ok_prescr_ust) THEN1996 DO i = 1, knon1997 print *,'ycdragm avant=',ycdragm(i)1998 vent= sqrt(yu(i,1)*yu(i,1)+yv(i,1)*yv(i,1))1999 ! ycdragm(i) = ust*ust/(1.+(yu(i,1)*yu(i,1)+yv(i,1)*yv(i,1)))2000 ! ycdragm(i) = ust*ust/((1.+sqrt(yu(i,1)*yu(i,1)+yv(i,1)*yv(i,1))) &2001 ! *sqrt(yu(i,1)*yu(i,1)+yv(i,1)*yv(i,1)))2002 ycdragm(i) = ust*ust/(1.+vent)/vent2003 ! print *,'ycdragm ust yu yv apres=',ycdragm(i),ust,yu(i,1),yv(i,1)2004 ENDDO2005 ENDIF2006 2007 IF (prt_level >=10) print *,'cdrag -> ycdragh ', ycdragh(1:knon)2008 ELSE !(iflag_split .eq.0)2009 2010 ! Faire disparaitre les lignes commentees fin 2015 (le temps des tests)2011 ! CALL clcdrag( knon, nsrf, ypaprs, ypplay, &2012 ! yu_x(:,1), yv_x(:,1), yt_x(:,1), yq_x(:,1), &2013 ! yts_x, yqsurf, yrugos, &2014 ! ycdragm_x, ycdragh_x )2015 ! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag2016 DO i = 1, knon2017 zgeo1_x(i) = RD * yt_x(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &2018 * (ypaprs(i,1)-ypplay(i,1))2019 speed_x(i) = SQRT(yu_x(i,1)**2+yv_x(i,1)**2)2020 ENDDO2021 2022 2023 CALL cdrag(knon, nsrf, &2024 speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),s_pblh_x,&2025 yts_x, yqsurf_x, yz0m, yz0h, yri0, 0, &2026 ycdragm_x, ycdragh_x, zri1_x, pref_x, rain_f, zxtsol, ypplay(:,1) )2027 2028 ! --- special Dice. JYG+MPL 251120132029 IF (ok_prescr_ust) THEN2030 DO i = 1, knon2031 ! print *,'ycdragm_x avant=',ycdragm_x(i)2032 vent= sqrt(yu_x(i,1)*yu_x(i,1)+yv_x(i,1)*yv_x(i,1))2033 ycdragm_x(i) = ust*ust/(1.+vent)/vent2034 ! print *,'ycdragm_x ust yu yv apres=',ycdragm_x(i),ust,yu_x(i,1),yv_x(i,1)2035 ENDDO2036 ENDIF2037 IF (prt_level >=10) print *,'clcdrag -> ycdragh_x ', ycdragh_x(1:knon)2038 !2039 ! Faire disparaitre les lignes commentees fin 2015 (le temps des tests)2040 ! CALL clcdrag( knon, nsrf, ypaprs, ypplay, &2041 ! yu_w(:,1), yv_w(:,1), yt_w(:,1), yq_w(:,1), &2042 ! yts_w, yqsurf, yz0m, &2043 ! ycdragm_w, ycdragh_w )2044 ! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag2045 DO i = 1, knon2046 zgeo1_w(i) = RD * yt_w(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &2047 * (ypaprs(i,1)-ypplay(i,1))2048 speed_w(i) = SQRT(yu_w(i,1)**2+yv_w(i,1)**2)2049 ENDDO2050 CALL cdrag(knon, nsrf, &2051 speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),s_pblh_w,&2052 yts_w, yqsurf_w, yz0m, yz0h, yri0, 0, &2053 ycdragm_w, ycdragh_w, zri1_w, pref_w, rain_f, zxtsol, ypplay(:,1) )2054 !2055 IF(ok_bug_zg_wk_pbl) THEN2056 zgeo1(1:knon) = wake_s(1:knon)*zgeo1_w(1:knon) + (1.-wake_s(1:knon))*zgeo1_x(1:knon)2057 ELSE2058 zgeo1(1:knon) = ywake_s(1:knon)*zgeo1_w(1:knon) + (1.-ywake_s(1:knon))*zgeo1_x(1:knon)2059 ENDIF2060 2061 ! --- special Dice. JYG+MPL 25112013 puis BOMEX2062 IF (ok_prescr_ust) THEN2063 DO i = 1, knon2064 ! print *,'ycdragm_w avant=',ycdragm_w(i)2065 vent= sqrt(yu_w(i,1)*yu_w(i,1)+yv_w(i,1)*yv_w(i,1))2066 ycdragm_w(i) = ust*ust/(1.+vent)/vent2067 ! print *,'ycdragm_w ust yu yv apres=',ycdragm_w(i),ust,yu_w(i,1),yv_w(i,1)2068 ENDDO2069 ENDIF2070 IF (prt_level >=10) print *,'clcdrag -> ycdragh_w ', ycdragh_w(1:knon)2071 !!!2072 ENDIF ! (iflag_split .eq.0)2073 !!!2074 2075 2076 !****************************************************************************************2077 ! 6b) Calculate coefficients for turbulent diffusion in the atmosphere, ycoefh et ycoefm.2078 !2079 !****************************************************************************************2080 2081 !!! jyg le 07/02/20122082 IF (iflag_split .eq.0) THEN2083 !!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/20122084 IF (prt_level >=10) THEN2085 print *,' args coef_diff_turb: yu ', yu(1:knon,:)2086 print *,' args coef_diff_turb: yv ', yv(1:knon,:)2087 print *,' args coef_diff_turb: yq ', yq(1:knon,:)2088 print *,' args coef_diff_turb: yt ', yt(1:knon,:)2089 print *,' args coef_diff_turb: yts ', yts(1:knon)2090 print *,' args coef_diff_turb: yz0m ', yz0m(1:knon)2091 print *,' args coef_diff_turb: yqsurf ', yqsurf(1:knon)2092 print *,' args coef_diff_turb: ycdragm ', ycdragm(1:knon)2093 print *,' args coef_diff_turb: ycdragh ', ycdragh(1:knon)2094 print *,' args coef_diff_turb: ytke ', ytke(1:knon,:)2095 ENDIF2096 2097 IF (iflag_pbl>=50) THEN2098 CALL call_atke(dtime,knon,klev,nsrf,ni,ycdragm(1:knon), ycdragh(1:knon),yus0(1:knon),yvs0(1:knon),yts(1:knon), &2099 yu(1:knon,:),yv(1:knon,:),yt(1:knon,:),yq(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:), &2100 ytke(1:knon,:),yeps(1:knon,:), ycoefm(1:knon,:), ycoefh(1:knon,:))2101 2102 ELSE2103 2104 CALL coef_diff_turb(dtime, nsrf, knon, ni, &2105 ypaprs, ypplay, yu, yv, yq, yt, yts, yqsurf, ycdragm, &2106 ycoefm, ycoefh, ytke, yeps, y_treedrg)2107 ! ycoefm, ycoefh, ytke)2108 !FC y_treedrg ajoute2109 IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN2110 ! In this case, coef_diff_turb is called for the Cd only2111 DO k = 2, klev2112 DO j = 1, knon2113 i = ni(j)2114 ycoefh(j,k) = zcoefh(i,k,nsrf)2115 ycoefm(j,k) = zcoefm(i,k,nsrf)2116 ENDDO2117 ENDDO2118 ENDIF2119 2120 ENDIF ! iflag_pbl >= 502121 2122 IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh ',ycoefh(1:knon,:)2123 2124 2125 ELSE !(iflag_split .eq.0)2126 2127 2128 IF (prt_level >=10) THEN2129 print *,' args coef_diff_turb: yu_x ', yu_x(1:knon,:)2130 print *,' args coef_diff_turb: yv_x ', yv_x(1:knon,:)2131 print *,' args coef_diff_turb: yq_x ', yq_x(1:knon,:)2132 print *,' args coef_diff_turb: yt_x ', yt_x(1:knon,:)2133 print *,' args coef_diff_turb: yts_x ', yts_x(1:knon)2134 print *,' args coef_diff_turb: yqsurf ', yqsurf(1:knon)2135 print *,' args coef_diff_turb: ycdragm_x ', ycdragm_x(1:knon)2136 print *,' args coef_diff_turb: ycdragh_x ', ycdragh_x(1:knon)2137 print *,' args coef_diff_turb: ytke_x ', ytke_x(1:knon,:)2138 ENDIF2139 2140 2141 IF (iflag_pbl>=50) THEN2142 2143 CALL call_atke(dtime,knon,klev,nsrf,ni,ycdragm_x(1:knon),ycdragh_x(1:knon),yus0(1:knon),yvs0(1:knon),yts_x(1:knon), &2144 yu_x(1:knon,:),yv_x(1:knon,:),yt_x(1:knon,:),yq_x(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:), &2145 ytke_x(1:knon,:),yeps_x(1:knon,:),ycoefm_x(1:knon,:), ycoefh_x(1:knon,:))2146 2147 ELSE2148 2149 CALL coef_diff_turb(dtime, nsrf, knon, ni, &2150 ypaprs, ypplay, yu_x, yv_x, yq_x, yt_x, yts_x, yqsurf_x, ycdragm_x, &2151 ycoefm_x, ycoefh_x, ytke_x,yeps_x,y_treedrg)2152 ! ycoefm_x, ycoefh_x, ytke_x)2153 !FC doit on le mettre ( on ne l utilise pas si il y a du spliting)2154 IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN2155 ! In this case, coef_diff_turb is called for the Cd only2156 DO k = 2, klev2157 DO j = 1, knon2158 i = ni(j)2159 ycoefh_x(j,k) = zcoefh(i,k,nsrf)2160 ycoefm_x(j,k) = zcoefm(i,k,nsrf)2161 ENDDO2162 ENDDO2163 ENDIF2164 2165 ENDIF ! iflag_pbl >= 502166 2167 IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_x ',ycoefh_x(1:knon,:)2168 !2169 IF (prt_level >=10) THEN2170 print *,' args coef_diff_turb: yu_w ', yu_w(1:knon,:)2171 print *,' args coef_diff_turb: yv_w ', yv_w(1:knon,:)2172 print *,' args coef_diff_turb: yq_w ', yq_w(1:knon,:)2173 print *,' args coef_diff_turb: yt_w ', yt_w(1:knon,:)2174 print *,' args coef_diff_turb: yts_w ', yts_w(1:knon)2175 print *,' args coef_diff_turb: yqsurf ', yqsurf(1:knon)2176 print *,' args coef_diff_turb: ycdragm_w ', ycdragm_w(1:knon)2177 print *,' args coef_diff_turb: ycdragh_w ', ycdragh_w(1:knon)2178 print *,' args coef_diff_turb: ytke_w ', ytke_w(1:knon,:)2179 ENDIF2180 2181 IF (iflag_pbl>=50) THEN2182 2183 CALL call_atke(dtime,knon,klev,nsrf,ni,ycdragm_w(1:knon),ycdragh_w(1:knon),yus0(1:knon),yvs0(1:knon),yts_w(1:knon), &2184 yu_w(1:knon,:),yv_w(1:knon,:),yt_w(1:knon,:),yq_w(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:), &2185 ytke_w(1:knon,:),yeps_w(1:knon,:),ycoefm_w(1:knon,:),ycoefh_w(1:knon,:))2186 2187 ELSE2188 2189 CALL coef_diff_turb(dtime, nsrf, knon, ni, &2190 ypaprs, ypplay, yu_w, yv_w, yq_w, yt_w, yts_w, yqsurf_w, ycdragm_w, &2191 ycoefm_w, ycoefh_w, ytke_w,yeps_w,y_treedrg)2192 ! ycoefm_w, ycoefh_w, ytke_w)2193 IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN2194 ! In this case, coef_diff_turb is called for the Cd only2195 DO k = 2, klev2196 DO j = 1, knon2197 i = ni(j)2198 ycoefh_w(j,k) = zcoefh(i,k,nsrf)2199 ycoefm_w(j,k) = zcoefm(i,k,nsrf)2200 ENDDO2201 ENDDO2202 ENDIF2203 2204 ENDIF ! iflag_pbl >= 502205 2206 2207 IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_w ',ycoefh_w(1:knon,:)2208 2209 !!!jyg le 10/04/20132210 !! En attendant de traiter le transport des traceurs dans les poches froides, formule2211 !! arbitraire pour ycoefh et ycoefm2212 DO k = 2,klev2213 DO j = 1,knon2214 ycoefh(j,k) = ycoefh_x(j,k) + ywake_s(j)*(ycoefh_w(j,k) - ycoefh_x(j,k))2215 ycoefm(j,k) = ycoefm_x(j,k) + ywake_s(j)*(ycoefm_w(j,k) - ycoefm_x(j,k))2216 ENDDO2217 ENDDO2218 2219 2220 ENDIF ! (iflag_split .eq.0)2221 2222 2223 !****************************************************************************************2224 !2225 ! 8) "La descente" - "The downhill"2226 !2227 ! climb_hq_down and climb_wind_down calculate the coefficients2228 ! Ccoef_X et Dcoef_X for X=[H, Q, U, V].2229 ! Only the coefficients at surface for H and Q are returned.2230 !2231 !****************************************************************************************2232 2233 ! - Calculate the coefficients Ccoef_H, Ccoef_Q, Dcoef_H and Dcoef_Q2234 !!! jyg le 07/02/20122235 IF (iflag_split .eq.0) THEN2236 !!!2237 !!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/20122238 CALL climb_hq_down(knon, ni, ycoefh, ypaprs, ypplay, &2239 ydelp, yt, yq, dtime, &2240 !!! jyg le 09/05/20112241 CcoefH, CcoefQ, DcoefH, DcoefQ, &2242 Kcoef_hq, gama_q, gama_h, &2243 !!!2244 AcoefH, AcoefQ, BcoefH, BcoefQ &2245 #ifdef ISO2246 & ,yxt, CcoefXT, DcoefXT, gama_xt, AcoefXT, BcoefXT &2247 #endif2248 & )2249 ELSE !(iflag_split .eq.0)2250 CALL climb_hq_down(knon, ni, ycoefh_x, ypaprs, ypplay, &2251 ydelp, yt_x, yq_x, dtime, &2252 !!! nrlmd le 02/05/20112253 CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x, &2254 Kcoef_hq_x, gama_q_x, gama_h_x, &2255 !!!2256 AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x &2257 #ifdef ISO2258 & ,yxt_x, CcoefXT_x, DcoefXT_x, gama_xt_x, AcoefXT_x, BcoefXT_x &2259 #endif2260 & )2261 !!!2262 IF (prt_level >=10) THEN2263 PRINT *,'pbl_surface (climb_hq_down.x->) AcoefH_x ',AcoefH_x2264 PRINT *,'pbl_surface (climb_hq_down.x->) AcoefQ_x ',AcoefQ_x2265 PRINT *,'pbl_surface (climb_hq_down.x->) BcoefH_x ',BcoefH_x2266 PRINT *,'pbl_surface (climb_hq_down.x->) BcoefQ_x ',BcoefQ_x2267 ENDIF2268 !2269 CALL climb_hq_down(knon, ni, ycoefh_w, ypaprs, ypplay, &2270 ydelp, yt_w, yq_w, dtime, &2271 !!! nrlmd le 02/05/20112272 CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w, &2273 Kcoef_hq_w, gama_q_w, gama_h_w, &2274 !!!2275 AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w &2276 #ifdef ISO2277 & ,yxt_w, CcoefXT_w, DcoefXT_w, gama_xt_w, AcoefXT_w, BcoefXT_w &2278 #endif2279 & )2280 !!!2281 IF (prt_level >=10) THEN2282 PRINT *,'pbl_surface (climb_hq_down.w->) AcoefH_w ',AcoefH_w2283 PRINT *,'pbl_surface (climb_hq_down.w->) AcoefQ_w ',AcoefQ_w2284 PRINT *,'pbl_surface (climb_hq_down.w->) BcoefH_w ',BcoefH_w2285 PRINT *,'pbl_surface (climb_hq_down.w->) BcoefQ_w ',BcoefQ_w2286 ENDIF2287 !!!2288 ENDIF ! (iflag_split .eq.0)2289 !!!2290 2291 ! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V2292 !!! jyg le 07/02/20122293 IF (iflag_split .eq.0) THEN2294 !!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/20122295 CALL climb_wind_down(knon, ni, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv, &2296 !!! jyg le 09/05/20112297 CcoefU, CcoefV, DcoefU, DcoefV, &2298 Kcoef_m, alf_1, alf_2, &2299 !!!2300 AcoefU, AcoefV, BcoefU, BcoefV)2301 ELSE ! (iflag_split .eq.0)2302 CALL climb_wind_down(knon, ni, dtime, ycoefm_x, ypplay, ypaprs, yt_x, ydelp, yu_x, yv_x, &2303 !!! nrlmd le 02/05/20112304 CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x, &2305 Kcoef_m_x, alf_1_x, alf_2_x, &2306 !!!2307 AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x)2308 !2309 CALL climb_wind_down(knon, ni, dtime, ycoefm_w, ypplay, ypaprs, yt_w, ydelp, yu_w, yv_w, &2310 !!! nrlmd le 02/05/20112311 CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w, &2312 Kcoef_m_w, alf_1_w, alf_2_w, &2313 !!!2314 AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w)2315 !!!2316 ENDIF ! (iflag_split .eq.0)2317 !!!2318 2319 ! For blowing snow:2320 IF (ok_bs) THEN2321 ! following Bintanja et al 2000, part II and Vionnet V PhD thesis2322 ! we assume that the eddy diffsivity coefficient for2323 ! suspended particles is a fraction of Kh2324 do k=1,klev2325 do j=1,knon2326 ycoefqbs(j,k)=ycoefh(j,k)*zeta_bs2327 enddo2328 enddo2329 CALL climb_qbs_down(knon, ni, ycoefqbs, ypaprs, ypplay, &2330 ydelp, yt, yqbs, dtime, &2331 CcoefQBS, DcoefQBS, &2332 Kcoef_qbs, gama_qbs, &2333 AcoefQBS, BcoefQBS)2334 ENDIF2335 2336 !****************************************************************************************2337 ! 9) Small calculations2338 !2339 !****************************************************************************************2340 2341 ! - Reference pressure is given the values at surface level2342 ypsref(:) = ypaprs(:,1)2343 2344 ! - CO2 field on 2D grid to be sent to ORCHIDEE2345 ! Transform to compressed field2346 IF (carbon_cycle_cpl) THEN2347 DO i=1,knon2348 r_co2_ppm(i) = co2_send(ni(i))2349 ENDDO2350 ELSE2351 r_co2_ppm(:) = co2_ppm ! Constant field2352 ENDIF2353 2354 !!! nrlmd le 02/05/2011 -----------------------On raccorde les 2 colonnes dans la couche 12355 !----------------------------------------------------------------------------------------2356 !!! jyg le 07/02/20122357 !!! jyg le 01/02/20172358 IF (iflag_split .eq. 0) THEN2359 yt1(:) = yt(:,1)2360 yq1(:) = yq(:,1)2361 #ifdef ISO2362 yxt1(:,:) = yxt(:,:,1)2363 #endif2364 2365 ELSE IF (iflag_split .ge. 1) THEN2366 #ifdef ISO2367 call abort_physic('pbl_surface_mod 2149','isos pas encore dans iflag_split=1',1)2368 #endif2369 2370 !2371 ! Cdragq computation2372 ! ------------------2373 !******************************************************************************2374 ! Cdragq computed from cdrag2375 ! The difference comes only from a factor (f_z0qh_oce) on z0, so that2376 ! it can be computed inside wx_pbl0_merge2377 ! More complicated appraches may require the propagation through2378 ! pbl_surface of an independant cdragq variable.2379 !******************************************************************************2380 !2381 IF ( f_z0qh_oce .ne. 1. .and. nsrf .eq.is_oce) THEN2382 ! Si on suit les formulations par exemple de Tessel, on2383 ! a z0h=0.4*nu/u*, z0q=0.62*nu/u*, d'ou f_z0qh_oce=0.62/0.4=1.552384 !! ycdragq_x(1:knon)=ycdragh_x(1:knon)* &2385 !! log(z1lay(1:knon)/yz0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*yz0h(1:knon)))2386 !! ycdragq_w(1:knon)=ycdragh_w(1:knon)* &2387 !! log(z1lay(1:knon)/yz0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*yz0h(1:knon)))2388 !2389 DO j = 1,knon2390 z1lay = zgeo1(j)/RG2391 fact_cdrag = log(z1lay/yz0h(j))/log(z1lay/(f_z0qh_oce*yz0h(j)))2392 ycdragq_x(j)=ycdragh_x(j)*fact_cdrag2393 ycdragq_w(j)=ycdragh_w(j)*fact_cdrag2394 !! Print *,'YYYYpbl0: fact_cdrag ', fact_cdrag2395 ENDDO ! j = 1,knon2396 !2397 !! Print *,'YYYYpbl0: z1lay, yz0h, f_z0qh_oce, ycdragh_w, ycdragq_w ', &2398 !! z1lay, yz0h(1:knon), f_z0qh_oce, ycdragh_w(1:knon), ycdragq_w(1:knon)2399 ELSE2400 ycdragq_x(1:knon)=ycdragh_x(1:knon)2401 ycdragq_w(1:knon)=ycdragh_w(1:knon)2402 ENDIF ! ( f_z0qh_oce .ne. 1. .and. nsrf .eq.is_oce)2403 !2404 CALL wx_pbl_prelim_0(knon, nsrf, dtime, ypplay, ypaprs, ywake_s, &2405 yts, y_delta_tsurf, ygustiness, &2406 yt_x, yt_w, yq_x, yq_w, &2407 yu_x, yu_w, yv_x, yv_w, &2408 ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &2409 ycdragm_x, ycdragm_w, &2410 AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &2411 AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &2412 BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &2413 BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &2414 Kech_h_x, Kech_h_w, Kech_h &2415 )2416 CALL wx_pbl_prelim_beta(knon, dtime, ywake_s, ybeta, &2417 BcoefQ_x, BcoefQ_w &2418 )2419 CALL wx_pbl0_merge(knon, ypplay, ypaprs, &2420 ywake_s, ydTs0, ydqs0, &2421 yt_x, yt_w, yq_x, yq_w, &2422 yu_x, yu_w, yv_x, yv_w, &2423 ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &2424 ycdragm_x, ycdragm_w, &2425 AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &2426 AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &2427 BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &2428 BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &2429 AcoefH_0, AcoefQ_0, AcoefU, AcoefV, &2430 BcoefH_0, BcoefQ_0, BcoefU, BcoefV, &2431 ycdragh, ycdragq, ycdragm, &2432 yt1, yq1, yu1, yv1 &2433 )2434 IF (iflag_split .eq. 2 .AND. nsrf .ne. is_oce) THEN2435 CALL wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, &2436 ywake_s, ybeta, ywake_cstar, ywake_dens, &2437 AcoefH_x, AcoefH_w, &2438 BcoefH_x, BcoefH_w, &2439 AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, &2440 AcoefH, AcoefQ, BcoefH, BcoefQ, &2441 HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &2442 phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &2443 yg_T, yg_Q, &2444 yGamma_dTs_phiT, yGamma_dQs_phiQ, &2445 ydTs_ins, ydqs_ins &2446 )2447 ELSE !2448 AcoefH(:) = AcoefH_0(:)2449 AcoefQ(:) = AcoefQ_0(:)2450 BcoefH(:) = BcoefH_0(:)2451 BcoefQ(:) = BcoefQ_0(:)2452 yg_T(:) = 0.2453 yg_Q(:) = 0.2454 yGamma_dTs_phiT(:) = 0.2455 yGamma_dQs_phiQ(:) = 0.2456 ydTs_ins(:) = 0.2457 ydqs_ins(:) = 0.2458 ENDIF ! (iflag_split .eq. 2)2459 ENDIF ! (iflag_split .eq.0)2460 !!!2461 IF (prt_level >=10) THEN2462 DO i = 1, min(1,knon)2463 PRINT *,'pbl_surface (merge->): yt(1,:) ',yt(i,:)2464 PRINT *,'pbl_surface (merge->): yq(1,:) ',yq(i,:)2465 PRINT *,'pbl_surface (merge->): yu(1,:) ',yu(i,:)2466 PRINT *,'pbl_surface (merge->): yv(1,:) ',yv(i,:)2467 PRINT *,'pbl_surface (merge->): AcoefH(1), AcoefQ(1), AcoefU(1), AcoefV(1) ', &2468 AcoefH(i), AcoefQ(i), AcoefU(i), AcoefV(i)2469 PRINT *,'pbl_surface (merge->): BcoefH(1), BcoefQ(1), BcoefU(1), BcoefV(1) ', &2470 BcoefH(i), BcoefQ(i), BcoefU(i), BcoefV(i)2471 ENDDO2472 2473 ENDIF2474 2475 ! Save initial value of z0h for use in evappot (z0h wiil be computed again in the surface models)2476 yz0h_old(1:knon) = yz0h(1:knon)2477 !2478 !****************************************************************************************2479 !2480 ! Calulate t2m and q2m for the case of calculation at land grid points2481 ! t2m and q2m are needed as input to ORCHIDEE2482 !2483 !****************************************************************************************2484 IF (nsrf == is_ter) THEN2485 2486 DO i = 1, knon2487 zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &2488 * (ypaprs(i,1)-ypplay(i,1))2489 ENDDO2490 2491 ! Calculate the temperature et relative humidity at 2m and the wind at 10m2492 IF (iflag_new_t2mq2m==1) THEN2493 CALL stdlevvarn(klon, knon, is_ter, zxli, &2494 yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &2495 yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &2496 yt2m, yq2m, yt10m, yq10m, yu10m, yustar, &2497 yn2mout(:, nsrf, :))2498 ELSE2499 CALL stdlevvar(klon, knon, is_ter, zxli, &2500 yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &2501 yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &2502 yt2m, yq2m, yt10m, yq10m, yu10m, yustar, ypblh, rain_f, zxtsol)2503 ENDIF2504 2505 ENDIF2506 2507 !****************************************************************************************2508 !2509 ! 10) Switch according to current surface2510 ! It is necessary to start with the continental surfaces because the ocean2511 ! needs their run-off.2512 !2513 !****************************************************************************************2514 SELECT CASE(nsrf)2515 2516 CASE(is_ter)2517 ! print*,"DEBUGTS",yts(knon/2),ylwdown(knon/2)2518 CALL surf_land(itap, dtime, date0, jour, knon, ni,&2519 rlon, rlat, yrmu0, &2520 debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, &2521 !!jyg yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&2522 yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt1, yq1,&2523 AcoefH, AcoefQ, BcoefH, BcoefQ, &2524 AcoefU, AcoefV, BcoefU, BcoefV, &2525 ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &2526 ylwdown, yq2m, yt2m, &2527 ysnow, yqsol, yagesno, ytsoil, &2528 yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,yfluxbs,&2529 yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &2530 y_flux_u1, y_flux_v1, &2531 yveget,ylai,yheight, tsurf_tersrf, tsoil_tersrf, qsurf_tersrf, tsurf_new_tersrf, &2532 cdragm_tersrf, cdragh_tersrf, &2533 swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf &2534 !GG2535 ! yveget,ylai,yheight,hice,tice,bilg_cumul, &2536 ! fcds, fcdi, dh_basal_growth, dh_basal_melt, dh_top_melt, dh_snow2sic, &2537 ! dtice_melt, dtice_snow2sic)2538 !GG2539 #ifdef ISO2540 & ,yxtrain_f, yxtsnow_f,yxt1, &2541 & yxtsnow,yxtsol,yxtevap,h1, &2542 & yrunoff_diag,yxtrunoff_diag,yRland_ice &2543 #endif2544 & )2545 2546 tsurf_tersrf(:,:) = tsurf_new_tersrf(:,:) ! for next time step2547 2548 !FC quid qd yveget ylai yheight ne sont pas definit2549 !FC yveget,ylai,yheight, &2550 IF (ifl_pbltree .ge. 1) THEN2551 CALL freinage(knon, yu, yv, yt, &2552 ! yveget,ylai, yheight,ypaprs,ypplay,y_d_u_frein,y_d_v_frein)2553 yveget,ylai, yheight,ypaprs,ypplay,y_treedrg, y_d_u_frein,y_d_v_frein)2554 ENDIF2555 2556 2557 ! Special DICE MPL 05082013 puis BOMEX2558 IF (ok_prescr_ust) THEN2559 DO j=1,knon2560 ! ysnow(:)=0.2561 ! yqsol(:)=0.2562 ! yagesno(:)=50.2563 ! ytsoil(:,:)=300.2564 ! yz0_new(:)=0.0012565 ! yevap(:)=flat/RLVTT2566 ! yfluxlat(:)=-flat2567 ! yfluxsens(:)=-fsens2568 ! yqsurf(:)=0.2569 ! ytsurf_new(:)=tg2570 ! y_dflux_t(:)=0.2571 ! y_dflux_q(:)=0.2572 y_flux_u1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yu(j,1)*ypplay(j,1)/RD/yt(j,1)2573 y_flux_v1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yv(j,1)*ypplay(j,1)/RD/yt(j,1)2574 ENDDO2575 ENDIF2576 2577 #ifdef ISOVERIF2578 DO j=1,knon2579 DO ixt=1,ntraciso2580 CALL iso_verif_noNaN(yxtevap(ixt,j), &2581 & 'pbl_surface 1056a: apres surf_land')2582 ENDDO2583 DO ixt=1,niso2584 CALL iso_verif_noNaN(yxtsol(ixt,j), &2585 & 'pbl_surface 1056b: apres surf_land')2586 ENDDO2587 ENDDO2588 #endif2589 #ifdef ISOVERIF2590 ! write(*,*) 'pbl_surface_mod 1038: sortie surf_land'2591 DO j=1,knon2592 IF (iso_eau >= 0) THEN2593 CALL iso_verif_egalite(yxtsnow(iso_eau,j), &2594 & ysnow(j),'pbl_surf_mod 1043')2595 ENDIF !if (iso_eau.gt.0) then2596 ENDDO !DO i=1,klon2597 #endif2598 2599 CASE(is_lic)2600 ! Martin2601 2602 IF (landice_opt .LT. 2) THEN2603 ! Land ice is treated by LMDZ and not by ORCHIDEE2604 CALL surf_landice(itap, dtime, knon, ni, &2605 rlon, rlat, debut, lafin, &2606 yrmu0, ylwdown, yalb, zgeo1, &2607 ysolsw, ysollw, yts, ypplay(:,1), &2608 ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt1, yq1,&2609 AcoefH, AcoefQ, BcoefH, BcoefQ, &2610 AcoefU, AcoefV, BcoefU, BcoefV, &2611 AcoefQBS, BcoefQBS, &2612 ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &2613 ysnow, yqsurf, yqsol,yqbs1, yagesno, &2614 ytsoil, yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yicesub_lic, yfluxsens,yfluxlat, &2615 yfluxbs, ytsurf_new, y_dflux_t, y_dflux_q, &2616 yzmea, yzsig, ycldt, &2617 ysnowhgt, yqsnow, ytoice, ysissnow, &2618 yalb3_new, yrunoff, &2619 y_flux_u1, y_flux_v1 &2620 #ifdef ISO2621 & ,yxtrain_f, yxtsnow_f,yxt1,yRland_ice &2622 & ,yxtsnow,yxtsol,yxtevap &2623 #endif2624 & )2625 2626 !jyg<2627 !! alb3_lic(:)=0.2628 !>jyg2629 DO j = 1, knon2630 i = ni(j)2631 alb3_lic(i) = yalb3_new(j)2632 snowhgt(i) = ysnowhgt(j)2633 qsnow(i) = yqsnow(j)2634 to_ice(i) = ytoice(j)2635 sissnow(i) = ysissnow(j)2636 runoff(i) = yrunoff(j)2637 icesub_lic(i) = yicesub_lic(j)*ypct(j)2638 ENDDO2639 ! Martin2640 ! Special DICE MPL 05082013 puis BOMEX MPL 201504102641 IF (ok_prescr_ust) THEN2642 DO j=1,knon2643 y_flux_u1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yu(j,1)*ypplay(j,1)/RD/yt(j,1)2644 y_flux_v1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yv(j,1)*ypplay(j,1)/RD/yt(j,1)2645 ENDDO2646 ENDIF2647 2648 #ifdef ISOVERIF2649 DO j=1,knon2650 DO ixt=1,ntraciso2651 CALL iso_verif_noNaN(yxtevap(ixt,j), &2652 & 'pbl_surface 1095a: apres surf_landice')2653 ENDDO2654 do ixt=1,niso2655 call iso_verif_noNaN(yxtsol(ixt,j), &2656 & 'pbl_surface 1095b: apres surf_landice')2657 enddo2658 enddo2659 #endif2660 #ifdef ISOVERIF2661 !write(*,*) 'pbl_surface_mod 1060: sortie surf_landice'2662 do j=1,knon2663 IF (iso_eau >= 0) THEN2664 CALL iso_verif_egalite(yxtsnow(iso_eau,j), &2665 & ysnow(j),'pbl_surf_mod 1064')2666 ENDIF !if (iso_eau >= 0) THEN2667 ENDDO !DO i=1,klon2668 #endif2669 2670 END IF2671 2672 CASE(is_oce)2673 2674 !GG2675 ! calculate length scale PBL2676 2677 if (iflag_leads == 1) then2678 ydthetadz = 999999.2679 ypphii = 999999.2680 ytheta = 999999.2681 2682 DO k = 1, klev2683 DO j = 1, knon2684 ytheta(j,k) = yt(j,k)*(ypplay(j,k)/1.e5)**(RD/RCPD)2685 ENDDO2686 ENDDO2687 2688 DO k = 2, klev2689 DO j = 1, knon2690 ydthetadz(j,k) = RG*( ytheta(j,k) - ytheta(j,k-1) ) / ( ypphi(j,k) - ypphi(j,k-1) )2691 ypphii(j,k) = (ypphi(j,k)+ypphi(j,k-1))/(RG*2.)2692 ENDDO2693 ENDDO2694 2695 DO j = 1, knon2696 ! print *, "ypphii(j,:)=", ypphii(j,:)2697 ! print *, "ypplay(j,:)=", ypplay(j,:)2698 ! print *, "ytheta(j,:)=", ytheta(j,:)2699 ! print *, "minloc(abs(ypphii(j,:)-300))=",2700 ! minloc(abs(ypphii(j,:)-300),1)2701 k= minloc(abs(ypphii(j,:)-300),1)2702 ydthetadz300(j)=ydthetadz(j,k)2703 ENDDO2704 end if2705 !GG2706 CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb_vis, &2707 ywindsp, yrmu0, yfder, yts, &2708 itap, dtime, jour, knon, ni, &2709 !!jyg ypplay(:,1), zgeo1/RG, ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&2710 ypplay(:,1), zgeo1(1:knon)/RG, ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt(:,1), yq(:,1),& ! ym missing init2711 AcoefH, AcoefQ, BcoefH, BcoefQ, &2712 AcoefU, AcoefV, BcoefU, BcoefV, &2713 ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &2714 ysnow, yqsurf, yagesno, &2715 yz0m, yz0h, SFRWL,yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&2716 ytsurf_new, y_dflux_t, y_dflux_q, slab_wfbils, &2717 y_flux_u1, y_flux_v1, ydelta_sst(:knon), ydelta_sal(:knon), &2718 yds_ns(:knon), ydt_ns(:knon), ydter(:knon), ydser(:knon), &2719 !GG ydt_ds(:knon), ytkt(:knon), ytks(:knon), ytaur(:knon), ysss)2720 ydt_ds(:knon), ytkt(:knon), ytks(:knon), ytaur(:knon), ysss, &2721 ydthetadz300,Ampl &2722 !GG2723 #ifdef ISO2724 & ,yxtrain_f, yxtsnow_f,yxt1,Roce, &2725 & yxtsnow,yxtevap,h1 &2726 #endif2727 & )2728 IF (prt_level >=10) THEN2729 print *,'arg de surf_ocean: ycdragh ',ycdragh(1:knon)2730 print *,'arg de surf_ocean: ycdragm ',ycdragm(1:knon)2731 print *,'arg de surf_ocean: yt ', yt(1:knon,:)2732 print *,'arg de surf_ocean: yq ', yq(1:knon,:)2733 print *,'arg de surf_ocean: yts ', yts(1:knon)2734 print *,'arg de surf_ocean: AcoefH ',AcoefH(1:knon)2735 print *,'arg de surf_ocean: AcoefQ ',AcoefQ(1:knon)2736 print *,'arg de surf_ocean: BcoefH ',BcoefH(1:knon)2737 print *,'arg de surf_ocean: BcoefQ ',BcoefQ(1:knon)2738 print *,'arg de surf_ocean: yevap ',yevap(1:knon)2739 print *,'arg de surf_ocean: yfluxsens ',yfluxsens(1:knon)2740 print *,'arg de surf_ocean: yfluxlat ',yfluxlat(1:knon)2741 print *,'arg de surf_ocean: ytsurf_new ',ytsurf_new(1:knon)2742 ENDIF2743 ! Special DICE MPL 05082013 puis BOMEX MPL 201504102744 IF (ok_prescr_ust) THEN2745 DO j=1,knon2746 y_flux_u1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yu(j,1)*ypplay(j,1)/RD/yt(j,1)2747 y_flux_v1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yv(j,1)*ypplay(j,1)/RD/yt(j,1)2748 ENDDO2749 ENDIF2750 2751 CASE(is_sic)2752 CALL surf_seaice( &2753 !albedo SB >>>2754 rlon, rlat, ysolsw, ysollw, yalb_vis, yfder, &2755 !albedo SB <<<2756 itap, dtime, jour, knon, ni, &2757 lafin, &2758 !!jyg yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&2759 yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&2760 AcoefH, AcoefQ, BcoefH, BcoefQ, &2761 AcoefU, AcoefV, BcoefU, BcoefV, &2762 ypsref, yu1, yv1, ygustiness, pctsrf, &2763 ysnow, yqsurf, yqsol, yagesno, ytsoil, &2764 !albedo SB >>>2765 yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&2766 !albedo SB <<<2767 ytsurf_new, y_dflux_t, y_dflux_q, &2768 !GG y_flux_u1, y_flux_v1)2769 y_flux_u1, y_flux_v1, &2770 hice,tice,bilg_cumul, &2771 fcds, fcdi, dh_basal_growth, dh_basal_melt, dh_top_melt, dh_snow2sic, &2772 dtice_melt, dtice_snow2sic &2773 !GG2774 #ifdef ISO2775 & ,yxtrain_f, yxtsnow_f,yxt1,Roce, &2776 & yxtsnow,yxtsol,yxtevap,Rland_ice &2777 #endif2778 & )2779 2780 ! Special DICE MPL 05082013 puis BOMEX MPL 201504102781 IF (ok_prescr_ust) THEN2782 DO j=1,knon2783 y_flux_u1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yu(j,1)*ypplay(j,1)/RD/yt(j,1)2784 y_flux_v1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yv(j,1)*ypplay(j,1)/RD/yt(j,1)2785 ENDDO2786 ENDIF2787 2788 #ifdef ISOVERIF2789 DO j=1,knon2790 DO ixt=1,ntraciso2791 CALL iso_verif_noNaN(yxtevap(ixt,j), &2792 & 'pbl_surface 1165a: apres surf_seaice')2793 ENDDO2794 DO ixt=1,niso2795 CALL iso_verif_noNaN(yxtsol(ixt,j), &2796 & 'pbl_surface 1165b: apres surf_seaice')2797 ENDDO2798 ENDDO2799 #endif2800 #ifdef ISOVERIF2801 !write(*,*) 'pbl_surface_mod 1077: sortie surf_seaice'2802 DO j=1,knon2803 IF (iso_eau >= 0) THEN2804 CALL iso_verif_egalite(yxtsnow(iso_eau,j), &2805 & ysnow(j),'pbl_surf_mod 1106')2806 ENDIF !IF (iso_eau >= 0) THEN2807 ENDDO !DO i=1,klon2808 #endif2809 2810 CASE DEFAULT2811 WRITE(lunout,*) 'Surface index = ', nsrf2812 abort_message = 'Surface index not valid'2813 CALL abort_physic(modname,abort_message,1)2814 END SELECT2815 2816 2817 !****************************************************************************************2818 ! 11) - Calcul the increment of surface temperature2819 !2820 !****************************************************************************************2821 2822 IF (evap0>=0.) THEN2823 yevap(1:knon)=evap02824 yevap(1:knon)=RLVTT*evap02825 ENDIF2826 2827 y_d_ts(1:knon) = ytsurf_new(1:knon) - yts(1:knon)2828 2829 !****************************************************************************************2830 !2831 ! 12) "La remontee" - "The uphill"2832 !2833 ! The fluxes (y_flux_X) and tendancy (y_d_X) are calculated2834 ! for X=H, Q, U and V, for all vertical levels.2835 !2836 !****************************************************************************************2837 !!2838 !!!2839 !!! jyg le 10/04/2013 et EV 10/20202840 2841 IF (ok_forc_tsurf) THEN2842 DO j=1,knon2843 ytsurf_new(j)=tg2844 y_d_ts(j) = ytsurf_new(j) - yts(j)2845 ENDDO2846 ENDIF ! ok_forc_tsurf2847 2848 !!!2849 IF (ok_flux_surf) THEN2850 IF (prt_level >=10) THEN2851 PRINT *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT2852 ENDIF2853 y_flux_t1(:) = fsens2854 y_flux_q1(:) = flat/RLVTT2855 yfluxlat(:) = flat2856 !2857 !! Test sur iflag_split retire le 2/02/2018, sans vraiment comprendre la raison de ce test. (jyg)2858 !! IF (iflag_split .eq.0) THEN2859 DO j=1,knon2860 Kech_h(j) = ycdragh(j) * (1.0+SQRT(yu(j,1)**2+yv(j,1)**2)) * &2861 ypplay(j,1)/(RD*yt(j,1))2862 ENDDO2863 !! ENDIF ! (iflag_split .eq.0)2864 2865 DO j = 1, knon2866 yt1_new=(1./RCPD)*(AcoefH(j)+BcoefH(j)*y_flux_t1(j)*dtime)2867 ytsurf_new(j)=yt1_new-y_flux_t1(j)/(Kech_h(j)*RCPD)2868 ! for cases forced in flux and for which forcing in Ts is needed2869 ! to prevent the latter to reach unrealistic value (even if not used,2870 ! Ts is calculated and hgardfou can appear during the calculation2871 ! of surface saturation humidity for example2872 if (ok_forc_tsurf) ytsurf_new(j)=tg2873 ENDDO2874 2875 DO j=1,knon2876 y_d_ts(j) = ytsurf_new(j) - yts(j)2877 ENDDO2878 2879 ELSE ! (ok_flux_surf)2880 DO j=1,knon2881 y_flux_t1(j) = yfluxsens(j)2882 y_flux_q1(j) = -yevap(j)2883 #ifdef ISO2884 y_flux_xt1(:,:) = -yxtevap(:,:)2885 #endif2886 ENDDO2887 ENDIF ! (ok_flux_surf)2888 2889 ! flux of blowing snow at the first level2890 IF (ok_bs) THEN2891 DO j=1,knon2892 y_flux_bs(j)=yfluxbs(j)2893 ENDDO2894 ENDIF2895 !2896 ! ------------------------------------------------------------------------------2897 ! 12a) Splitting2898 ! ------------------------------------------------------------------------------2899 2900 IF (iflag_split .GE. 1) THEN2901 #ifdef ISO2902 call abort_physic('pbl_surface_mod 2607','isos pas encore dans iflag_split=1',1)2903 #endif2904 !2905 !2906 IF (nsrf .ne. is_oce) THEN2907 !2908 ! Compute potential evaporation and aridity factor (jyg, 20200328)2909 ybeta_prev(:) = ybeta(:)2910 DO j = 1, knon2911 yqa(j) = AcoefQ(j) - BcoefQ(j)*yevap(j)*dtime2912 ENDDO2913 !2914 CALL wx_evappot(knon, yqa, yTsurf_new, yevap_pot)2915 !2916 ybeta(1:knon) = min(yevap(1:knon)/yevap_pot(1:knon), 1.)2917 2918 IF (prt_level >=10) THEN2919 DO j=1,knon2920 print*,'y_flux_t1,yfluxlat,wakes' &2921 & , y_flux_t1(j), yfluxlat(j), ywake_s(j)2922 print*,'beta_prev, beta, ytsurf_new', ybeta_prev(j), ybeta(j), ytsurf_new(j)2923 print*,'inertia,facteur,cstar', inertia, facteur,wake_cstar(j)2924 ENDDO2925 ENDIF ! (prt_level >=10)2926 !2927 ! Second call to wx_pbl0_merge and wx_pbl_dts_merge in order to take into account2928 ! the update of the aridity coeficient beta.2929 !2930 CALL wx_pbl_prelim_beta(knon, dtime, ywake_s, ybeta, &2931 BcoefQ_x, BcoefQ_w &2932 )2933 CALL wx_pbl0_merge(knon, ypplay, ypaprs, &2934 ywake_s, ydTs0, ydqs0, &2935 yt_x, yt_w, yq_x, yq_w, &2936 yu_x, yu_w, yv_x, yv_w, &2937 ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &2938 ycdragm_x, ycdragm_w, &2939 AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &2940 AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &2941 BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &2942 BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &2943 AcoefH_0, AcoefQ_0, AcoefU, AcoefV, &2944 BcoefH_0, BcoefQ_0, BcoefU, BcoefV, &2945 ycdragh, ycdragq, ycdragm, &2946 yt1, yq1, yu1, yv1 &2947 )2948 IF (iflag_split .eq. 2) THEN2949 CALL wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, &2950 ywake_s, ybeta, ywake_cstar, ywake_dens, &2951 AcoefH_x, AcoefH_w, &2952 BcoefH_x, BcoefH_w, &2953 AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, &2954 AcoefH, AcoefQ, BcoefH, BcoefQ, &2955 HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &2956 phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &2957 yg_T, yg_Q, &2958 yGamma_dTs_phiT, yGamma_dQs_phiQ, &2959 ydTs_ins, ydqs_ins &2960 )2961 ELSE !2962 AcoefH(:) = AcoefH_0(:)2963 AcoefQ(:) = AcoefQ_0(:)2964 BcoefH(:) = BcoefH_0(:)2965 BcoefQ(:) = BcoefQ_0(:)2966 yg_T(:) = 0.2967 yg_Q(:) = 0.2968 yGamma_dTs_phiT(:) = 0.2969 yGamma_dQs_phiQ(:) = 0.2970 ydTs_ins(:) = 0.2971 ydqs_ins(:) = 0.2972 ENDIF ! (iflag_split .eq. 2)2973 !2974 ELSE ! (nsrf .ne. is_oce)2975 ybeta(1:knon) = 1.2976 yevap_pot(1:knon) = yevap(1:knon)2977 AcoefH(:) = AcoefH_0(:)2978 AcoefQ(:) = AcoefQ_0(:)2979 BcoefH(:) = BcoefH_0(:)2980 BcoefQ(:) = BcoefQ_0(:)2981 yg_T(:) = 0.2982 yg_Q(:) = 0.2983 yGamma_dTs_phiT(:) = 0.2984 yGamma_dQs_phiQ(:) = 0.2985 ydTs_ins(:) = 0.2986 ydqs_ins(:) = 0.2987 ENDIF ! (nsrf .ne. is_oce)2988 !2989 CALL wx_pbl_split(knon, nsrf, dtime, ywake_s, ybeta, iflag_split, &2990 yg_T, yg_Q, &2991 yGamma_dTs_phiT, yGamma_dQs_phiQ, &2992 ydTs_ins, ydqs_ins, &2993 y_flux_t1, y_flux_q1, y_flux_u1, y_flux_v1, &2994 !!!! HTRn_b, dd_HTRn, HTphiT_b, dd_HTphiT, &2995 phiQ0_b, phiT0_b, &2996 y_flux_t1_x, y_flux_t1_w, &2997 y_flux_q1_x, y_flux_q1_w, &2998 y_flux_u1_x, y_flux_u1_w, &2999 y_flux_v1_x, y_flux_v1_w, &3000 yfluxlat_x, yfluxlat_w, &3001 y_delta_qsats, &3002 y_delta_tsurf_new, y_delta_qsurf &3003 )3004 !3005 CALL wx_pbl_check(knon, dtime, ypplay, ypaprs, ywake_s, ybeta, iflag_split, &3006 yTs, y_delta_tsurf, &3007 yqsurf, yTsurf_new, &3008 y_delta_tsurf_new, y_delta_qsats, &3009 AcoefH_x, AcoefH_w, &3010 BcoefH_x, BcoefH_w, &3011 AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, &3012 AcoefH, AcoefQ, BcoefH, BcoefQ, &3013 y_flux_t1, y_flux_q1, &3014 y_flux_t1_x, y_flux_t1_w, &3015 y_flux_q1_x, y_flux_q1_w)3016 !3017 IF (nsrf .ne. is_oce) THEN3018 CALL wx_pbl_dts_check(knon, dtime, ypplay, ypaprs, ywake_s, ybeta, iflag_split, &3019 yTs, y_delta_tsurf, &3020 yqsurf, yTsurf_new, &3021 y_delta_qsats, y_delta_tsurf_new, y_delta_qsurf, &3022 AcoefH_x, AcoefH_w, &3023 BcoefH_x, BcoefH_w, &3024 AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, &3025 AcoefH, AcoefQ, BcoefH, BcoefQ, &3026 HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &3027 phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &3028 yg_T, yg_Q, &3029 yGamma_dTs_phiT, yGamma_dQs_phiQ, &3030 ydTs_ins, ydqs_ins, &3031 y_flux_t1, y_flux_q1, &3032 y_flux_t1_x, y_flux_t1_w, &3033 y_flux_q1_x, y_flux_q1_w )3034 ENDIF ! (nsrf .ne. is_oce)3035 !3036 ELSE ! (iflag_split .ge. 1)3037 ybeta(1:knon) = 1.3038 yevap_pot(1:knon) = yevap(1:knon)3039 ENDIF ! (iflag_split .ge. 1)3040 !3041 IF (prt_level >= 10) THEN3042 print *,'pbl_surface, ybeta , yevap, yevap_pot ', &3043 ybeta(1:knon) , yevap(1:knon), yevap_pot(1:knon)3044 ENDIF ! (prt_level >= 10)3045 !3046 !>jyg3047 !3048 3049 !!jyg!! A reprendre apres reflexion ===============================================3050 !!jyg!!3051 !!jyg!! DO j=1,knon3052 !!jyg!!!!! nrlmd le 13/06/20113053 !!jyg!!3054 !!jyg!!!----Diffusion dans le sol dans le cas continental seulement3055 !!jyg!! IF (nsrf.eq.is_ter) THEN3056 !!jyg!!!----Calcul du coefficient delta_coeff3057 !!jyg!! tau_eq(j)=(ywake_s(j)/2.)*(1./max(wake_cstar(j),0.01))*sqrt(0.4/(3.14*max(wake_dens(j),8e-12)))3058 !!jyg!!3059 !!jyg!!! delta_coef(j)=dtime/(inertia*sqrt(tau_eq(j)))3060 !!jyg!! delta_coef(j)=facteur*sqrt(tau_eq(j))/inertia3061 !!jyg!!! delta_coef(j)=0.3062 !!jyg!! ELSE3063 !!jyg!! delta_coef(j)=0.3064 !!jyg!! ENDIF3065 !!jyg!!3066 !!jyg!!!----Calcul de delta_tsurf3067 !!jyg!! y_delta_tsurf(j)=delta_coef(j)*y_delta_flux_t1(j)3068 !!jyg!!3069 !!jyg!!!----Si il n'y a pas des poches...3070 !!jyg!! IF (wake_cstar(j).le.0.01) THEN3071 !!jyg!! y_delta_tsurf(j)=0.3072 !!jyg!! y_delta_flux_t1(j)=0.3073 !!jyg!! ENDIF3074 !!jyg!!3075 !!jyg!!!-----Calcul de ybeta (evap_r\'eelle/evap_potentielle)3076 !!jyg!!!!!!! jyg le 23/02/20123077 !!jyg!!!!!!!3078 !!jyg!!!! ybeta(j)=y_flux_q1(j) / &3079 !!jyg!!!! & (Kech_h(j)*(yq(j,1)-yqsatsurf(j)))3080 !!jyg!!!!!! ybeta(j)=-1.*yevap(j) / &3081 !!jyg!!!!!! & (ywake_s(j)*Kech_h_w(j)*(yq_w(j,1)-yqsatsurf_w(j))+(1.-ywake_s(j))*Kech_h_x(j)*(yq_x(j,1)-yqsatsurf_x(j)))3082 !!jyg!!!!!!! fin jyg3083 !!jyg!!!!!3084 !!jyg!!3085 !!jyg!! ENDDO3086 !!jyg!!3087 !!jyg!!!!! fin nrlmd le 13/06/20113088 !!jyg!!3089 IF (iflag_split .ge. 1) THEN3090 IF (prt_level >=10) THEN3091 DO j = 1, knon3092 print*,'Chx,Chw,Ch', ycdragh_x(j), ycdragh_w(j), ycdragh(j)3093 print*,'Khx,Khw,Kh', Kech_h_x(j), Kech_h_w(j), Kech_h(j)3094 print*,'t1x, t1w, t1, t1_ancien', &3095 & yt_x(j,1), yt_w(j,1), yt(j,1), t(j,1)3096 print*,'delta_coef,delta_flux,delta_tsurf,tau', delta_coef(j), y_delta_flux_t1(j), y_delta_tsurf(j), tau_eq(j)3097 ENDDO3098 3099 DO j=1,knon3100 print*,'fluxT_x, fluxT_w, y_flux_t1, fluxQ_x, fluxQ_w, yfluxlat, wakes' &3101 & , 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)3102 print*,'beta, ytsurf_new ', ybeta(j), ytsurf_new(j)3103 print*,'inertia, facteur, cstar', inertia, facteur,wake_cstar(j)3104 ENDDO3105 ENDIF ! (prt_level >=10)3106 3107 !!! jyg le 07/02/20123108 ENDIF ! (iflag_split .ge.1)3109 !!!3110 3111 !!! jyg le 07/02/20123112 IF (iflag_split .eq.0) THEN3113 !!!3114 !!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/20123115 CALL climb_hq_up(knon, ni, dtime, yt, yq, &3116 y_flux_q1, y_flux_t1, ypaprs, ypplay, &3117 !!! jyg le 07/02/20123118 AcoefH, AcoefQ, BcoefH, BcoefQ, &3119 CcoefH, CcoefQ, DcoefH, DcoefQ, &3120 Kcoef_hq, gama_q, gama_h, &3121 !!!3122 y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:) &3123 #ifdef ISO3124 & ,yxt,y_flux_xt1 &3125 & ,AcoefXT,BcoefXT,CcoefXT,DcoefXT,gama_xt &3126 & ,y_flux_xt(:,:,:),y_d_xt(:,:,:) &3127 #endif3128 & )3129 ELSE !(iflag_split .eq.0)3130 CALL climb_hq_up(knon, ni, dtime, yt_x, yq_x, &3131 y_flux_q1_x, y_flux_t1_x, ypaprs, ypplay, &3132 !!! nrlmd le 02/05/20113133 AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x, &3134 CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x, &3135 Kcoef_hq_x, gama_q_x, gama_h_x, &3136 !!!3137 y_flux_q_x(:,:), y_flux_t_x(:,:), y_d_q_x(:,:), y_d_t_x(:,:) &3138 #ifdef ISO3139 & ,yxt_x,y_flux_xt1_x &3140 & ,AcoefXT_x,BcoefXT_x,CcoefXT_x,DcoefXT_x,gama_xt_x &3141 & ,y_flux_xt_x(:,:,:),y_d_xt_x(:,:,:) &3142 #endif3143 & )3144 !3145 CALL climb_hq_up(knon, ni, dtime, yt_w, yq_w, &3146 y_flux_q1_w, y_flux_t1_w, ypaprs, ypplay, &3147 !!! nrlmd le 02/05/20113148 AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w, &3149 CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w, &3150 Kcoef_hq_w, gama_q_w, gama_h_w, &3151 !!!3152 y_flux_q_w(:,:), y_flux_t_w(:,:), y_d_q_w(:,:), y_d_t_w(:,:) &3153 #ifdef ISO3154 & ,yxt_w,y_flux_xt1_w &3155 & ,AcoefXT_w,BcoefXT_w,CcoefXT_w,DcoefXT_w,gama_xt_w &3156 & ,y_flux_xt_w(:,:,:),y_d_xt_w(:,:,:) &3157 #endif3158 & )3159 !!!3160 ENDIF ! (iflag_split .eq.0)3161 !!!3162 3163 !!! jyg le 07/02/20123164 IF (iflag_split .eq.0) THEN3165 !!!3166 !!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/20123167 CALL climb_wind_up(knon, ni, dtime, yu, yv, y_flux_u1, y_flux_v1, &3168 !!! jyg le 07/02/20123169 AcoefU, AcoefV, BcoefU, BcoefV, &3170 CcoefU, CcoefV, DcoefU, DcoefV, &3171 Kcoef_m, &3172 !!!3173 y_flux_u, y_flux_v, y_d_u, y_d_v)3174 y_d_t_diss(:,:)=0.3175 IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN3176 CALL yamada_c(knon, knon,dtime,ypaprs,ypplay &3177 & ,yu,yv,yt,y_d_u,y_d_v,y_d_t,ycdragm,ytke,ycoefm,ycoefh,ycoefq,y_d_t_diss,yustar &3178 & ,iflag_pbl)3179 ENDIF3180 ! print*,'yamada_c OK'3181 3182 ELSE !(iflag_split .eq.0)3183 CALL climb_wind_up(knon, ni, dtime, yu_x, yv_x, y_flux_u1_x, y_flux_v1_x, &3184 !!! nrlmd le 02/05/20113185 AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x, &3186 CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x, &3187 Kcoef_m_x, &3188 !!!3189 y_flux_u_x, y_flux_v_x, y_d_u_x, y_d_v_x)3190 !3191 y_d_t_diss_x(:,:)=0.3192 IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN3193 CALL yamada_c(knon, knon,dtime,ypaprs,ypplay &3194 & ,yu_x,yv_x,yt_x,y_d_u_x,y_d_v_x,y_d_t_x,ycdragm_x,ytke_x,ycoefm_x,ycoefh_x &3195 ,ycoefq_x,y_d_t_diss_x,yustar_x &3196 & ,iflag_pbl)3197 ENDIF3198 ! print*,'yamada_c OK'3199 3200 CALL climb_wind_up(knon, ni, dtime, yu_w, yv_w, y_flux_u1_w, y_flux_v1_w, &3201 !!! nrlmd le 02/05/20113202 AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w, &3203 CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w, &3204 Kcoef_m_w, &3205 !!!3206 y_flux_u_w, y_flux_v_w, y_d_u_w, y_d_v_w)3207 !!!3208 y_d_t_diss_w(:,:)=0.3209 IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN3210 CALL yamada_c(knon, knon,dtime,ypaprs,ypplay &3211 & ,yu_w,yv_w,yt_w,y_d_u_w,y_d_v_w,y_d_t_w,ycdragm_w,ytke_w,ycoefm_w,ycoefh_w &3212 ,ycoefq_w,y_d_t_diss_w,yustar_w &3213 & ,iflag_pbl)3214 ENDIF3215 ! print*,'yamada_c OK'3216 !3217 IF (prt_level >=10) THEN3218 print *, 'After climbing up, lfuxlat_x, fluxlat_w ', &3219 yfluxlat_x(1:knon), yfluxlat_w(1:knon)3220 ENDIF3221 !3222 ENDIF ! (iflag_split .eq.0)3223 3224 IF (ok_bs) THEN3225 CALL climb_qbs_up(knon, ni, dtime, yqbs, &3226 y_flux_bs, ypaprs, ypplay, &3227 AcoefQBS, BcoefQBS, &3228 CcoefQBS, DcoefQBS, &3229 Kcoef_qbs, gama_qbs, &3230 y_flux_qbs(:,:), y_d_qbs(:,:))3231 ENDIF3232 3233 !!!3234 !!3235 !! DO j = 1, knon3236 !! y_dflux_t(j) = y_dflux_t(j) * ypct(j)3237 !! y_dflux_q(j) = y_dflux_q(j) * ypct(j)3238 !! ENDDO3239 !!3240 !****************************************************************************************3241 ! 13) Transform variables for output format :3242 ! - Decompress3243 ! - Multiply with pourcentage of current surface3244 ! - Cumulate in global variable3245 !3246 !****************************************************************************************3247 3248 3249 !!! jyg le 07/02/20123250 IF (iflag_split.EQ.0) THEN3251 !!!3252 DO k = 1, klev3253 DO j = 1, knon3254 i = ni(j)3255 y_d_t_diss(j,k) = y_d_t_diss(j,k) * ypct(j)3256 y_d_t(j,k) = y_d_t(j,k) * ypct(j)3257 y_d_q(j,k) = y_d_q(j,k) * ypct(j)3258 y_d_u(j,k) = y_d_u(j,k) * ypct(j)3259 y_d_v(j,k) = y_d_v(j,k) * ypct(j)3260 !FC3261 IF (nsrf .EQ. is_ter .and. ifl_pbltree .GE. 1) THEN3262 ! if (y_d_u_frein(j,k).ne.0. ) then3263 ! print*, nsrf,'IS_TER ++', y_d_u_frein(j,k)*ypct(j),y_d_u(j,k),j,k3264 ! ENDIF3265 y_d_u(j,k) =y_d_u(j,k) + y_d_u_frein(j,k)*ypct(j)3266 y_d_v(j,k) =y_d_v(j,k) + y_d_v_frein(j,k)*ypct(j)3267 treedrg(i,k,nsrf)=y_treedrg(j,k)3268 ELSE3269 treedrg(i,k,nsrf)=0.3270 ENDIF3271 !FC3272 flux_t(i,k,nsrf) = y_flux_t(j,k)3273 flux_q(i,k,nsrf) = y_flux_q(j,k)3274 flux_u(i,k,nsrf) = y_flux_u(j,k)3275 flux_v(i,k,nsrf) = y_flux_v(j,k)3276 3277 #ifdef ISO3278 DO ixt=1,ntraciso3279 y_d_xt(ixt,j,k) = y_d_xt(ixt,j,k) * ypct(j)3280 flux_xt(ixt,i,k,nsrf) = y_flux_xt(ixt,j,k)3281 ENDDO ! DO ixt=1,ntraciso3282 h1_diag(i)=h1(j)3283 #endif3284 3285 ENDDO3286 ENDDO3287 3288 #ifdef ISO3289 #ifdef ISOVERIF3290 if (iso_eau.gt.0) then3291 call iso_verif_egalite_vect2D( &3292 y_d_xt,y_d_q, &3293 'pbl_surface_mod 2600',ntraciso,klon,klev)3294 endif3295 #endif3296 #endif3297 3298 ELSE !(iflag_split .eq.0)3299 3300 ! Tendances hors poches3301 DO k = 1, klev3302 DO j = 1, knon3303 i = ni(j)3304 y_d_t_diss_x(j,k) = y_d_t_diss_x(j,k) * ypct(j)3305 y_d_t_x(j,k) = y_d_t_x(j,k) * ypct(j)3306 y_d_q_x(j,k) = y_d_q_x(j,k) * ypct(j)3307 y_d_u_x(j,k) = y_d_u_x(j,k) * ypct(j)3308 y_d_v_x(j,k) = y_d_v_x(j,k) * ypct(j)3309 3310 flux_t_x(i,k,nsrf) = y_flux_t_x(j,k)3311 flux_q_x(i,k,nsrf) = y_flux_q_x(j,k)3312 flux_u_x(i,k,nsrf) = y_flux_u_x(j,k)3313 flux_v_x(i,k,nsrf) = y_flux_v_x(j,k)3314 3315 #ifdef ISO3316 DO ixt=1,ntraciso3317 y_d_xt_x(ixt,j,k) = y_d_xt_x(ixt,j,k) * ypct(j)3318 flux_xt_x(ixt,i,k,nsrf) = y_flux_xt_x(ixt,j,k)3319 ENDDO ! DO ixt=1,ntraciso3320 #endif3321 ENDDO3322 ENDDO3323 3324 ! Tendances dans les poches3325 DO k = 1, klev3326 DO j = 1, knon3327 i = ni(j)3328 y_d_t_diss_w(j,k) = y_d_t_diss_w(j,k) * ypct(j)3329 y_d_t_w(j,k) = y_d_t_w(j,k) * ypct(j)3330 y_d_q_w(j,k) = y_d_q_w(j,k) * ypct(j)3331 y_d_u_w(j,k) = y_d_u_w(j,k) * ypct(j)3332 y_d_v_w(j,k) = y_d_v_w(j,k) * ypct(j)3333 3334 flux_t_w(i,k,nsrf) = y_flux_t_w(j,k)3335 flux_q_w(i,k,nsrf) = y_flux_q_w(j,k)3336 flux_u_w(i,k,nsrf) = y_flux_u_w(j,k)3337 flux_v_w(i,k,nsrf) = y_flux_v_w(j,k)3338 3339 #ifdef ISO3340 DO ixt=1,ntraciso3341 y_d_xt_w(ixt,j,k) = y_d_xt_w(ixt,j,k) * ypct(j)3342 flux_xt_w(ixt,i,k,nsrf) = y_flux_xt_w(ixt,j,k)3343 ENDDO ! do ixt=1,ntraciso3344 #endif3345 3346 ENDDO3347 ENDDO3348 3349 ! Flux, tendances et Tke moyenne dans la maille3350 DO k = 1, klev3351 DO j = 1, knon3352 i = ni(j)3353 flux_t(i,k,nsrf) = flux_t_x(i,k,nsrf)+ywake_s(j)*(flux_t_w(i,k,nsrf)-flux_t_x(i,k,nsrf))3354 flux_q(i,k,nsrf) = flux_q_x(i,k,nsrf)+ywake_s(j)*(flux_q_w(i,k,nsrf)-flux_q_x(i,k,nsrf))3355 flux_u(i,k,nsrf) = flux_u_x(i,k,nsrf)+ywake_s(j)*(flux_u_w(i,k,nsrf)-flux_u_x(i,k,nsrf))3356 flux_v(i,k,nsrf) = flux_v_x(i,k,nsrf)+ywake_s(j)*(flux_v_w(i,k,nsrf)-flux_v_x(i,k,nsrf))3357 #ifdef ISO3358 DO ixt=1,ntraciso3359 flux_xt(ixt,i,k,nsrf) = flux_xt_x(ixt,i,k,nsrf)+ywake_s(j)*(flux_xt_w(ixt,i,k,nsrf)-flux_xt_x(ixt,i,k,nsrf))3360 ENDDO ! do ixt=1,ntraciso3361 #endif3362 ENDDO3363 ENDDO3364 DO j=1,knon3365 yfluxlat(j)=yfluxlat_x(j)+ywake_s(j)*(yfluxlat_w(j)-yfluxlat_x(j))3366 ENDDO3367 IF (prt_level >=10) THEN3368 print *,' nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf) ', &3369 nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf)3370 ENDIF3371 3372 DO k = 1, klev3373 DO j = 1, knon3374 y_d_t_diss(j,k) = y_d_t_diss_x(j,k)+ywake_s(j)*(y_d_t_diss_w(j,k) -y_d_t_diss_x(j,k))3375 y_d_t(j,k) = y_d_t_x(j,k)+ywake_s(j)*(y_d_t_w(j,k) -y_d_t_x(j,k))3376 y_d_q(j,k) = y_d_q_x(j,k)+ywake_s(j)*(y_d_q_w(j,k) -y_d_q_x(j,k))3377 y_d_u(j,k) = y_d_u_x(j,k)+ywake_s(j)*(y_d_u_w(j,k) -y_d_u_x(j,k))3378 y_d_v(j,k) = y_d_v_x(j,k)+ywake_s(j)*(y_d_v_w(j,k) -y_d_v_x(j,k))3379 ENDDO3380 ENDDO3381 3382 ENDIF ! (iflag_split .eq.0)3383 3384 3385 ! tendencies of blowing snow3386 IF (ok_bs) THEN3387 DO k = 1, klev3388 DO j = 1, knon3389 i = ni(j)3390 y_d_qbs(j,k)=y_d_qbs(j,k) * ypct(j)3391 flux_qbs(i,k,nsrf) = y_flux_qbs(j,k)3392 ENDDO3393 ENDDO3394 ENDIF3395 3396 3397 DO j = 1, knon3398 i = ni(j)3399 evap(i,nsrf) = - flux_q(i,1,nsrf) !jyg3400 if (ok_bs) then ; snowerosion(i,nsrf)=flux_qbs(i,1,nsrf); endif3401 beta(i,nsrf) = ybeta(j) !jyg3402 d_ts(i,nsrf) = y_d_ts(j)3403 !albedo SB >>>3404 DO k=1,nsw3405 alb_dir(i,k,nsrf) = yalb_dir_new(j,k)3406 alb_dif(i,k,nsrf) = yalb_dif_new(j,k)3407 ENDDO3408 !albedo SB <<<3409 snow(i,nsrf) = ysnow(j)3410 qsurf(i,nsrf) = yqsurf(j)3411 z0m(i,nsrf) = yz0m(j)3412 z0h(i,nsrf) = yz0h(j)3413 fluxlat(i,nsrf) = yfluxlat(j)3414 agesno(i,nsrf) = yagesno(j)3415 cdragh(i) = cdragh(i) + ycdragh(j)*ypct(j)3416 cdragm(i) = cdragm(i) + ycdragm(j)*ypct(j)3417 dflux_t(i) = dflux_t(i) + y_dflux_t(j)*ypct(j)3418 dflux_q(i) = dflux_q(i) + y_dflux_q(j)*ypct(j)3419 #ifdef ISO3420 DO ixt=1,niso3421 xtsnow(ixt,i,nsrf) = yxtsnow(ixt,j)3422 ENDDO3423 DO ixt=1,ntraciso3424 xtevap(ixt,i,nsrf) = - flux_xt(ixt,i,1,nsrf)3425 dflux_xt(ixt,i) = dflux_xt(ixt,i) + y_dflux_xt(ixt,j)*ypct(j)3426 ENDDO3427 IF (nsrf == is_lic) THEN3428 DO ixt=1,niso3429 Rland_ice(ixt,i) = yRland_ice(ixt,j)3430 ENDDO3431 ENDIF !IF (nsrf == is_lic) THEN3432 #ifdef ISOVERIF3433 IF (iso_eau.gt.0) THEN3434 call iso_verif_egalite_choix(Rland_ice(iso_eau,i),1.0, &3435 & 'pbl_surf_mod 1230',errmax,errmaxrel)3436 ENDIF !if (iso_eau.gt.0) then3437 #endif3438 #endif3439 ENDDO3440 3441 ! print*,'Dans pbl OK2'3442 3443 !!! jyg le 07/02/20123444 IF (iflag_split .ge.1) THEN3445 !!!3446 !!! nrlmd le 02/05/20113447 DO j = 1, knon3448 i = ni(j)3449 fluxlat_x(i,nsrf) = yfluxlat_x(j)3450 fluxlat_w(i,nsrf) = yfluxlat_w(j)3451 !!!3452 !!! nrlmd le 13/06/20113453 !!jyg20170131 delta_tsurf(i,nsrf)=y_delta_tsurf(j)*ypct(j)3454 !!jyg20210118 delta_tsurf(i,nsrf)=y_delta_tsurf(j)3455 delta_tsurf(i,nsrf)=y_delta_tsurf_new(j)3456 !3457 delta_qsurf(i) = delta_qsurf(i) + y_delta_qsurf(j)*ypct(j)3458 !3459 cdragh_x(i) = cdragh_x(i) + ycdragh_x(j)*ypct(j)3460 cdragh_w(i) = cdragh_w(i) + ycdragh_w(j)*ypct(j)3461 cdragm_x(i) = cdragm_x(i) + ycdragm_x(j)*ypct(j)3462 cdragm_w(i) = cdragm_w(i) + ycdragm_w(j)*ypct(j)3463 kh(i) = kh(i) + Kech_h(j)*ypct(j)3464 kh_x(i) = kh_x(i) + Kech_h_x(j)*ypct(j)3465 kh_w(i) = kh_w(i) + Kech_h_w(j)*ypct(j)3466 !!!3467 ENDDO3468 !!!3469 ENDIF ! (iflag_split .ge.1)3470 !!!3471 !!! nrlmd le 02/05/20113472 !!jyg le 20/02/20113473 !! tke_x(:,:,nsrf)=0.3474 !! tke_w(:,:,nsrf)=0.3475 !!jyg le 20/02/20113476 !! DO k = 1, klev+13477 !! DO j = 1, knon3478 !! i = ni(j)3479 !! wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k)3480 !! tke(i,k,nsrf) = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)3481 !! ENDDO3482 !! ENDDO3483 !!jyg le 20/02/20113484 !! DO k = 1, klev+13485 !! DO j = 1, knon3486 !! i = ni(j)3487 !! tke(i,k,nsrf)=(1.-ywake_s(j))*tke_x(i,k,nsrf)+ywake_s(j)*tke_w(i,k,nsrf)3488 !! ENDDO3489 !! ENDDO3490 !!!3491 IF (iflag_split .eq.0) THEN3492 wake_dltke(:,:,nsrf) = 0.3493 DO k = 1, klev+13494 DO j = 1, knon3495 i = ni(j)3496 !jyg<3497 !! tke(i,k,nsrf) = ytke(j,k)3498 !! tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)*ypct(j)3499 tke_x(i,k,nsrf) = ytke(j,k)3500 tke_x(i,k,is_ave) = tke_x(i,k,is_ave) + ytke(j,k)*ypct(j)3501 eps_x(i,k,nsrf) = yeps(j,k)3502 eps_x(i,k,is_ave) = eps_x(i,k,is_ave) + yeps(j,k)*ypct(j)3503 !>jyg3504 ENDDO3505 ENDDO3506 3507 ELSE ! (iflag_split .eq.0)3508 DO k = 1, klev+13509 DO j = 1, knon3510 i = ni(j)3511 wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k)3512 !jyg<3513 !! tke(i,k,nsrf) = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)3514 !! tke(i,k,is_ave) = tke(i,k,is_ave) + tke(i,k,nsrf)*ypct(j)3515 tke_x(i,k,nsrf) = ytke_x(j,k)3516 tke_x(i,k,is_ave) = tke_x(i,k,is_ave) + tke_x(i,k,nsrf)*ypct(j)3517 eps_x(i,k,nsrf) = yeps_x(j,k)3518 eps_x(i,k,is_ave) = eps_x(i,k,is_ave) + eps_x(i,k,nsrf)*ypct(j)3519 wake_dltke(i,k,is_ave) = wake_dltke(i,k,is_ave) + wake_dltke(i,k,nsrf)*ypct(j)3520 3521 3522 !>jyg3523 ENDDO3524 ENDDO3525 ENDIF ! (iflag_split .eq.0)3526 !!!3527 DO k = 2, klev3528 DO j = 1, knon3529 i = ni(j)3530 zcoefh(i,k,nsrf) = ycoefh(j,k)3531 zcoefm(i,k,nsrf) = ycoefm(j,k)3532 zcoefh(i,k,is_ave) = zcoefh(i,k,is_ave) + ycoefh(j,k)*ypct(j)3533 zcoefm(i,k,is_ave) = zcoefm(i,k,is_ave) + ycoefm(j,k)*ypct(j)3534 ENDDO3535 ENDDO3536 3537 ! print*,'Dans pbl OK3'3538 3539 IF ( nsrf .EQ. is_ter ) THEN3540 DO j = 1, knon3541 i = ni(j)3542 qsol(i) = yqsol(j)3543 #ifdef ISO3544 runoff_diag(i)=yrunoff_diag(j)3545 DO ixt=1,niso3546 xtsol(ixt,i) = yxtsol(ixt,j)3547 xtrunoff_diag(ixt,i)=yxtrunoff_diag(ixt,j)3548 ENDDO3549 #endif3550 ENDDO3551 ENDIF3552 3553 !jyg<3554 !! ftsoil(:,:,nsrf) = 0.3555 !>jyg3556 DO k = 1, nsoilmx3557 DO j = 1, knon3558 i = ni(j)3559 ftsoil(i, k, nsrf) = ytsoil(j,k)3560 ENDDO3561 ENDDO3562 3563 #ifdef ISO3564 #ifdef ISOVERIF3565 !write(*,*) 'pbl_surface 2858'3566 DO i = 1, klon3567 DO ixt=1,niso3568 call iso_verif_noNaN(xtsol(ixt,i),'pbl_surface 1405')3569 ENDDO3570 ENDDO3571 #endif3572 #ifdef ISOVERIF3573 IF (iso_eau.gt.0) THEN3574 call iso_verif_egalite_vect2D( &3575 y_d_xt,y_d_q, &3576 'pbl_surface_mod 1261',ntraciso,klon,klev)3577 ENDIF !if (iso_eau.gt.0) then3578 #endif3579 #endif3580 !!! jyg le 07/02/20123581 IF (iflag_split .ge.1) THEN3582 !!!3583 !!! nrlmd+jyg le 02/05/2011 et le 20/02/20123584 DO k = 1, klev3585 DO j = 1, knon3586 i = ni(j)3587 d_t_diss_x(i,k) = d_t_diss_x(i,k) + y_d_t_diss_x(j,k)3588 d_t_x(i,k) = d_t_x(i,k) + y_d_t_x(j,k)3589 d_q_x(i,k) = d_q_x(i,k) + y_d_q_x(j,k)3590 d_u_x(i,k) = d_u_x(i,k) + y_d_u_x(j,k)3591 d_v_x(i,k) = d_v_x(i,k) + y_d_v_x(j,k)3592 !3593 d_t_diss_w(i,k) = d_t_diss_w(i,k) + y_d_t_diss_w(j,k)3594 d_t_w(i,k) = d_t_w(i,k) + y_d_t_w(j,k)3595 d_q_w(i,k) = d_q_w(i,k) + y_d_q_w(j,k)3596 d_u_w(i,k) = d_u_w(i,k) + y_d_u_w(j,k)3597 d_v_w(i,k) = d_v_w(i,k) + y_d_v_w(j,k)3598 #ifdef ISO3599 DO ixt=1,ntraciso3600 d_xt_x(ixt,i,k) = d_xt_x(ixt,i,k) + y_d_xt_x(ixt,j,k)3601 d_xt_w(ixt,i,k) = d_xt_w(ixt,i,k) + y_d_xt_w(ixt,j,k)3602 ENDDO ! DO ixt=1,ntraciso3603 #endif3604 3605 !3606 !! d_wake_dlt(i,k) = d_wake_dlt(i,k) + y_d_t_w(i,k)-y_d_t_x(i,k)3607 !! d_wake_dlq(i,k) = d_wake_dlq(i,k) + y_d_q_w(i,k)-y_d_q_x(i,k)3608 ENDDO3609 ENDDO3610 !!!3611 ENDIF ! (iflag_split .ge.1)3612 !!!3613 3614 DO k = 1, klev3615 DO j = 1, knon3616 i = ni(j)3617 d_t_diss(i,k) = d_t_diss(i,k) + y_d_t_diss(j,k)3618 d_t(i,k) = d_t(i,k) + y_d_t(j,k)3619 d_q(i,k) = d_q(i,k) + y_d_q(j,k)3620 #ifdef ISO3621 DO ixt=1,ntraciso3622 d_xt(ixt,i,k) = d_xt(ixt,i,k) + y_d_xt(ixt,j,k)3623 ENDDO !DO ixt=1,ntraciso3624 #endif3625 d_u(i,k) = d_u(i,k) + y_d_u(j,k)3626 d_v(i,k) = d_v(i,k) + y_d_v(j,k)3627 ENDDO3628 ENDDO3629 3630 3631 IF (ok_bs) THEN3632 DO k = 1, klev3633 DO j = 1, knon3634 i = ni(j)3635 d_qbs(i,k) = d_qbs(i,k) + y_d_qbs(j,k)3636 ENDDO3637 ENDDO3638 ENDIF3639 3640 #ifdef ISO3641 #ifdef ISOVERIF3642 ! write(*,*) 'd_q,d_xt(iso_eau,554,19)=',d_q(554,19),d_xt(iso_eau,554,19)3643 ! write(*,*) 'pbl_surface 2929: d_q,d_xt(iso_eau,2,1)=',d_q(2,1),d_xt(iso_eau,2,1)3644 ! write(*,*) 'y_d_q,y_d_xt(iso_eau,2,1)=',y_d_q(2,1),y_d_xt(iso_eau,2,1)3645 ! write(*,*) 'iso_eau.gt.0=',iso_eau.gt.03646 call iso_verif_noNaN_vect2D( &3647 & d_xt, &3648 & 'pbl_surface 1385',ntraciso,klon,klev)3649 IF (iso_eau >= 0) THEN3650 call iso_verif_egalite_vect2D( &3651 y_d_xt,y_d_q, &3652 'pbl_surface_mod 2945',ntraciso,klon,klev)3653 call iso_verif_egalite_vect2D( &3654 d_xt,d_q, &3655 'pbl_surface_mod 1276',ntraciso,klon,klev)3656 ENDIF !IF (iso_eau >= 0) THEN3657 #endif3658 #endif3659 3660 ! print*,'Dans pbl OK4'3661 3662 IF (prt_level >=10) THEN3663 print *, 'pbl_surface tendencies for w: d_t_w, d_t_x, d_t ', &3664 d_t_w(1:knon,1), d_t_x(1:knon,1), d_t(1:knon,1)3665 ENDIF3666 3667 if (nsrf == is_oce .and. activate_ocean_skin >= 1) then3668 delta_sal = missing_val3669 ds_ns = missing_val3670 dt_ns = missing_val3671 delta_sst = missing_val3672 dter = missing_val3673 dser = missing_val3674 tkt = missing_val3675 tks = missing_val3676 taur = missing_val3677 sss = missing_val3678 3679 delta_sal(ni(:knon)) = ydelta_sal(:knon)3680 ds_ns(ni(:knon)) = yds_ns(:knon)3681 dt_ns(ni(:knon)) = ydt_ns(:knon)3682 delta_sst(ni(:knon)) = ydelta_sst(:knon)3683 dter(ni(:knon)) = ydter(:knon)3684 dser(ni(:knon)) = ydser(:knon)3685 tkt(ni(:knon)) = ytkt(:knon)3686 tks(ni(:knon)) = ytks(:knon)3687 taur(ni(:knon)) = ytaur(:knon)3688 sss(ni(:knon)) = ysss(:knon)3689 3690 if (activate_ocean_skin == 2 .and. type_ocean == "couple") then3691 dt_ds = missing_val3692 dt_ds(ni(:knon)) = ydt_ds(:knon)3693 end if3694 end if3695 3696 3697 !****************************************************************************************3698 ! 14) Calculate the temperature and relative humidity at 2m and the wind at 10m3699 ! Call HBTM3700 !3701 !****************************************************************************************3702 !!!3703 !3704 #undef T2m3705 #define T2m3706 #ifdef T2m3707 ! Calculations of diagnostic t,q at 2m and u, v at 10m3708 3709 ! print*,'Dans pbl OK41'3710 ! print*,'tair1,yt(:,1),y_d_t(:,1)'3711 ! print*, tair1,yt(:,1),y_d_t(:,1)3712 !!! jyg le 07/02/20123713 IF (iflag_split .eq.0) THEN3714 DO j=1, knon3715 uzon(j) = yu(j,1) + y_d_u(j,1)3716 vmer(j) = yv(j,1) + y_d_v(j,1)3717 tair1(j) = yt(j,1) + y_d_t(j,1) + y_d_t_diss(j,1)3718 qair1(j) = yq(j,1) + y_d_q(j,1)3719 zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &3720 * (ypaprs(j,1)-ypplay(j,1))3721 tairsol(j) = yts(j) + y_d_ts(j)3722 qairsol(j) = yqsurf(j)3723 ENDDO3724 ELSE ! (iflag_split .eq.0)3725 DO j=1, knon3726 uzon_x(j) = yu_x(j,1) + y_d_u_x(j,1)3727 vmer_x(j) = yv_x(j,1) + y_d_v_x(j,1)3728 tair1_x(j) = yt_x(j,1) + y_d_t_x(j,1) + y_d_t_diss_x(j,1)3729 qair1_x(j) = yq_x(j,1) + y_d_q_x(j,1)3730 zgeo1_x(j) = RD * tair1_x(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &3731 * (ypaprs(j,1)-ypplay(j,1))3732 tairsol(j) = yts(j) + y_d_ts(j)3733 !! tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf(j)3734 tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf_new(j)3735 qairsol(j) = yqsurf(j)3736 ENDDO3737 DO j=1, knon3738 uzon_w(j) = yu_w(j,1) + y_d_u_w(j,1)3739 vmer_w(j) = yv_w(j,1) + y_d_v_w(j,1)3740 tair1_w(j) = yt_w(j,1) + y_d_t_w(j,1) + y_d_t_diss_w(j,1)3741 qair1_w(j) = yq_w(j,1) + y_d_q_w(j,1)3742 zgeo1_w(j) = RD * tair1_w(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &3743 * (ypaprs(j,1)-ypplay(j,1))3744 tairsol_w(j) = tairsol(j) + (1.- ywake_s(j))*y_delta_tsurf(j)3745 qairsol(j) = yqsurf(j)3746 ENDDO3747 !!!3748 ENDIF ! (iflag_split .eq.0)3749 !!!3750 DO j=1, knon3751 ! i = ni(j)3752 ! yz0h_oupas(j) = yz0m(j)3753 ! IF(nsrf.EQ.is_oce) THEN3754 ! yz0h_oupas(j) = z0m(i,nsrf)3755 ! ENDIF3756 psfce(j)=ypaprs(j,1)3757 patm(j)=ypplay(j,1)3758 ENDDO3759 3760 IF (iflag_pbl_surface_t2m_bug==1) THEN3761 yz0h_oupas(1:knon)=yz0m(1:knon)3762 ELSE3763 yz0h_oupas(1:knon)=yz0h(1:knon)3764 ENDIF3765 3766 ! print*,'Dans pbl OK42A'3767 ! print*,'tair1,yt(:,1),y_d_t(:,1)'3768 ! print*, tair1,yt(:,1),y_d_t(:,1)3769 3770 ! Calculate the temperature and relative humidity at 2m and the wind at 10m3771 !!! jyg le 07/02/20123772 IF (iflag_split .eq.0) THEN3773 IF (iflag_new_t2mq2m==1) THEN3774 CALL stdlevvarn(klon, knon, nsrf, zxli, &3775 uzon, vmer, tair1, qair1, zgeo1, &3776 tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &3777 yt2m, yq2m, yt10m, yq10m, yu10m, yustar, &3778 yn2mout(:, nsrf, :))3779 ELSE3780 CALL stdlevvar(klon, knon, nsrf, zxli, &3781 uzon, vmer, tair1, qair1, zgeo1, &3782 tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &3783 yt2m, yq2m, yt10m, yq10m, yu10m, yustar, ypblh, rain_f, zxtsol)3784 ENDIF3785 ELSE !(iflag_split .eq.0)3786 IF (iflag_new_t2mq2m==1) THEN3787 CALL stdlevvarn(klon, knon, nsrf, zxli, &3788 uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &3789 tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, &3790 yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x, &3791 yn2mout_x(:, nsrf, :))3792 CALL stdlevvarn(klon, knon, nsrf, zxli, &3793 uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &3794 tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &3795 yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w, &3796 yn2mout_w(:, nsrf, :))3797 ELSE3798 CALL stdlevvar(klon, knon, nsrf, zxli, &3799 uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &3800 tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, &3801 yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x, ypblh_x, rain_f, zxtsol)3802 CALL stdlevvar(klon, knon, nsrf, zxli, &3803 uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &3804 tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &3805 yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w, ypblh_w, rain_f, zxtsol)3806 ENDIF3807 !!!3808 ENDIF ! (iflag_split .eq.0)3809 !!!3810 !!! jyg le 07/02/20123811 IF (iflag_split .eq.0) THEN3812 DO j=1, knon3813 i = ni(j)3814 t2m(i,nsrf)=yt2m(j)3815 q2m(i,nsrf)=yq2m(j)3816 ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman3817 ustar(i,nsrf)=yustar(j)3818 u10m(i,nsrf)=(yu10m(j) * uzon(j))/max(SQRT(uzon(j)**2+vmer(j)**2), smallestreal)3819 v10m(i,nsrf)=(yu10m(j) * vmer(j))/max(SQRT(uzon(j)**2+vmer(j)**2), smallestreal)3820 !3821 DO k = 1, 63822 n2mout(i,nsrf,k) = yn2mout(j,nsrf,k)3823 END DO3824 !3825 ENDDO3826 ELSE !(iflag_split .eq.0)3827 DO j=1, knon3828 i = ni(j)3829 t2m_x(i,nsrf)=yt2m_x(j)3830 q2m_x(i,nsrf)=yq2m_x(j)3831 ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman3832 ustar_x(i,nsrf)=yustar_x(j)3833 u10m_x(i,nsrf)=(yu10m_x(j) * uzon_x(j))/max(SQRT(uzon_x(j)**2+vmer_x(j)**2), smallestreal)3834 v10m_x(i,nsrf)=(yu10m_x(j) * vmer_x(j))/max(SQRT(uzon_x(j)**2+vmer_x(j)**2), smallestreal)3835 !3836 DO k = 1, 63837 n2mout_x(i,nsrf,k) = yn2mout_x(j,nsrf,k)3838 END DO3839 !3840 ENDDO3841 DO j=1, knon3842 i = ni(j)3843 t2m_w(i,nsrf)=yt2m_w(j)3844 q2m_w(i,nsrf)=yq2m_w(j)3845 ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman3846 ustar_w(i,nsrf)=yustar_w(j)3847 u10m_w(i,nsrf)=(yu10m_w(j) * uzon_w(j))/max(SQRT(uzon_w(j)**2+vmer_w(j)**2), smallestreal)3848 v10m_w(i,nsrf)=(yu10m_w(j) * vmer_w(j))/max(SQRT(uzon_w(j)**2+vmer_w(j)**2), smallestreal)3849 !3850 ustar(i,nsrf) = ustar_x(i,nsrf) + wake_s(i)*(ustar_w(i,nsrf)-ustar_x(i,nsrf))3851 u10m(i,nsrf) = u10m_x(i,nsrf) + wake_s(i)*(u10m_w(i,nsrf)-u10m_x(i,nsrf))3852 v10m(i,nsrf) = v10m_x(i,nsrf) + wake_s(i)*(v10m_w(i,nsrf)-v10m_x(i,nsrf))3853 !3854 DO k = 1, 63855 n2mout_w(i,nsrf,k) = yn2mout_w(j,nsrf,k)3856 END DO3857 !3858 ENDDO3859 !!!3860 ENDIF ! (iflag_split .eq.0)3861 !!!3862 3863 ! print*,'Dans pbl OK43'3864 !IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique3865 !IM Ajoute dependance type surface3866 IF (thermcep) THEN3867 !!! jyg le 07/02/20123868 IF (iflag_split .eq.0) THEN3869 DO j = 1, knon3870 i=ni(j)3871 zdelta1 = MAX(0.,SIGN(1., rtt-yt2m(j) ))3872 zx_qs1 = r2es * FOEEW(yt2m(j),zdelta1)/paprs(i,1)3873 zx_qs1 = MIN(0.5,zx_qs1)3874 zcor1 = 1./(1.-RETV*zx_qs1)3875 zx_qs1 = zx_qs1*zcor13876 3877 rh2m(i) = rh2m(i) + yq2m(j)/zx_qs1 * pctsrf(i,nsrf)3878 qsat2m(i) = qsat2m(i) + zx_qs1 * pctsrf(i,nsrf)3879 ENDDO3880 ELSE ! (iflag_split .eq.0)3881 DO j = 1, knon3882 i=ni(j)3883 zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_x(j) ))3884 zx_qs1 = r2es * FOEEW(yt2m_x(j),zdelta1)/paprs(i,1)3885 zx_qs1 = MIN(0.5,zx_qs1)3886 zcor1 = 1./(1.-RETV*zx_qs1)3887 zx_qs1 = zx_qs1*zcor13888 3889 rh2m_x(i) = rh2m_x(i) + yq2m_x(j)/zx_qs1 * pctsrf(i,nsrf)3890 qsat2m_x(i) = qsat2m_x(i) + zx_qs1 * pctsrf(i,nsrf)3891 ENDDO3892 DO j = 1, knon3893 i=ni(j)3894 zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_w(j) ))3895 zx_qs1 = r2es * FOEEW(yt2m_w(j),zdelta1)/paprs(i,1)3896 zx_qs1 = MIN(0.5,zx_qs1)3897 zcor1 = 1./(1.-RETV*zx_qs1)3898 zx_qs1 = zx_qs1*zcor13899 3900 rh2m_w(i) = rh2m_w(i) + yq2m_w(j)/zx_qs1 * pctsrf(i,nsrf)3901 qsat2m_w(i) = qsat2m_w(i) + zx_qs1 * pctsrf(i,nsrf)3902 ENDDO3903 !!!3904 ENDIF ! (iflag_split .eq.0)3905 !!!3906 ENDIF3907 !3908 IF (prt_level >=10) THEN3909 print *, 'T2m, q2m, RH2m ', &3910 t2m(1:knon,:), q2m(1:knon,:), rh2m(1:knon)3911 ENDIF3912 3913 ! print*,'OK pbl 5'3914 !3915 !!! jyg le 07/02/20123916 IF (iflag_split .eq.0) THEN3917 CALL hbtm(knon, ypaprs, ypplay, &3918 yt2m,yt10m,yq2m,yq10m,yustar,ywstar, &3919 y_flux_t,y_flux_q,yu,yv,yt,yq, &3920 ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, &3921 ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)3922 IF (prt_level >=10) THEN3923 print *,' Arg. de HBTM: yt2m ',yt2m(1:knon)3924 print *,' Arg. de HBTM: yt10m ',yt10m(1:knon)3925 print *,' Arg. de HBTM: yq2m ',yq2m(1:knon)3926 print *,' Arg. de HBTM: yq10m ',yq10m(1:knon)3927 print *,' Arg. de HBTM: yustar ',yustar(1:knon)3928 print *,' Arg. de HBTM: y_flux_t ',y_flux_t(1:knon,:)3929 print *,' Arg. de HBTM: y_flux_q ',y_flux_q(1:knon,:)3930 print *,' Arg. de HBTM: yu ',yu(1:knon,:)3931 print *,' Arg. de HBTM: yv ',yv(1:knon,:)3932 print *,' Arg. de HBTM: yt ',yt(1:knon,:)3933 print *,' Arg. de HBTM: yq ',yq(1:knon,:)3934 ENDIF3935 ELSE ! (iflag_split .eq.0)3936 CALL HBTM(knon, ypaprs, ypplay, &3937 yt2m_x,yt10m_x,yq2m_x,yq10m_x,yustar_x,ywstar_x, &3938 y_flux_t_x,y_flux_q_x,yu_x,yv_x,yt_x,yq_x, &3939 ypblh_x,ycapCL_x,yoliqCL_x,ycteiCL_x,ypblT_x, &3940 ytherm_x,ytrmb1_x,ytrmb2_x,ytrmb3_x,ylcl_x)3941 IF (prt_level >=10) THEN3942 print *,' Arg. de HBTM: yt2m_x ',yt2m_x(1:knon)3943 print *,' Arg. de HBTM: yt10m_x ',yt10m_x(1:knon)3944 print *,' Arg. de HBTM: yq2m_x ',yq2m_x(1:knon)3945 print *,' Arg. de HBTM: yq10m_x ',yq10m_x(1:knon)3946 print *,' Arg. de HBTM: yustar_x ',yustar_x(1:knon)3947 print *,' Arg. de HBTM: y_flux_t_x ',y_flux_t_x(1:knon,:)3948 print *,' Arg. de HBTM: y_flux_q_x ',y_flux_q_x(1:knon,:)3949 print *,' Arg. de HBTM: yu_x ',yu_x(1:knon,:)3950 print *,' Arg. de HBTM: yv_x ',yv_x(1:knon,:)3951 print *,' Arg. de HBTM: yt_x ',yt_x(1:knon,:)3952 print *,' Arg. de HBTM: yq_x ',yq_x(1:knon,:)3953 ENDIF3954 CALL HBTM(knon, ypaprs, ypplay, &3955 yt2m_w,yt10m_w,yq2m_w,yq10m_w,yustar_w,ywstar_w, &3956 y_flux_t_w,y_flux_q_w,yu_w,yv_w,yt_w,yq_w, &3957 ypblh_w,ycapCL_w,yoliqCL_w,ycteiCL_w,ypblT_w, &3958 ytherm_w,ytrmb1_w,ytrmb2_w,ytrmb3_w,ylcl_w)3959 !!!3960 ENDIF ! (iflag_split .eq.0)3961 !!!3962 3963 !!! jyg le 07/02/20123964 IF (iflag_split .eq.0) THEN3965 !!!3966 DO j=1, knon3967 i = ni(j)3968 pblh(i,nsrf) = ypblh(j)3969 wstar(i,nsrf) = ywstar(j)3970 plcl(i,nsrf) = ylcl(j)3971 capCL(i,nsrf) = ycapCL(j)3972 oliqCL(i,nsrf) = yoliqCL(j)3973 cteiCL(i,nsrf) = ycteiCL(j)3974 pblT(i,nsrf) = ypblT(j)3975 therm(i,nsrf) = ytherm(j)3976 trmb1(i,nsrf) = ytrmb1(j)3977 trmb2(i,nsrf) = ytrmb2(j)3978 trmb3(i,nsrf) = ytrmb3(j)3979 ENDDO3980 IF (prt_level >=10) THEN3981 print *, 'After HBTM: pblh ', pblh(1:knon,:)3982 print *, 'After HBTM: plcl ', plcl(1:knon,:)3983 print *, 'After HBTM: cteiCL ', cteiCL(1:knon,:)3984 ENDIF3985 ELSE !(iflag_split .eq.0)3986 DO j=1, knon3987 i = ni(j)3988 pblh_x(i,nsrf) = ypblh_x(j)3989 wstar_x(i,nsrf) = ywstar_x(j)3990 plcl_x(i,nsrf) = ylcl_x(j)3991 capCL_x(i,nsrf) = ycapCL_x(j)3992 oliqCL_x(i,nsrf) = yoliqCL_x(j)3993 cteiCL_x(i,nsrf) = ycteiCL_x(j)3994 pblT_x(i,nsrf) = ypblT_x(j)3995 therm_x(i,nsrf) = ytherm_x(j)3996 trmb1_x(i,nsrf) = ytrmb1_x(j)3997 trmb2_x(i,nsrf) = ytrmb2_x(j)3998 trmb3_x(i,nsrf) = ytrmb3_x(j)3999 ENDDO4000 IF (prt_level >=10) THEN4001 print *, 'After HBTM: pblh_x ', pblh_x(1:knon,:)4002 print *, 'After HBTM: plcl_x ', plcl_x(1:knon,:)4003 print *, 'After HBTM: cteiCL_x ', cteiCL_x(1:knon,:)4004 ENDIF4005 DO j=1, knon4006 i = ni(j)4007 pblh_w(i,nsrf) = ypblh_w(j)4008 wstar_w(i,nsrf) = ywstar_w(j)4009 plcl_w(i,nsrf) = ylcl_w(j)4010 capCL_w(i,nsrf) = ycapCL_w(j)4011 oliqCL_w(i,nsrf) = yoliqCL_w(j)4012 cteiCL_w(i,nsrf) = ycteiCL_w(j)4013 pblT_w(i,nsrf) = ypblT_w(j)4014 therm_w(i,nsrf) = ytherm_w(j)4015 trmb1_w(i,nsrf) = ytrmb1_w(j)4016 trmb2_w(i,nsrf) = ytrmb2_w(j)4017 trmb3_w(i,nsrf) = ytrmb3_w(j)4018 ENDDO4019 IF (prt_level >=10) THEN4020 print *, 'After HBTM: pblh_w ', pblh_w(1:knon,:)4021 print *, 'After HBTM: plcl_w ', plcl_w(1:knon,:)4022 print *, 'After HBTM: cteiCL_w ', cteiCL_w(1:knon,:)4023 ENDIF4024 !!!4025 ENDIF ! (iflag_split .eq.0)4026 !!!4027 4028 ! print*,'OK pbl 6'4029 #else4030 ! T2m not defined4031 ! No calculation4032 PRINT*,' Warning !!! No T2m calculation. Output is set to zero.'4033 #endif4034 4035 !****************************************************************************************4036 ! 15) End of loop over different surfaces4037 !4038 !****************************************************************************************4039 ENDDO loop_nbsrf4040 !4041 !----------------------------------------------------------------------------------------4042 ! Reset iflag_split4043 !4044 iflag_split=iflag_split_ref4045 4046 #ifdef ISO4047 #ifdef ISOVERIF4048 ! write(*,*) 'pbl_surface tmp 3249: d_q,d_xt(iso_eau,2,1)=',d_q(2,1),d_xt(iso_eau,2,1)4049 IF (iso_eau >= 0) THEN4050 call iso_verif_egalite_vect2D( &4051 d_xt,d_q, &4052 'pbl_surface_mod 1276',ntraciso,klon,klev)4053 ENDIF !IF (iso_eau >= 0) THEN4054 #endif4055 #endif4056 4057 !****************************************************************************************4058 ! 16) Calculate the mean value over all sub-surfaces for some variables4059 !4060 !****************************************************************************************4061 4062 z0m(:,nbsrf+1) = 0.04063 z0h(:,nbsrf+1) = 0.04064 DO nsrf = 1, nbsrf4065 DO i = 1, klon4066 z0m(i,nbsrf+1) = z0m(i,nbsrf+1) + z0m(i,nsrf)*pctsrf(i,nsrf)4067 z0h(i,nbsrf+1) = z0h(i,nbsrf+1) + z0h(i,nsrf)*pctsrf(i,nsrf)4068 ENDDO4069 ENDDO4070 4071 ! print*,'OK pbl 7'4072 zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.04073 zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.04074 zxfluxt_x(:,:) = 0.0 ; zxfluxq_x(:,:) = 0.04075 zxfluxu_x(:,:) = 0.0 ; zxfluxv_x(:,:) = 0.04076 zxfluxt_w(:,:) = 0.0 ; zxfluxq_w(:,:) = 0.04077 zxfluxu_w(:,:) = 0.0 ; zxfluxv_w(:,:) = 0.04078 #ifdef ISO4079 zxfluxxt(:,:,:) = 0.04080 zxfluxxt_x(:,:,:) = 0.04081 zxfluxxt_w(:,:,:) = 0.04082 #endif4083 4084 4085 !!! jyg le 07/02/20124086 IF (iflag_split .ge.1) THEN4087 !!!4088 !!! nrlmd & jyg les 02/05/2011, 05/02/20124089 4090 DO nsrf = 1, nbsrf4091 DO k = 1, klev4092 DO i = 1, klon4093 zxfluxt_x(i,k) = zxfluxt_x(i,k) + flux_t_x(i,k,nsrf) * pctsrf(i,nsrf)4094 zxfluxq_x(i,k) = zxfluxq_x(i,k) + flux_q_x(i,k,nsrf) * pctsrf(i,nsrf)4095 zxfluxu_x(i,k) = zxfluxu_x(i,k) + flux_u_x(i,k,nsrf) * pctsrf(i,nsrf)4096 zxfluxv_x(i,k) = zxfluxv_x(i,k) + flux_v_x(i,k,nsrf) * pctsrf(i,nsrf)4097 !4098 zxfluxt_w(i,k) = zxfluxt_w(i,k) + flux_t_w(i,k,nsrf) * pctsrf(i,nsrf)4099 zxfluxq_w(i,k) = zxfluxq_w(i,k) + flux_q_w(i,k,nsrf) * pctsrf(i,nsrf)4100 zxfluxu_w(i,k) = zxfluxu_w(i,k) + flux_u_w(i,k,nsrf) * pctsrf(i,nsrf)4101 zxfluxv_w(i,k) = zxfluxv_w(i,k) + flux_v_w(i,k,nsrf) * pctsrf(i,nsrf)4102 #ifdef ISO4103 DO ixt=1,ntraciso4104 zxfluxxt_x(ixt,i,k) = zxfluxxt_x(ixt,i,k) + flux_xt_x(ixt,i,k,nsrf) * pctsrf(i,nsrf)4105 zxfluxxt_w(ixt,i,k) = zxfluxxt_w(ixt,i,k) + flux_xt_w(ixt,i,k,nsrf) * pctsrf(i,nsrf)4106 ENDDO ! DO ixt=1,ntraciso4107 #endif4108 ENDDO4109 ENDDO4110 ENDDO4111 4112 DO i = 1, klon4113 zxsens_x(i) = - zxfluxt_x(i,1)4114 zxsens_w(i) = - zxfluxt_w(i,1)4115 ENDDO4116 !!!4117 ENDIF ! (iflag_split .ge.1)4118 !!!4119 4120 DO nsrf = 1, nbsrf4121 DO k = 1, klev4122 DO i = 1, klon4123 zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf(i,nsrf)4124 zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf(i,nsrf)4125 zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf)4126 zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf)4127 #ifdef ISO4128 DO ixt=1,niso4129 zxfluxxt(ixt,i,k) = zxfluxxt(ixt,i,k) + flux_xt(ixt,i,k,nsrf) * pctsrf(i,nsrf)4130 ENDDO ! DO ixt=1,niso4131 #endif4132 ENDDO4133 ENDDO4134 ENDDO4135 4136 DO i = 1, klon4137 zxsens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol4138 zxevap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol4139 fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)4140 ENDDO4141 4142 ! if blowing snow4143 if (ok_bs) then4144 DO nsrf = 1, nbsrf4145 DO k = 1, klev4146 DO i = 1, klon4147 zxfluxqbs(i,k) = zxfluxqbs(i,k) + flux_qbs(i,k,nsrf) * pctsrf(i,nsrf)4148 ENDDO4149 ENDDO4150 ENDDO4151 4152 DO i = 1, klon4153 zxsnowerosion(i) = zxfluxqbs(i,1) ! blowings snow flux at the surface4154 END DO4155 endif4156 4157 #ifdef ISO4158 DO i = 1, klon4159 DO ixt=1,ntraciso4160 zxxtevap(ixt,i) = - zxfluxxt(ixt,i,1)4161 ENDDO4162 ENDDO4163 #endif4164 4165 !!!4166 4167 !4168 ! Incrementer la temperature du sol4169 !4170 zxtsol(:) = 0.0 ; zxfluxlat(:) = 0.04171 zt2m(:) = 0.0 ; zq2m(:) = 0.0 ; zn2mout(:,:) = 04172 zustar(:)=0.0 ; zu10m(:) = 0.0 ; zv10m(:) = 0.04173 s_pblh(:) = 0.0 ; s_plcl(:) = 0.04174 !!! jyg le 07/02/20124175 s_pblh_x(:) = 0.0 ; s_plcl_x(:) = 0.04176 s_pblh_w(:) = 0.0 ; s_plcl_w(:) = 0.04177 !!!4178 s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.04179 s_cteiCL(:) = 0.0; s_pblT(:) = 0.04180 s_therm(:) = 0.0 ; s_trmb1(:) = 0.04181 s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.04182 wstar(:,is_ave)=0.4183 4184 ! print*,'OK pbl 9'4185 4186 !!! nrlmd le 02/05/20114187 zxfluxlat_x(:) = 0.0 ; zxfluxlat_w(:) = 0.04188 !!!4189 4190 DO nsrf = 1, nbsrf4191 DO i = 1, klon4192 ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf)4193 4194 wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &4195 + flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)4196 4197 wfevap(i,nsrf) = evap(i,nsrf)*pctsrf(i,nsrf)4198 4199 zxtsol(i) = zxtsol(i) + ts(i,nsrf) * pctsrf(i,nsrf)4200 zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf(i,nsrf)4201 ENDDO4202 ENDDO4203 !4204 !<al1 order 2 correction to zxtsol, for radiation computations (main atm effect of Ts)4205 IF (iflag_order2_sollw == 1) THEN4206 meansqT(:) = 0. ! as working buffer4207 DO nsrf = 1, nbsrf4208 DO i = 1, klon4209 meansqT(i) = meansqT(i)+(ts(i,nsrf)-zxtsol(i))**2 *pctsrf(i,nsrf)4210 ENDDO4211 ENDDO4212 zxtsol(:) = zxtsol(:)+1.5*meansqT(:)/zxtsol(:)4213 ENDIF ! iflag_order2_sollw == 14214 !>al14215 4216 !!! jyg le 07/02/20124217 IF (iflag_split .eq.0) THEN4218 DO nsrf = 1, nbsrf4219 DO i = 1, klon4220 zt2m(i) = zt2m(i) + t2m(i,nsrf) * pctsrf(i,nsrf)4221 zq2m(i) = zq2m(i) + q2m(i,nsrf) * pctsrf(i,nsrf)4222 !4223 DO k = 1, 64224 zn2mout(i,k) = zn2mout(i,k) + n2mout(i,nsrf,k) * pctsrf(i,nsrf)4225 ENDDO4226 !4227 zustar(i) = zustar(i) + ustar(i,nsrf) * pctsrf(i,nsrf)4228 wstar(i,is_ave)=wstar(i,is_ave)+wstar(i,nsrf)*pctsrf(i,nsrf)4229 zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf(i,nsrf)4230 zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf(i,nsrf)4231 4232 s_pblh(i) = s_pblh(i) + pblh(i,nsrf) * pctsrf(i,nsrf)4233 s_plcl(i) = s_plcl(i) + plcl(i,nsrf) * pctsrf(i,nsrf)4234 s_capCL(i) = s_capCL(i) + capCL(i,nsrf) * pctsrf(i,nsrf)4235 s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf(i,nsrf)4236 s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf(i,nsrf)4237 s_pblT(i) = s_pblT(i) + pblT(i,nsrf) * pctsrf(i,nsrf)4238 s_therm(i) = s_therm(i) + therm(i,nsrf) * pctsrf(i,nsrf)4239 s_trmb1(i) = s_trmb1(i) + trmb1(i,nsrf) * pctsrf(i,nsrf)4240 s_trmb2(i) = s_trmb2(i) + trmb2(i,nsrf) * pctsrf(i,nsrf)4241 s_trmb3(i) = s_trmb3(i) + trmb3(i,nsrf) * pctsrf(i,nsrf)4242 ENDDO4243 ENDDO4244 ELSE !(iflag_split .eq.0)4245 DO nsrf = 1, nbsrf4246 DO i = 1, klon4247 !!! nrlmd le 02/05/20114248 zxfluxlat_x(i) = zxfluxlat_x(i) + fluxlat_x(i,nsrf) * pctsrf(i,nsrf)4249 zxfluxlat_w(i) = zxfluxlat_w(i) + fluxlat_w(i,nsrf) * pctsrf(i,nsrf)4250 !!!4251 !!! jyg le 08/02/20124252 !! Pour le moment, on sort les valeurs dans (x) et (w) de pblh et de plcl ;4253 !! pour zt2m, on fait la moyenne surfacique sur les sous-surfaces ;4254 !! pour qsat2m, on fait la moyenne surfacique sur (x) et (w) ;4255 !! pour les autres variables, on sort les valeurs de la region (x).4256 zt2m(i) = zt2m(i) + (t2m_x(i,nsrf)+wake_s(i)*(t2m_w(i,nsrf)-t2m_x(i,nsrf))) * pctsrf(i,nsrf)4257 zq2m(i) = zq2m(i) + q2m_x(i,nsrf) * pctsrf(i,nsrf)4258 !4259 DO k = 1, 64260 zn2mout(i,k) = zn2mout(i,k) + n2mout_x(i,nsrf,k) * pctsrf(i,nsrf)4261 ENDDO4262 !4263 zustar(i) = zustar(i) + ustar_x(i,nsrf) * pctsrf(i,nsrf)4264 wstar(i,is_ave)=wstar(i,is_ave)+wstar_x(i,nsrf)*pctsrf(i,nsrf)4265 zu10m(i) = zu10m(i) + u10m_x(i,nsrf) * pctsrf(i,nsrf)4266 zv10m(i) = zv10m(i) + v10m_x(i,nsrf) * pctsrf(i,nsrf)4267 !4268 s_pblh(i) = s_pblh(i) + pblh_x(i,nsrf) * pctsrf(i,nsrf)4269 s_pblh_x(i) = s_pblh_x(i) + pblh_x(i,nsrf) * pctsrf(i,nsrf)4270 s_pblh_w(i) = s_pblh_w(i) + pblh_w(i,nsrf) * pctsrf(i,nsrf)4271 !4272 s_plcl(i) = s_plcl(i) + plcl_x(i,nsrf) * pctsrf(i,nsrf)4273 s_plcl_x(i) = s_plcl_x(i) + plcl_x(i,nsrf) * pctsrf(i,nsrf)4274 s_plcl_w(i) = s_plcl_w(i) + plcl_w(i,nsrf) * pctsrf(i,nsrf)4275 !4276 s_capCL(i) = s_capCL(i) + capCL_x(i,nsrf) * pctsrf(i,nsrf)4277 s_oliqCL(i) = s_oliqCL(i) + oliqCL_x(i,nsrf)* pctsrf(i,nsrf)4278 s_cteiCL(i) = s_cteiCL(i) + cteiCL_x(i,nsrf)* pctsrf(i,nsrf)4279 s_pblT(i) = s_pblT(i) + pblT_x(i,nsrf) * pctsrf(i,nsrf)4280 s_therm(i) = s_therm(i) + therm_x(i,nsrf) * pctsrf(i,nsrf)4281 s_trmb1(i) = s_trmb1(i) + trmb1_x(i,nsrf) * pctsrf(i,nsrf)4282 s_trmb2(i) = s_trmb2(i) + trmb2_x(i,nsrf) * pctsrf(i,nsrf)4283 s_trmb3(i) = s_trmb3(i) + trmb3_x(i,nsrf) * pctsrf(i,nsrf)4284 ENDDO4285 ENDDO4286 DO i = 1, klon4287 qsat2m(i)= qsat2m_x(i)+ wake_s(i)*(qsat2m_x(i)-qsat2m_w(i))4288 ENDDO4289 !!!4290 ENDIF ! (iflag_split .eq.0)4291 !!!4292 4293 IF (check) THEN4294 amn=MIN(ts(1,is_ter),1000.)4295 amx=MAX(ts(1,is_ter),-1000.)4296 DO i=2, klon4297 amn=MIN(ts(i,is_ter),amn)4298 amx=MAX(ts(i,is_ter),amx)4299 ENDDO4300 PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx4301 ENDIF4302 4303 !jg ?4304 !!$!4305 !!$! If a sub-surface does not exsist for a grid point, the mean value for all4306 !!$! sub-surfaces is distributed.4307 !!$!4308 !!$ DO nsrf = 1, nbsrf4309 !!$ DO i = 1, klon4310 !!$ IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN4311 !!$ ts(i,nsrf) = zxtsol(i)4312 !!$ t2m(i,nsrf) = zt2m(i)4313 !!$ q2m(i,nsrf) = zq2m(i)4314 !!$ u10m(i,nsrf) = zu10m(i)4315 !!$ v10m(i,nsrf) = zv10m(i)4316 !!$4317 !!$! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour4318 !!$ pblh(i,nsrf) = s_pblh(i)4319 !!$ plcl(i,nsrf) = s_plcl(i)4320 !!$ capCL(i,nsrf) = s_capCL(i)4321 !!$ oliqCL(i,nsrf) = s_oliqCL(i)4322 !!$ cteiCL(i,nsrf) = s_cteiCL(i)4323 !!$ pblT(i,nsrf) = s_pblT(i)4324 !!$ therm(i,nsrf) = s_therm(i)4325 !!$ trmb1(i,nsrf) = s_trmb1(i)4326 !!$ trmb2(i,nsrf) = s_trmb2(i)4327 !!$ trmb3(i,nsrf) = s_trmb3(i)4328 !!$ ENDIF4329 !!$ ENDDO4330 !!$ ENDDO4331 4332 4333 DO i = 1, klon4334 fder(i) = - 4.0*RSIGMA*zxtsol(i)**34335 ENDDO4336 4337 zxqsurf(:) = 0.04338 zxsnow(:) = 0.04339 #ifdef ISO4340 zxxtsnow(:,:) = 0.04341 #endif4342 4343 DO nsrf = 1, nbsrf4344 DO i = 1, klon4345 zxqsurf(i) = zxqsurf(i) + MAX(qsurf(i,nsrf),0.0) * pctsrf(i,nsrf)4346 zxsnow(i) = zxsnow(i) + snow(i,nsrf) * pctsrf(i,nsrf)4347 #ifdef ISO4348 DO ixt=1,niso4349 zxxtsnow(ixt,i) = zxxtsnow(ixt,i) + xtsnow(ixt,i,nsrf) * pctsrf(i,nsrf)4350 ENDDO ! DO ixt=1,niso4351 #endif4352 ENDDO4353 ENDDO4354 4355 ! Premier niveau de vent sortie dans physiq.F4356 zu1(:) = u(:,1)4357 zv1(:) = v(:,1)4358 4359 END SUBROUTINE pbl_surface4360 4361 #endif4362 4363 4364 4365 4366 326 4367 327 !
Note: See TracChangeset
for help on using the changeset viewer.
