Index: LMDZ6/branches/PBLSURF_GPUPORT/libf/phylmd/pbl_surface_mod.F90 =================================================================== --- LMDZ6/branches/PBLSURF_GPUPORT/libf/phylmd/pbl_surface_mod.F90 (revision 5868) +++ LMDZ6/branches/PBLSURF_GPUPORT/libf/phylmd/pbl_surface_mod.F90 (revision 5869) @@ -324,4044 +324,4 @@ END SUBROUTINE pbl_surface_init_iso -#endif - - -#ifndef toSupress - -! -!**************************************************************************************** -! - - SUBROUTINE pbl_surface( & - dtime, date0, itap, jour, & - debut, lafin, & - rlon, rlat, rugoro, rmu0, & - !GG lwdown_m, cldt, & - lwdown_m, pphi, cldt, & - !GG - rain_f, snow_f, bs_f, solsw_m, solswfdiff_m, sollw_m, & - gustiness, & - t, q, qbs, u, v, & -!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012 -!! t_x, q_x, t_w, q_w, & - wake_dlt, wake_dlq, & - wake_cstar, wake_s, & -!!! - pplay, paprs, pctsrf, & - ts,SFRWL, alb_dir, alb_dif,ustar, u10m, v10m,wstar, & - cdragh, cdragm, zu1, zv1, & -!jyg< (26/09/2019) - beta, & -!>jyg - alb_dir_m, alb_dif_m, zxsens, zxevap, zxsnowerosion, & - icesub_lic, alb3_lic, runoff, snowhgt, qsnow, to_ice, sissnow, & - zxtsol, zxfluxlat, zt2m, qsat2m, zn2mout, & - d_t, d_q, d_qbs, d_u, d_v, d_t_diss, & -!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012 - d_t_w, d_q_w, & - d_t_x, d_q_x, & -!! d_wake_dlt,d_wake_dlq, & - zxsens_x, zxfluxlat_x,zxsens_w,zxfluxlat_w, & -!!! -!!! nrlmd le 13/06/2011 - delta_tsurf,wake_dens,cdragh_x,cdragh_w, & - cdragm_x,cdragm_w,kh,kh_x,kh_w, & -!!! - zcoefh, zcoefm, slab_wfbils, & - qsol, zq2m, s_pblh, s_plcl, & -!!! -!!! jyg le 08/02/2012 - s_pblh_x, s_plcl_x, s_pblh_w, s_plcl_w, & -!!! - s_capCL, s_oliqCL, s_cteiCL, s_pblT, & - s_therm, s_trmb1, s_trmb2, s_trmb3, & - zustar,zu10m, zv10m, fder_print, & - zxqsurf, delta_qsurf, & - rh2m, zxfluxu, zxfluxv, & - z0m, z0h, agesno, sollw, solsw, & - d_ts, evap, fluxlat, t2m, & - wfbils, wfevap, & - flux_t, flux_u, flux_v, & - dflux_t, dflux_q, zxsnow, & -!jyg< -!! zxfluxt, zxfluxq, q2m, flux_q, tke, & - zxfluxt, zxfluxq, zxfluxqbs, q2m, flux_q, flux_qbs, tke_x, eps_x, & -!>jyg -!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012 -!! tke_x, tke_w & - wake_dltke, & -!GG treedrg & - treedrg,hice ,tice, bilg_cumul, & - fcds, fcdi, dh_basal_growth, dh_basal_melt, & - dh_top_melt, dh_snow2sic, & - dtice_melt, dtice_snow2sic , & -!GG -!FC -!AM heterogeneous continental sub-surfaces - tsurf_tersrf, tsoil_tersrf, qsurf_tersrf, tsurf_new_tersrf, & - cdragm_tersrf, cdragh_tersrf, & - swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf & -!!! -#ifdef ISO - & ,xtrain_f, xtsnow_f,xt, & - & wake_dlxt,zxxtevap,xtevap, & - & d_xt,d_xt_w,d_xt_x, & - & xtsol,dflux_xt,zxxtsnow,zxfluxxt,flux_xt, & - & h1_diag,runoff_diag,xtrunoff_diag & -#endif - & ) -!**************************************************************************************** -! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818 -! Objet: interface de "couche limite" (diffusion verticale) -! -!AA REM: -!AA----- -!AA Tout ce qui a trait au traceurs est dans phytrac maintenant -!AA pour l'instant le calcul de la couche limite pour les traceurs -!AA se fait avec cltrac et ne tient pas compte de la differentiation -!AA des sous-fraction de sol. -!AA REM bis : -!AA---------- -!AA Pour pouvoir extraire les coefficient d'echanges et le vent -!AA dans la premiere couche, 3 champs supplementaires ont ete crees -!AA zcoefh, zu1 et zv1. Pour l'instant nous avons moyenne les valeurs -!AA de ces trois champs sur les 4 subsurfaces du modele. Dans l'avenir -!AA si les informations des subsurfaces doivent etre prises en compte -!AA il faudra sortir ces memes champs en leur ajoutant une dimension, -!AA c'est a dire nbsrf (nbre de subsurface). -! -! Arguments: -! -! dtime----input-R- interval du temps (secondes) -! itap-----input-I- numero du pas de temps -! date0----input-R- jour initial -! t--------input-R- temperature (K) -! q--------input-R- vapeur d'eau (kg/kg) -! u--------input-R- vitesse u -! v--------input-R- vitesse v -! wake_dlt-input-R- temperatre difference between (w) and (x) (K) -! wake_dlq-input-R- humidity difference between (w) and (x) (kg/kg) -!wake_cstar-input-R- wake gust front speed (m/s) -! wake_s---input-R- wake fractionnal area -! ts-------input-R- temperature du sol (en Kelvin) -! paprs----input-R- pression a intercouche (Pa) -! pplay----input-R- pression au milieu de couche (Pa) -! rlat-----input-R- latitude en degree -! z0m, z0h ----input-R- longeur de rugosite (en m) -! Martin -! cldt-----input-R- total cloud fraction -! Martin -!GG -! pphi-----input-R- geopotentiel de chaque couche (g z) (reference sol) -!GG -! -! d_t------output-R- le changement pour "t" -! d_q------output-R- le changement pour "q" -! d_u------output-R- le changement pour "u" -! d_v------output-R- le changement pour "v" -! d_ts-----output-R- le changement pour "ts" -! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2) -! (orientation positive vers le bas) -! tke_x---input/output-R- tke in the (x) region (kg/m**2/s) -! wake_dltke-input/output-R- tke difference between (w) and (x) (kg/m**2/s) -! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s) -! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal -! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal -! dflux_t--output-R- derive du flux sensible -! dflux_q--output-R- derive du flux latent -! zu1------output-R- le vent dans la premiere couche -! zv1------output-R- le vent dans la premiere couche -! trmb1----output-R- deep_cape -! trmb2----output-R- inhibition -! trmb3----output-R- Point Omega -! cteiCL---output-R- Critere d'instab d'entrainmt des nuages de CL -! plcl-----output-R- Niveau de condensation -! pblh-----output-R- HCL -! pblT-----output-R- T au nveau HCL -! treedrg--output-R- tree drag (m) -! qsurf_tersrf--output-R- surface specific humidity of continental sub-surfaces -! cdragm_tersrf--output-R- momentum drag coefficient of continental sub-surfaces -! cdragh_tersrf--output-R- heat drag coefficient of continental sub-surfaces -! tsurf_new_tersrf--output-R- surface temperature of continental sub-surfaces -! swnet_tersrf--output-R- net shortwave radiation of continental sub-surfaces -! lwnet_tersrf--output-R- net longwave radiation of continental sub-surfaces -! fluxsens_tersrf--output-R- sensible heat flux of continental sub-surfaces -! fluxlat_tersrf--output-R- latent heat flux of continental sub-surfaces - - USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, carbon_cycle_tr, level_coupling_esm - USE carbon_cycle_mod, ONLY : co2_send, nbcf_out, fields_out, yfields_out, cfname_out - use hbtm_mod, only: hbtm - USE indice_sol_mod - USE time_phylmdz_mod, ONLY : day_ini,annee_ref,itau_phy - USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat, grid1dto2d_glo - USE print_control_mod, ONLY : prt_level,lunout -#ifdef ISO - USE isotopes_mod, ONLY: Rdefault,iso_eau -#ifdef ISOVERIF - USE isotopes_verif_mod -#endif -#ifdef ISOTRAC - USE isotrac_mod, only: index_iso -#endif -#endif -USE dimpft_mod_h - USE flux_arp_mod_h - USE compbl_mod_h - USE yoethf_mod_h - USE clesphys_mod_h - USE ioipsl_getin_p_mod, ONLY : getin_p - use phys_state_var_mod, only: ds_ns, dt_ns, delta_sst, delta_sal, dter, & - dser, dt_ds, zsig, zmea, & - frac_tersrf, z0m_tersrf, ratio_z0m_z0h_tersrf, albedo_tersrf !AM - use phys_output_var_mod, only: tkt, tks, taur, sss - use lmdz_blowing_snow_ini, only : zeta_bs - use wxios_mod, ONLY: missing_val_xios => missing_val, using_xios - USE netcdf, only: missing_val_netcdf => nf90_fill_real - USE dimsoil_mod_h, ONLY: nsoilmx - USE surf_param_mod, ONLY: eff_surf_param !AM - - USE yomcst_mod_h -IMPLICIT NONE - - INCLUDE "FCTTRE.h" -!FC - -!**************************************************************************************** - REAL, INTENT(IN) :: dtime ! time interval (s) - REAL, INTENT(IN) :: date0 ! initial day - INTEGER, INTENT(IN) :: itap ! time step - INTEGER, INTENT(IN) :: jour ! current day of the year - LOGICAL, INTENT(IN) :: debut ! true if first run step - LOGICAL, INTENT(IN) :: lafin ! true if last run step - REAL, DIMENSION(klon), INTENT(IN) :: rlon ! longitudes in degrees - REAL, DIMENSION(klon), INTENT(IN) :: rlat ! latitudes in degrees - REAL, DIMENSION(klon), INTENT(IN) :: rugoro ! rugosity length - REAL, DIMENSION(klon), INTENT(IN) :: rmu0 ! cosine of solar zenith angle - REAL, DIMENSION(klon), INTENT(IN) :: rain_f ! rain fall - REAL, DIMENSION(klon), INTENT(IN) :: snow_f ! snow fall - REAL, DIMENSION(klon), INTENT(IN) :: bs_f ! blowing snow fall - REAL, DIMENSION(klon), INTENT(IN) :: solsw_m ! net shortwave radiation at mean surface - REAL, DIMENSION(klon), INTENT(IN) :: solswfdiff_m ! diffuse fraction fordownward shortwave radiation at mean surface - REAL, DIMENSION(klon), INTENT(IN) :: sollw_m ! net longwave radiation at mean surface - REAL, DIMENSION(klon,klev), INTENT(IN) :: t ! temperature (K) - REAL, DIMENSION(klon,klev), INTENT(IN) :: q ! water vapour (kg/kg) - REAL, DIMENSION(klon,klev), INTENT(IN) :: qbs ! blowing snow specific content (kg/kg) - REAL, DIMENSION(klon,klev), INTENT(IN) :: u ! u speed - REAL, DIMENSION(klon,klev), INTENT(IN) :: v ! v speed - REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! mid-layer pression (Pa) - REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! pression between layers (Pa) - REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf ! sub-surface fraction -! Martin - REAL, DIMENSION(klon), INTENT(IN) :: lwdown_m ! downward longwave radiation at mean s - REAL, DIMENSION(klon), INTENT(IN) :: gustiness ! gustiness - -!GG - REAL, DIMENSION(klon,klev), INTENT(IN) :: pphi ! geopotential (m2/s2) -!GG - REAL, DIMENSION(klon), INTENT(IN) :: cldt ! total cloud - -#ifdef ISO - REAL, DIMENSION(ntraciso,klon,klev), INTENT(IN) :: xt ! water vapour (kg/kg) - REAL, DIMENSION(ntraciso,klon), INTENT(IN) :: xtrain_f ! rain fall - REAL, DIMENSION(ntraciso,klon), INTENT(IN) :: xtsnow_f ! snow fall -#endif - -!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012 -!! REAL, DIMENSION(klon,klev), INTENT(IN) :: t_x ! Temp\'erature hors poche froide -!! REAL, DIMENSION(klon,klev), INTENT(IN) :: t_w ! Temp\'erature dans la poches froide -!! REAL, DIMENSION(klon,klev), INTENT(IN) :: q_x ! -!! REAL, DIMENSION(klon,klev), INTENT(IN) :: q_w ! Pareil pour l'humidit\'e - REAL, DIMENSION(klon,klev), INTENT(IN) :: wake_dlt !temperature difference between (w) and (x) (K) - REAL, DIMENSION(klon,klev), INTENT(IN) :: wake_dlq !humidity difference between (w) and (x) (K) - REAL, DIMENSION(klon), INTENT(IN) :: wake_s ! Fraction de poches froides - REAL, DIMENSION(klon), INTENT(IN) :: wake_cstar! Vitesse d'expansion des poches froides - REAL, DIMENSION(klon), INTENT(IN) :: wake_dens -!!! -#ifdef ISO - REAL, DIMENSION(ntraciso,klon,klev), INTENT(IN) :: wake_dlxt -#endif -! Input/Output variables -!**************************************************************************************** -!jyg< - REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: beta ! Aridity factor -!>jyg - REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: ts ! temperature at surface (K) - REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: delta_tsurf !surface temperature difference between - !wake and off-wake regions -!albedo SB >>> - REAL, DIMENSIOn(6),intent(in) :: SFRWL - REAL, DIMENSION(klon, nsw, nbsrf), INTENT(INOUT) :: alb_dir,alb_dif -!albedo SB <<< -!jyg Pourquoi ustar et wstar sont-elles INOUT ? - REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: ustar ! u* (m/s) - REAL, DIMENSION(klon, nbsrf+1), INTENT(INOUT) :: wstar ! w* (m/s) - REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: u10m ! u speed at 10m - REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: v10m ! v speed at 10m -!jyg< -!! REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke - REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke_x -!>jyg - -!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012 - REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: wake_dltke ! TKE_w - TKE_x -!!! - -! Output variables -!**************************************************************************************** - REAL, DIMENSION(klon,klev+1,nbsrf+1), INTENT(OUT) :: eps_x ! TKE dissipation rate - - REAL, DIMENSION(klon), INTENT(OUT) :: cdragh ! drag coefficient for T and Q - REAL, DIMENSION(klon), INTENT(OUT) :: cdragm ! drag coefficient for wind - REAL, DIMENSION(klon), INTENT(OUT) :: zu1 ! u wind speed in first layer - REAL, DIMENSION(klon), INTENT(OUT) :: zv1 ! v wind speed in first layer -!albedo SB >>> - REAL, DIMENSION(klon, nsw), INTENT(OUT) :: alb_dir_m,alb_dif_m -!albedo SB <<< - ! Martin - REAL, DIMENSION(klon), INTENT(OUT) :: alb3_lic - ! Martin - REAL, DIMENSION(klon), INTENT(OUT) :: zxsens ! sensible heat flux at surface with inversed sign - ! (=> positive sign upwards) - REAL, DIMENSION(klon), INTENT(OUT) :: zxevap ! water vapour flux at surface, positiv upwards - REAL, DIMENSION(klon), INTENT(OUT) :: zxsnowerosion ! blowing snow flux at surface - REAL, DIMENSION(klon), INTENT(OUT) :: icesub_lic ! ice (no snow!) sublimation over ice sheet - REAL, DIMENSION(klon), INTENT(OUT) :: zxtsol ! temperature at surface, mean for each grid point -!!! jyg le ??? - REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t_w ! ! - REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q_w ! ! Tendances dans les poches - REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t_x ! ! - REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q_x ! ! Tendances hors des poches -!!! jyg - REAL, DIMENSION(klon), INTENT(OUT) :: zxfluxlat ! latent flux, mean for each grid point - REAL, DIMENSION(klon), INTENT(OUT) :: zt2m ! temperature at 2m, mean for each grid point - INTEGER, DIMENSION(klon, 6), INTENT(OUT) :: zn2mout ! number of times the 2m temperature is out of the [tsol,temp] - REAL, DIMENSION(klon), INTENT(OUT) :: qsat2m - REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_t ! change in temperature - REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_t_diss ! change in temperature - REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_q ! change in water vapour - REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_u ! change in u speed - REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_v ! change in v speed - REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_qbs ! change in blowing snow specific content - - - REAL, INTENT(OUT):: zcoefh(klon, klev+1, nbsrf + 1) ! (klon, klev, nbsrf + 1) => only use (klon, klev, nbsrf+1) - ! coef for turbulent diffusion of T and Q, mean for each grid point - - REAL, INTENT(OUT):: zcoefm(klon, klev+1, nbsrf + 1) ! (klon, klev, nbsrf + 1) => only use (klon, klev, nbsrf+1) - ! coef for turbulent diffusion of U and V (?), mean for each grid point -#ifdef ISO - REAL, DIMENSION(ntraciso,klon), INTENT(OUT) :: zxxtevap ! water vapour flux at surface, positiv upwards - REAL, DIMENSION(ntraciso,klon, klev), INTENT(OUT) :: d_xt ! change in water vapour - REAL, DIMENSION(klon), INTENT(OUT) :: runoff_diag - REAL, DIMENSION(niso,klon), INTENT(OUT) :: xtrunoff_diag - REAL, DIMENSION(ntraciso,klon,klev), INTENT(OUT) :: d_xt_w - REAL, DIMENSION(ntraciso,klon,klev), INTENT(OUT) :: d_xt_x -#endif - - - -!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012 - REAL, DIMENSION(klon), INTENT(OUT) :: zxsens_x ! Flux sensible hors poche - REAL, DIMENSION(klon), INTENT(OUT) :: zxsens_w ! Flux sensible dans la poche - REAL, DIMENSION(klon), INTENT(OUT) :: zxfluxlat_x! Flux latent hors poche - REAL, DIMENSION(klon), INTENT(OUT) :: zxfluxlat_w! Flux latent dans la poche -!! REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_wake_dlt -!! REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_wake_dlq - -! Output only for diagnostics - REAL, DIMENSION(klon), INTENT(OUT) :: cdragh_x - REAL, DIMENSION(klon), INTENT(OUT) :: cdragh_w - REAL, DIMENSION(klon), INTENT(OUT) :: cdragm_x - REAL, DIMENSION(klon), INTENT(OUT) :: cdragm_w - REAL, DIMENSION(klon), INTENT(OUT) :: kh - REAL, DIMENSION(klon), INTENT(OUT) :: kh_x - REAL, DIMENSION(klon), INTENT(OUT) :: kh_w -!!! - REAL, DIMENSION(klon), INTENT(OUT) :: slab_wfbils! heat balance at surface only for slab at ocean points - REAL, DIMENSION(klon), INTENT(OUT) :: qsol ! water height in the soil (mm) - REAL, DIMENSION(klon), INTENT(OUT) :: zq2m ! water vapour at 2m, mean for each grid point - REAL, DIMENSION(klon), INTENT(OUT) :: s_pblh ! height of the planetary boundary layer(HPBL) -!!! jyg le 08/02/2012 - REAL, DIMENSION(klon), INTENT(OUT) :: s_pblh_x ! height of the PBL in the off-wake region - REAL, DIMENSION(klon), INTENT(OUT) :: s_pblh_w ! height of the PBL in the wake region -!!! - REAL, DIMENSION(klon), INTENT(OUT) :: s_plcl ! condensation level -!!! jyg le 08/02/2012 - REAL, DIMENSION(klon), INTENT(OUT) :: s_plcl_x ! condensation level in the off-wake region - REAL, DIMENSION(klon), INTENT(OUT) :: s_plcl_w ! condensation level in the wake region -!!! - REAL, DIMENSION(klon), INTENT(OUT) :: s_capCL ! CAPE of PBL - REAL, DIMENSION(klon), INTENT(OUT) :: s_oliqCL ! liquid water intergral of PBL - REAL, DIMENSION(klon), INTENT(OUT) :: s_cteiCL ! cloud top instab. crit. of PBL - REAL, DIMENSION(klon), INTENT(OUT) :: s_pblT ! temperature at PBLH - REAL, DIMENSION(klon), INTENT(OUT) :: s_therm ! thermal virtual temperature excess - REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb1 ! deep cape, mean for each grid point - REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb2 ! inhibition, mean for each grid point - REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb3 ! point Omega, mean for each grid point - REAL, DIMENSION(klon), INTENT(OUT) :: zustar ! u* - REAL, DIMENSION(klon), INTENT(OUT) :: zu10m ! u speed at 10m, mean for each grid point - REAL, DIMENSION(klon), INTENT(OUT) :: zv10m ! v speed at 10m, mean for each grid point - REAL, DIMENSION(klon), INTENT(OUT) :: fder_print ! fder for printing (=fder(i) + dflux_t(i) + dflux_q(i)) - REAL, DIMENSION(klon), INTENT(OUT) :: zxqsurf ! humidity at surface, mean for each grid point - REAL, DIMENSION(klon), INTENT(OUT) :: delta_qsurf! humidity difference at surface, mean for each grid point - REAL, DIMENSION(klon), INTENT(OUT) :: rh2m ! relative humidity at 2m - REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxu ! u wind tension, mean for each grid point - REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxv ! v wind tension, mean for each grid point - REAL, DIMENSION(klon, nbsrf+1), INTENT(INOUT) :: z0m,z0h ! rugosity length (m) - REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: agesno ! age of snow at surface - REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: solsw ! net shortwave radiation at surface - REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: sollw ! net longwave radiation at surface - REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: d_ts ! change in temperature at surface - REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: evap ! evaporation at surface - REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: fluxlat ! latent flux - REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: t2m ! temperature at 2 meter height - REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfbils ! heat balance at surface - REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfevap ! water balance (evap) at surface weighted by srf - REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_t ! sensible heat flux (CpT) J/m**2/s (W/m**2) - ! positve orientation downwards - REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_u ! u wind tension (kg m/s)/(m**2 s) or Pascal - REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_v ! v wind tension (kg m/s)/(m**2 s) or Pascal -!FC - REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: treedrg ! tree drag (m) -!AM heterogeneous continental sub-surfaces - REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: tsurf_tersrf ! surface temperature of continental sub-surfaces (K) - REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: qsurf_tersrf ! surface specific humidity of continental sub-surfaces (kg/kg) - REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: tsurf_new_tersrf ! surface temperature of continental sub-surfaces (K) - REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: cdragm_tersrf ! momentum drag coefficient of continental sub-surfaces (-) - REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: cdragh_tersrf ! heat drag coefficient of continental sub-surfaces (-) - REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: swnet_tersrf ! net shortwave radiation of continental sub-surfaces (W/m2) - REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: lwnet_tersrf ! net longwave radiation of continental sub-surfaces (W/m2) - REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: fluxsens_tersrf ! sensible heat flux of continental sub-surfaces (W/m2) - REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: fluxlat_tersrf ! latent heat flux of continental sub-surfaces (W/m2) - REAL, DIMENSION(klon, nsoilmx, nbtersrf), INTENT(INOUT) :: tsoil_tersrf ! soil temperature of continental sub-surfaces (K) -#ifdef ISO - REAL, DIMENSION(niso,klon), INTENT(OUT) :: xtsol ! water height in the soil (mm) - REAL, DIMENSION(ntraciso,klon, nbsrf) :: xtevap ! evaporation at surface - REAL, DIMENSION(klon), INTENT(OUT) :: h1_diag ! just diagnostic, not useful -#endif - - -! Output not needed - REAL, DIMENSION(klon), INTENT(OUT) :: dflux_t ! change of sensible heat flux - REAL, DIMENSION(klon), INTENT(OUT) :: dflux_q ! change of water vapour flux - REAL, DIMENSION(klon), INTENT(OUT) :: zxsnow ! snow at surface, mean for each grid point - REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxt ! sensible heat flux, mean for each grid point - REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxq ! water vapour flux, mean for each grid point - REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxqbs ! blowing snow flux, mean for each grid point - REAL, DIMENSION(klon, nbsrf),INTENT(OUT) :: q2m ! water vapour at 2 meter height - REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_q ! water vapour flux(latent flux) (kg/m**2/s) - REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_qbs ! blowind snow vertical flux (kg/m**2 - -#ifdef ISO - REAL, DIMENSION(ntraciso,klon), INTENT(OUT) :: dflux_xt ! change of water vapour flux - REAL, DIMENSION(niso,klon), INTENT(OUT) :: zxxtsnow ! snow at surface, mean for each grid point - REAL, DIMENSION(ntraciso,klon, klev), INTENT(OUT) :: zxfluxxt ! water vapour flux, mean for each grid point - REAL, DIMENSION(ntraciso,klon, klev, nbsrf), INTENT(OUT) :: flux_xt ! water vapour flux(latent flux) (kg/m**2/s) -#endif - -! Martin -! inlandsis - REAL, DIMENSION(klon), INTENT(OUT) :: qsnow ! snow water content - REAL, DIMENSION(klon), INTENT(OUT) :: snowhgt ! snow height - REAL, DIMENSION(klon), INTENT(OUT) :: to_ice ! snow passed to ice - REAL, DIMENSION(klon), INTENT(OUT) :: sissnow ! snow in snow model - REAL, DIMENSION(klon), INTENT(OUT) :: runoff ! runoff on land ice -! Martin -!GG - REAL, DIMENSION(klon), INTENT(INOUT) :: hice ! hice - REAL, DIMENSION(klon), INTENT(INOUT) :: tice ! tice - REAL, DIMENSION(klon), INTENT(INOUT) :: bilg_cumul ! flux cumulated - REAL, DIMENSION(klon), INTENT(INOUT) :: fcds - REAL, DIMENSION(klon), INTENT(INOUT) :: fcdi - REAL, DIMENSION(klon), INTENT(INOUT) :: dh_basal_growth - REAL, DIMENSION(klon), INTENT(INOUT) :: dh_basal_melt - REAL, DIMENSION(klon), INTENT(INOUT) :: dh_top_melt - REAL, DIMENSION(klon), INTENT(INOUT) :: dh_snow2sic - REAL, DIMENSION(klon), INTENT(INOUT) :: dtice_melt - REAL, DIMENSION(klon), INTENT(INOUT) :: dtice_snow2sic -!GG - -! Other local variables -!**************************************************************************************** -! >> PC - INTEGER :: ierr - INTEGER :: n -! << PC - INTEGER :: iflag_split, iflag_split_ref - INTEGER :: i, k, nsrf - INTEGER :: knon, j - INTEGER :: idayref - INTEGER , DIMENSION(klon) :: ni - REAL :: yt1_new - REAL :: zx_alf1, zx_alf2 !valeur ambiante par extrapola - REAL :: amn, amx - REAL :: f1 ! fraction de longeurs visibles parmi tout SW intervalle - REAL, DIMENSION(klon) :: r_co2_ppm ! taux CO2 atmosphere - REAL, DIMENSION(klon) :: yts, yz0m, yz0h, ypct - REAL, DIMENSION(klon) :: yz0h_old -!albedo SB >>> - REAL, DIMENSION(klon) :: yalb,yalb_vis -!albedo SB <<< - REAL, DIMENSION(klon) :: yt1, yq1, yu1, yv1, yqbs1 - REAL, DIMENSION(klon) :: yqa - REAL, DIMENSION(klon) :: ysnow, yqsurf, yagesno, yqsol - REAL, DIMENSION(klon) :: yrain_f, ysnow_f, ybs_f -#ifdef ISO - REAL, DIMENSION(ntraciso,klon) :: yxt1 - REAL, DIMENSION(niso,klon) :: yxtsnow, yxtsol - REAL, DIMENSION(ntraciso,klon) :: yxtrain_f, yxtsnow_f - REAL, DIMENSION(klon) :: yrunoff_diag - REAL, DIMENSION(niso,klon) :: yxtrunoff_diag - REAL, DIMENSION(niso,klon) :: yRland_ice -#endif - REAL, DIMENSION(klon) :: ysolsw, ysollw - REAL, DIMENSION(klon) :: yfder - REAL, DIMENSION(klon) :: yrugoro - REAL, DIMENSION(klon) :: yfluxlat - REAL, DIMENSION(klon) :: yfluxbs - REAL, DIMENSION(klon) :: y_d_ts - REAL, DIMENSION(klon) :: y_flux_t1, y_flux_q1 - REAL, DIMENSION(klon) :: y_dflux_t, y_dflux_q -#ifdef ISO - REAL, DIMENSION(ntraciso,klon) :: y_flux_xt1 - REAL, DIMENSION(ntraciso,klon) :: y_dflux_xt -#endif - REAL, DIMENSION(klon) :: y_flux_u1, y_flux_v1 - REAL, DIMENSION(klon) :: y_flux_bs, y_flux0 - REAL, DIMENSION(klon) :: yt2m, yq2m, yu10m - INTEGER, DIMENSION(klon, nbsrf, 6) :: yn2mout, yn2mout_x, yn2mout_w - INTEGER, DIMENSION(klon, nbsrf, 6) :: n2mout, n2mout_x, n2mout_w - REAL, DIMENSION(klon) :: yustar - REAL, DIMENSION(klon) :: ywstar - REAL, DIMENSION(klon) :: ywindsp - REAL, DIMENSION(klon) :: yt10m, yq10m - REAL, DIMENSION(klon) :: ypblh - REAL, DIMENSION(klon) :: ylcl - REAL, DIMENSION(klon) :: ycapCL - REAL, DIMENSION(klon) :: yoliqCL - REAL, DIMENSION(klon) :: ycteiCL - REAL, DIMENSION(klon) :: ypblT - REAL, DIMENSION(klon) :: ytherm - REAL, DIMENSION(klon) :: ytrmb1 - REAL, DIMENSION(klon) :: ytrmb2 - REAL, DIMENSION(klon) :: ytrmb3 - REAL, DIMENSION(klon) :: uzon, vmer - REAL, DIMENSION(klon) :: tair1, qair1, tairsol - REAL, DIMENSION(klon) :: psfce, patm - REAL, DIMENSION(klon) :: qairsol, zgeo1, speed, zri1, pref !speed, zri1, pref, added by Fuxing WANG, 04/03/2015 - REAL, DIMENSION(klon) :: yz0h_oupas - REAL, DIMENSION(klon) :: yfluxsens - REAL, DIMENSION(klon) :: AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0 - REAL, DIMENSION(klon) :: AcoefH, AcoefQ, BcoefH, BcoefQ -#ifdef ISO - REAL, DIMENSION(ntraciso,klon) :: AcoefXT, BcoefXT -#endif - REAL, DIMENSION(klon) :: AcoefU, AcoefV, BcoefU, BcoefV - REAL, DIMENSION(klon) :: AcoefQBS, BcoefQBS - REAL, DIMENSION(klon) :: ypsref - REAL, DIMENSION(klon) :: yevap, yevap_pot, ytsurf_new, yalb3_new, yicesub_lic -!albedo SB >>> - REAL, DIMENSION(klon,nsw) :: yalb_dir_new, yalb_dif_new -!albedo SB <<< - REAL, DIMENSION(klon) :: ztsol - REAL, DIMENSION(klon) :: meansqT ! mean square deviation of subsurface temperatures - REAL, DIMENSION(klon) :: alb_m ! mean albedo for whole SW interval - REAL, DIMENSION(klon,klev) :: y_d_t, y_d_q, y_d_t_diss, y_d_qbs - REAL, DIMENSION(klon,klev) :: y_d_u, y_d_v - REAL, DIMENSION(klon,klev) :: y_flux_t, y_flux_q, y_flux_qbs - REAL, DIMENSION(klon,klev) :: y_flux_u, y_flux_v - REAL, DIMENSION(klon,klev) :: ycoefh,ycoefm,ycoefq,ycoefqbs - REAL, DIMENSION(klon) :: ycdragh, ycdragq, ycdragm - REAL, DIMENSION(klon,klev) :: yu, yv - REAL, DIMENSION(klon,klev) :: yt, yq, yqbs -#ifdef ISO - REAL, DIMENSION(ntraciso,klon) :: yxtevap - REAL, DIMENSION(ntraciso,klon,klev) :: y_d_xt - REAL, DIMENSION(ntraciso,klon,klev) :: y_flux_xt - REAL, DIMENSION(ntraciso,klon,klev) :: yxt -#endif - REAL, DIMENSION(klon,klev) :: ypplay, ydelp - REAL, DIMENSION(klon,klev) :: delp - REAL, DIMENSION(klon,klev+1) :: ypaprs - REAL, DIMENSION(klon,klev+1) :: ytke, yeps - REAL, DIMENSION(klon,nsoilmx) :: ytsoil -!FC - REAL, DIMENSION(klon,nvm_lmdz) :: yveget - REAL, DIMENSION(klon,nvm_lmdz) :: ylai - REAL, DIMENSION(klon,nvm_lmdz) :: yheight - REAL, DIMENSION(klon,klev) :: y_d_u_frein - REAL, DIMENSION(klon,klev) :: y_d_v_frein - REAL, DIMENSION(klon,klev) :: y_treedrg -!FC - - CHARACTER(len=80) :: abort_message - CHARACTER(len=20) :: modname = 'pbl_surface' - LOGICAL, PARAMETER :: zxli=.FALSE. ! utiliser un jeu de fonctions simples - LOGICAL, PARAMETER :: check=.FALSE. - -!!! nrlmd le 02/05/2011 -!!! jyg le 07/02/2012 - REAL, DIMENSION(klon) :: ywake_s, ywake_cstar, ywake_dens -!!! - REAL, DIMENSION(klon,klev+1) :: ytke_x, ytke_w, yeps_x, yeps_w - REAL, DIMENSION(klon,klev+1) :: ywake_dltke - REAL, DIMENSION(klon,klev) :: yu_x, yv_x, yu_w, yv_w - REAL, DIMENSION(klon,klev) :: yt_x, yq_x, yt_w, yq_w - REAL, DIMENSION(klon,klev) :: ycoefh_x, ycoefm_x, ycoefh_w, ycoefm_w - REAL, DIMENSION(klon,klev) :: ycoefq_x, ycoefq_w - REAL, DIMENSION(klon) :: ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w - REAL, DIMENSION(klon) :: ycdragm_x, ycdragm_w - REAL, DIMENSION(klon) :: AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x - REAL, DIMENSION(klon) :: AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w - REAL, DIMENSION(klon) :: AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x - REAL, DIMENSION(klon) :: AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w - REAL, DIMENSION(klon) :: y_flux_t1_x, y_flux_q1_x, y_flux_t1_w, y_flux_q1_w - REAL, DIMENSION(klon) :: y_flux_u1_x, y_flux_v1_x, y_flux_u1_w, y_flux_v1_w - REAL, DIMENSION(klon,klev) :: y_flux_t_x, y_flux_q_x, y_flux_t_w, y_flux_q_w - REAL, DIMENSION(klon,klev) :: y_flux_u_x, y_flux_v_x, y_flux_u_w, y_flux_v_w - REAL, DIMENSION(klon) :: yfluxlat_x, yfluxlat_w - REAL, DIMENSION(klon,klev) :: y_d_t_x, y_d_q_x, y_d_t_w, y_d_q_w - REAL, DIMENSION(klon,klev) :: y_d_t_diss_x, y_d_t_diss_w - REAL, DIMENSION(klon,klev) :: d_t_diss_x, d_t_diss_w - REAL, DIMENSION(klon,klev) :: y_d_u_x, y_d_v_x, y_d_u_w, y_d_v_w - REAL, DIMENSION(klon, klev, nbsrf) :: flux_t_x, flux_q_x, flux_t_w, flux_q_w - REAL, DIMENSION(klon, klev, nbsrf) :: flux_u_x, flux_v_x, flux_u_w, flux_v_w - REAL, DIMENSION(klon, nbsrf) :: fluxlat_x, fluxlat_w - REAL, DIMENSION(klon, klev) :: zxfluxt_x, zxfluxq_x, zxfluxt_w, zxfluxq_w - REAL, DIMENSION(klon, klev) :: zxfluxu_x, zxfluxv_x, zxfluxu_w, zxfluxv_w - REAL :: zx_qs_surf, zcor_surf, zdelta_surf -!jyg< - REAL, DIMENSION(klon) :: ybeta - REAL, DIMENSION(klon) :: ybeta_prev -!>jyg - REAL, DIMENSION(klon, klev) :: d_u_x - REAL, DIMENSION(klon, klev) :: d_u_w - REAL, DIMENSION(klon, klev) :: d_v_x - REAL, DIMENSION(klon, klev) :: d_v_w - - REAL, DIMENSION(klon,klev) :: CcoefH, CcoefQ, DcoefH, DcoefQ - REAL, DIMENSION(klon,klev) :: CcoefU, CcoefV, DcoefU, DcoefV - REAL, DIMENSION(klon,klev) :: CcoefQBS, DcoefQBS - REAL, DIMENSION(klon,klev) :: CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x - REAL, DIMENSION(klon,klev) :: CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w - REAL, DIMENSION(klon,klev) :: CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x - REAL, DIMENSION(klon,klev) :: CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w - REAL, DIMENSION(klon,klev) :: Kcoef_hq, Kcoef_m, gama_h, gama_q - REAL, DIMENSION(klon,klev) :: gama_qbs, Kcoef_qbs - REAL, DIMENSION(klon,klev) :: Kcoef_hq_x, Kcoef_m_x, gama_h_x, gama_q_x - REAL, DIMENSION(klon,klev) :: Kcoef_hq_w, Kcoef_m_w, gama_h_w, gama_q_w - REAL, DIMENSION(klon) :: alf_1, alf_2, alf_1_x, alf_2_x, alf_1_w, alf_2_w -#ifdef ISO - REAL, DIMENSION(ntraciso,klon,klev) :: yxt_x, yxt_w - REAL, DIMENSION(ntraciso,klon) :: y_flux_xt1_x , y_flux_xt1_w - REAL, DIMENSION(ntraciso,klon,klev) :: y_flux_xt_x,y_d_xt_x,zxfluxxt_x - REAL, DIMENSION(ntraciso,klon,klev) :: y_flux_xt_w,y_d_xt_w,zxfluxxt_w - REAL, DIMENSION(ntraciso,klon,klev,nbsrf) :: flux_xt_x, flux_xt_w - REAL, DIMENSION(ntraciso,klon) :: AcoefXT_x, BcoefXT_x - REAL, DIMENSION(ntraciso,klon) :: AcoefXT_w, BcoefXT_w - REAL, DIMENSION(ntraciso,klon,klev) :: CcoefXT, DcoefXT - REAL, DIMENSION(ntraciso,klon,klev) :: CcoefXT_x, DcoefXT_x - REAL, DIMENSION(ntraciso,klon,klev) :: CcoefXT_w, DcoefXT_w - REAL, DIMENSION(ntraciso,klon,klev) :: gama_xt,gama_xt_x,gama_xt_w -#endif -!!! -!!!jyg le 08/02/2012 - REAL, DIMENSION(klon, nbsrf) :: windsp -! - REAL, DIMENSION(klon, nbsrf) :: t2m_x - REAL, DIMENSION(klon, nbsrf) :: q2m_x - REAL, DIMENSION(klon) :: rh2m_x - REAL, DIMENSION(klon) :: qsat2m_x - REAL, DIMENSION(klon, nbsrf) :: u10m_x - REAL, DIMENSION(klon, nbsrf) :: v10m_x - REAL, DIMENSION(klon, nbsrf) :: ustar_x - REAL, DIMENSION(klon, nbsrf) :: wstar_x -! - REAL, DIMENSION(klon, nbsrf) :: pblh_x - REAL, DIMENSION(klon, nbsrf) :: plcl_x - REAL, DIMENSION(klon, nbsrf) :: capCL_x - REAL, DIMENSION(klon, nbsrf) :: oliqCL_x - REAL, DIMENSION(klon, nbsrf) :: cteiCL_x - REAL, DIMENSION(klon, nbsrf) :: pblt_x - REAL, DIMENSION(klon, nbsrf) :: therm_x - REAL, DIMENSION(klon, nbsrf) :: trmb1_x - REAL, DIMENSION(klon, nbsrf) :: trmb2_x - REAL, DIMENSION(klon, nbsrf) :: trmb3_x -! - REAL, DIMENSION(klon, nbsrf) :: t2m_w - REAL, DIMENSION(klon, nbsrf) :: q2m_w - REAL, DIMENSION(klon) :: rh2m_w - REAL, DIMENSION(klon) :: qsat2m_w - REAL, DIMENSION(klon, nbsrf) :: u10m_w - REAL, DIMENSION(klon, nbsrf) :: v10m_w - REAL, DIMENSION(klon, nbsrf) :: ustar_w - REAL, DIMENSION(klon, nbsrf) :: wstar_w -! - REAL, DIMENSION(klon, nbsrf) :: pblh_w - REAL, DIMENSION(klon, nbsrf) :: plcl_w - REAL, DIMENSION(klon, nbsrf) :: capCL_w - REAL, DIMENSION(klon, nbsrf) :: oliqCL_w - REAL, DIMENSION(klon, nbsrf) :: cteiCL_w - REAL, DIMENSION(klon, nbsrf) :: pblt_w - REAL, DIMENSION(klon, nbsrf) :: therm_w - REAL, DIMENSION(klon, nbsrf) :: trmb1_w - REAL, DIMENSION(klon, nbsrf) :: trmb2_w - REAL, DIMENSION(klon, nbsrf) :: trmb3_w -! - REAL, DIMENSION(klon) :: yt2m_x - REAL, DIMENSION(klon) :: yq2m_x - REAL, DIMENSION(klon) :: yt10m_x - REAL, DIMENSION(klon) :: yq10m_x - REAL, DIMENSION(klon) :: yu10m_x - REAL, DIMENSION(klon) :: yv10m_x - REAL, DIMENSION(klon) :: yustar_x - REAL, DIMENSION(klon) :: ywstar_x -! - REAL, DIMENSION(klon) :: ypblh_x - REAL, DIMENSION(klon) :: ylcl_x - REAL, DIMENSION(klon) :: ycapCL_x - REAL, DIMENSION(klon) :: yoliqCL_x - REAL, DIMENSION(klon) :: ycteiCL_x - REAL, DIMENSION(klon) :: ypblt_x - REAL, DIMENSION(klon) :: ytherm_x - REAL, DIMENSION(klon) :: ytrmb1_x - REAL, DIMENSION(klon) :: ytrmb2_x - REAL, DIMENSION(klon) :: ytrmb3_x -! - REAL, DIMENSION(klon) :: yt2m_w - REAL, DIMENSION(klon) :: yq2m_w - REAL, DIMENSION(klon) :: yt10m_w - REAL, DIMENSION(klon) :: yq10m_w - REAL, DIMENSION(klon) :: yu10m_w - REAL, DIMENSION(klon) :: yv10m_w - REAL, DIMENSION(klon) :: yustar_w - REAL, DIMENSION(klon) :: ywstar_w -! - REAL, DIMENSION(klon) :: ypblh_w - REAL, DIMENSION(klon) :: ylcl_w - REAL, DIMENSION(klon) :: ycapCL_w - REAL, DIMENSION(klon) :: yoliqCL_w - REAL, DIMENSION(klon) :: ycteiCL_w - REAL, DIMENSION(klon) :: ypblt_w - REAL, DIMENSION(klon) :: ytherm_w - REAL, DIMENSION(klon) :: ytrmb1_w - REAL, DIMENSION(klon) :: ytrmb2_w - REAL, DIMENSION(klon) :: ytrmb3_w -! - 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/2015 - REAL, DIMENSION(klon) :: zgeo1_x, tair1_x, qair1_x, tairsol_x -! - 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/2015 - REAL, DIMENSION(klon) :: zgeo1_w, tair1_w, qair1_w, tairsol_w - REAL, DIMENSION(klon) :: yus0, yvs0 - -!!! jyg le 25/03/2013 -!! Variables intermediaires pour le raccord des deux colonnes \`a la surface -!jyg< -!! REAL :: dd_Ch -!! REAL :: dd_Cm -!! REAL :: dd_Kh -!! REAL :: dd_Km -!! REAL :: dd_u -!! REAL :: dd_v -!! REAL :: dd_t -!! REAL :: dd_q -!! REAL :: dd_AH -!! REAL :: dd_AQ -!! REAL :: dd_AU -!! REAL :: dd_AV -!! REAL :: dd_BH -!! REAL :: dd_BQ -!! REAL :: dd_BU -!! REAL :: dd_BV -!! -!! REAL :: dd_KHp -!! REAL :: dd_KQp -!! REAL :: dd_KUp -!! REAL :: dd_KVp -!>jyg - -!!! -!!! nrlmd le 13/06/2011 - REAL, DIMENSION(klon) :: y_delta_flux_t1, y_delta_flux_q1, y_delta_flux_u1, y_delta_flux_v1 - REAL, DIMENSION(klon) :: y_delta_tsurf, y_delta_tsurf_new - REAL, DIMENSION(klon) :: delta_coef, tau_eq - REAL, DIMENSION(klon) :: HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn - REAL, DIMENSION(klon) :: phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0 - REAL, DIMENSION(klon) :: y_delta_qsurf - REAL, DIMENSION(klon) :: y_delta_qsats - REAL, DIMENSION(klon) :: yg_T, yg_Q - REAL, DIMENSION(klon) :: yGamma_dTs_phiT, yGamma_dQs_phiQ - REAL, DIMENSION(klon) :: ydTs_ins, ydqs_ins -! - REAL, PARAMETER :: facteur=2./sqrt(3.14) - REAL, PARAMETER :: inertia=2000. - REAL, DIMENSION(klon) :: ydtsurf_th - REAL :: zdelta_surf_x,zdelta_surf_w,zx_qs_surf_x,zx_qs_surf_w - REAL :: zcor_surf_x,zcor_surf_w - REAL :: mod_wind_x, mod_wind_w - REAL :: rho1 - REAL, DIMENSION(klon) :: Kech_h ! Coefficient d'echange pour l'energie - REAL, DIMENSION(klon) :: Kech_h_x, Kech_h_w - REAL, DIMENSION(klon) :: Kech_m - REAL, DIMENSION(klon) :: Kech_m_x, Kech_m_w - REAL, DIMENSION(klon) :: yts_x, yts_w - REAL, DIMENSION(klon) :: yqsatsrf0_x, yqsatsrf0_w - REAL, DIMENSION(klon) :: yqsurf_x, yqsurf_w -!jyg< -!! REAL, DIMENSION(klon) :: Kech_Hp, Kech_H_xp, Kech_H_wp -!! REAL, DIMENSION(klon) :: Kech_Qp, Kech_Q_xp, Kech_Q_wp -!! REAL, DIMENSION(klon) :: Kech_Up, Kech_U_xp, Kech_U_wp -!! REAL, DIMENSION(klon) :: Kech_Vp, Kech_V_xp, Kech_V_wp -!>jyg - - REAL :: fact_cdrag - REAL :: z1lay - - REAL :: vent -! -! For debugging with IOIPSL - INTEGER, DIMENSION(nbp_lon*nbp_lat) :: ndexbg - REAL :: zjulian - REAL, DIMENSION(klon) :: tabindx - REAL, DIMENSION(nbp_lon,nbp_lat) :: zx_lon, zx_lat - REAL, DIMENSION(nbp_lon,nbp_lat) :: debugtab - - - REAL, DIMENSION(klon,nbsrf) :: pblh ! height of the planetary boundary layer - REAL, DIMENSION(klon,nbsrf) :: plcl ! condensation level - REAL, DIMENSION(klon,nbsrf) :: capCL - REAL, DIMENSION(klon,nbsrf) :: oliqCL - REAL, DIMENSION(klon,nbsrf) :: cteiCL - REAL, DIMENSION(klon,nbsrf) :: pblT - REAL, DIMENSION(klon,nbsrf) :: therm - REAL, DIMENSION(klon,nbsrf) :: trmb1 ! deep cape - REAL, DIMENSION(klon,nbsrf) :: trmb2 ! inhibition - REAL, DIMENSION(klon,nbsrf) :: trmb3 ! point Omega - REAL, DIMENSION(klon,nbsrf) :: zx_rh2m, zx_qsat2m - REAL, DIMENSION(klon,nbsrf) :: zx_t1 - REAL, DIMENSION(klon, nbsrf) :: alb ! mean albedo for whole SW interval - REAL, DIMENSION(klon,nbsrf) :: snowerosion - REAL, DIMENSION(klon) :: ylwdown ! jg : temporary (ysollwdown) - REAL, DIMENSION(klon) :: ygustiness ! jg : temporary (ysollwdown) - - REAL :: zx_qs1, zcor1, zdelta1 - - ! Martin - REAL, DIMENSION(klon, nbsrf) :: sollwd ! net longwave radiation at surface - REAL, DIMENSION(klon) :: ytoice - REAL, DIMENSION(klon) :: ysnowhgt, yqsnow, ysissnow, yrunoff - REAL, DIMENSION(klon) :: yzmea - REAL, DIMENSION(klon) :: yzsig - REAL, DIMENSION(klon) :: ycldt - REAL, DIMENSION(klon) :: yrmu0 - ! Martin - REAL, DIMENSION(klon) :: yri0 - - REAL, DIMENSION(klon):: ydelta_sst, ydelta_sal, yds_ns, ydt_ns, ydter, & - ydser, ydt_ds, ytkt, ytks, ytaur, ysss - ! compression of delta_sst, delta_sal, ds_ns, dt_ns, dter, dser, - ! dt_ds, tkt, tks, taur, sss on ocean points - REAL :: missing_val - - ! GG - REAL, DIMENSION(klon,klev) :: ytheta - REAL, DIMENSION(klon,klev) :: ypphii - REAL, DIMENSION(klon,klev) :: ypphi - REAL, DIMENSION(klon,klev) :: ydthetadz - REAL, DIMENSION(klon) :: ydthetadz300 - REAL, DIMENSION(klon) :: Ampl - ! GG - - ! AM ! - REAL, DIMENSION(klon) :: albedo_eff - REAL, DIMENSION(klon, nbtersrf) :: yfrac_tersrf, yz0m_tersrf, yz0h_tersrf -#ifdef ISO - REAL, DIMENSION(klon) :: h1 - INTEGER :: ixt -!#ifdef ISOVERIF -! integer iso_verif_positif_nostop -!#endif -#endif - -!**************************************************************************************** -! End of declarations -!**************************************************************************************** - IF (using_xios) THEN - missing_val=missing_val_xios - ELSE - missing_val=missing_val_netcdf - ENDIF - - IF (prt_level >=10) print *,' -> pbl_surface, itap ',itap -! -!!jyg iflag_split = mod(iflag_pbl_split,2) -!!jyg iflag_split = mod(iflag_pbl_split,10) -! -! Flags controlling the splitting of the turbulent boundary layer: -! iflag_split_ref = 0 ==> no splitting -! = 1 ==> splitting without coupling with surface temperature -! = 2 ==> splitting with coupling with surface temperature over land -! = 3 ==> splitting over ocean; no splitting over land -! iflag_split: actual flag controlling the splitting. -! iflag_split = iflag_split_ref outside the sub-surface loop -! = iflag_split_ref if iflag_split_ref = 0, 1, or 2 -! = 0 over land if iflga_split_ref = 3 -! = 1 over ocean if iflga_split_ref = 3 - - iflag_split_ref = mod(iflag_pbl_split,10) - iflag_split = iflag_split_ref - -#ifdef ISO -#ifdef ISOVERIF - DO i=1,klon - DO ixt=1,niso - CALL iso_verif_noNaN(xtsol(ixt,i),'pbl_surface 608') - ENDDO - ENDDO -#endif -#ifdef ISOVERIF - DO i=1,klon - IF (iso_eau >= 0) THEN - CALL iso_verif_egalite_choix(Rland_ice(iso_eau,i),1.0, & - & 'pbl_surf_mod 585',errmax,errmaxrel) - CALL iso_verif_egalite_choix(xtsnow_f(iso_eau,i),snow_f(i), & - & 'pbl_surf_mod 594',errmax,errmaxrel) - IF (iso_verif_egalite_choix_nostop(xtsol(iso_eau,i),qsol(i), & - & 'pbl_surf_mod 596',errmax,errmaxrel) == 1) THEN - WRITE(*,*) 'i=',i - STOP - ENDIF - DO nsrf=1,nbsrf - CALL iso_verif_egalite_choix(xtsnow(iso_eau,i,nsrf),snow(i,nsrf), & - & 'pbl_surf_mod 598',errmax,errmaxrel) - ENDDO - ENDIF !IF (iso_eau >= 0) THEN - ENDDO !DO i=1,knon - DO k=1,klev - DO i=1,klon - IF (iso_eau >= 0) THEN - CALL iso_verif_egalite_choix(xt(iso_eau,i,k),q(i,k), & - & 'pbl_surf_mod 595',errmax,errmaxrel) - ENDIF !IF (iso_eau >= 0) THEN - ENDDO !DO i=1,knon - ENDDO !DO k=1,klev -#endif -#endif - - - -!**************************************************************************************** -! Force soil water content to qsol0 if qsol0>0 and VEGET=F (use bucket -! instead of ORCHIDEE) - IF (qsol0>=0.) THEN - PRINT*,'WARNING : On impose qsol=',qsol0 - qsol(:)=qsol0 -#ifdef ISO - DO ixt=1,niso - xtsol(ixt,:)=qsol0*Rdefault(ixt) - ENDDO -#ifdef ISOTRAC - DO ixt=1+niso,ntraciso - xtsol(ixt,:)=qsol0*Rdefault(index_iso(ixt)) - ENDDO -#endif -#endif - ENDIF -!**************************************************************************************** - -!**************************************************************************************** -! 2) Initialization to zero -!**************************************************************************************** -! -! 2a) Initialization of all argument variables with INTENT(OUT) -!**************************************************************************************** - cdragh(:)=0. ; cdragm(:)=0. - zu1(:)=0. ; zv1(:)=0. - yus0(:)=0. ; yvs0(:)=0. -!albedo SB >>> - alb_dir_m=0. ; alb_dif_m=0. ; alb3_lic(:)=0. -!albedo SB <<< - zxsens(:)=0. ; zxevap(:)=0. ; zxtsol(:)=0. ; zxsnowerosion(:)=0. - d_t_w(:,:)=0. ; d_q_w(:,:)=0. ; d_t_x(:,:)=0. ; d_q_x(:,:)=0. - zxfluxlat(:)=0. - zt2m(:)=0. ; zq2m(:)=0. ; qsat2m(:)=0. ; rh2m(:)=0. - zn2mout(:,:)=0 ; - d_t(:,:)=0. ; d_t_diss(:,:)=0. ; d_q(:,:)=0. ; d_qbs(:,:)=0. ; d_u(:,:)=0. ; d_v(:,:)=0. - zcoefh(:,:,:)=0. ; zcoefm(:,:,:)=0. - zxsens_x(:)=0. ; zxsens_w(:)=0. ; zxfluxlat_x(:)=0. ; zxfluxlat_w(:)=0. - cdragh_x(:)=0. ; cdragh_w(:)=0. ; cdragm_x(:)=0. ; cdragm_w(:)=0. - kh(:)=0. ; kh_x(:)=0. ; kh_w(:)=0. - slab_wfbils(:)=0. - s_pblh(:)=0. ; s_pblh_x(:)=0. ; s_pblh_w(:)=0. - s_plcl(:)=0. ; s_plcl_x(:)=0. ; s_plcl_w(:)=0. - s_capCL(:)=0. ; s_oliqCL(:)=0. ; s_cteiCL(:)=0. ; s_pblT(:)=0. - s_therm(:)=0. - s_trmb1(:)=0. ; s_trmb2(:)=0. ; s_trmb3(:)=0. - zustar(:)=0. - zu10m(:)=0. ; zv10m(:)=0. - fder_print(:)=0. - zxqsurf(:)=0. - delta_qsurf(:) = 0. - zxfluxu(:,:)=0. ; zxfluxv(:,:)=0. - solsw(:,:)=0. ; sollw(:,:)=0. - d_ts(:,:)=0. - evap(:,:)=0. - snowerosion(:,:)=0. - fluxlat(:,:)=0. - wfbils(:,:)=0. ; wfevap(:,:)=0. ; - flux_t(:,:,:)=0. ; flux_q(:,:,:)=0. ; flux_u(:,:,:)=0. ; flux_v(:,:,:)=0. - flux_qbs(:,:,:)=0. - dflux_t(:)=0. ; dflux_q(:)=0. - zxsnow(:)=0. - zxfluxt(:,:)=0. ; zxfluxq(:,:)=0.; zxfluxqbs(:,:)=0. - qsnow(:)=0. ; snowhgt(:)=0. ; to_ice(:)=0. ; sissnow(:)=0. - runoff(:)=0. ; icesub_lic(:)=0. -#ifdef ISO -zxxtevap(:,:)=0. - d_xt(:,:,:)=0. - d_xt_x(:,:,:)=0. - d_xt_w(:,:,:)=0. - flux_xt(:,:,:,:)=0. -! xtsnow(:,:,:)=0.! attention, xtsnow est l'équivalent de snow et non de qsnow - xtevap(:,:,:)=0. -#endif - IF (iflag_pbl<20.or.iflag_pbl>=30) THEN - zcoefh(:,:,:) = 0.0 - zcoefh(:,1,:) = 999999. ! zcoefh(:,k=1) should never be used - zcoefm(:,:,:) = 0.0 - zcoefm(:,1,:) = 999999. ! - ELSE - zcoefm(:,:,is_ave)=0. - zcoefh(:,:,is_ave)=0. - ENDIF -!! -! The components "is_ave" of tke_x and wake_deltke are "OUT" variables -!jyg< -!! tke(:,:,is_ave)=0. - tke_x(:,:,is_ave)=0. - eps_x(:,:,is_ave)=0. - - wake_dltke(:,:,is_ave)=0. -!>jyg -!!! jyg le 23/02/2013 - t2m(:,:) = 999999. ! t2m and q2m are meaningfull only over sub-surfaces - q2m(:,:) = 999999. ! actually present in the grid cell. -!!! - rh2m(:) = 0. ; qsat2m(:) = 0. -!!! -!!! jyg le 10/02/2012 - rh2m_x(:) = 0. ; qsat2m_x(:) = 0. ; rh2m_w(:) = 0. ; qsat2m_w(:) = 0. - -! 2b) Initialization of all local variables that will be compressed later -!**************************************************************************************** -!! cdragh = 0.0 ; cdragm = 0.0 ; dflux_t = 0.0 ; dflux_q = 0.0 - ypct = 0.0 ; yts = 0.0 ; ysnow = 0.0 -!! zv1 = 0.0 ; yqsurf = 0.0 -!albedo SB >>> - yqsurf = 0.0 ; yalb = 0.0 ; yalb_vis = 0.0 -!albedo SB <<< - yrain_f = 0.0 ; ysnow_f = 0.0 ; ybs_f=0.0 ; yfder = 0.0 ; ysolsw = 0.0 - ysollw = 0.0 ; yz0m = 0.0 ; yz0h = 0.0 ; yz0h_oupas = 0.0 ; yu1 = 0.0 - yv1 = 0.0 ; ypaprs = 0.0 ; ypplay = 0.0 ; yqbs1 = 0.0 - ydelp = 0.0 ; yu = 0.0 ; yv = 0.0 ; yt = 0.0 - yq = 0.0 ; y_dflux_t = 0.0 ; y_dflux_q = 0.0 - yqbs(:,:)=0.0 - yrugoro = 0.0 ; ywindsp = 0.0 -!! d_ts = 0.0 ; yfluxlat=0.0 ; flux_t = 0.0 ; flux_q = 0.0 - yfluxlat=0.0 ; y_flux0(:)=0.0 -!! flux_u = 0.0 ; flux_v = 0.0 ; d_t = 0.0 ; d_q = 0.0 -!! d_t_diss= 0.0 ;d_u = 0.0 ; d_v = 0.0 - yqsol = 0.0 - - ytke=0. - yeps=0. - yri0(:)=0. -!FC - y_treedrg=0. - - ! Martin - ysnowhgt = 0.0; yqsnow = 0.0 ; yrunoff = 0.0 ; ytoice =0.0 - yalb3_new = 0.0 ; ysissnow = 0.0 - ycldt = 0.0 ; yrmu0 = 0.0 - ! Martin - y_d_qbs(:,:)=0.0 - -!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012 - ytke_x=0. ; ytke_w=0. ; ywake_dltke=0. - yeps_x=0. ; yeps_w=0. - y_d_t_x=0. ; y_d_t_w=0. ; y_d_q_x=0. ; y_d_q_w=0. -!! d_t_w=0. ; d_q_w=0. -!! d_t_x=0. ; d_q_x=0. -!! d_wake_dlt=0. ; d_wake_dlq=0. - yfluxlat_x=0. ; yfluxlat_w=0. - ywake_s=0. ; ywake_cstar=0. ;ywake_dens=0. -!!! -!!! nrlmd le 13/06/2011 - tau_eq=0. ; delta_coef=0. - y_delta_flux_t1=0. - ydtsurf_th=0. - yts_x(:)=0. ; yts_w(:)=0. - y_delta_tsurf(:)=0. ; y_delta_qsurf(:)=0. - yqsurf_x(:)=0. ; yqsurf_w(:)=0. - yg_T(:) = 0. ; yg_Q(:) = 0. - yGamma_dTs_phiT(:) = 0. ; yGamma_dQs_phiQ(:) = 0. - ydTs_ins(:) = 0. ; ydqs_ins(:) = 0. - -!!! - ytsoil = 999999. -!FC - y_d_u_frein(:,:)=0. - y_d_v_frein(:,:)=0. -!FC - -#ifdef ISO - yxtrain_f = 0.0 ; yxtsnow_f = 0.0 - yxtsnow = 0.0 - yxt = 0.0 - yxtsol = 0.0 - flux_xt = 0.0 - yRland_ice = 0.0 -! d_xt = 0.0 - y_dflux_xt = 0.0 - dflux_xt=0.0 - y_d_xt_x=0. ; y_d_xt_w=0. -#endif - -! >> PC -!the yfields_out variable is defined in (klon,nbcf_out) even if it is used on -!the ORCHIDEE grid and as such should be defined in yfields_out(knon,nbcf_out) but -!the knon variable is not known at that level of pbl_surface_mod - -!the yfields_in variable is defined in (klon,nbcf_in) even if it is used on the -!ORCHIDEE grid and as such should be defined in yfields_in(knon,nbcf_in) but the -!knon variable is not known at that level of pbl_surface_mod - - yfields_out(:,:) = 0. -! << PC - -!GG - ypphi = 0.0 -!GG - - -! 2c) Initialization of all local variables computed within the subsurface loop and used later on -!**************************************************************************************** - d_t_diss_x(:,:) = 0. ; d_t_diss_w(:,:) = 0. - d_u_x(:,:)=0. ; d_u_w(:,:)=0. - d_v_x(:,:)=0. ; d_v_w(:,:)=0. - flux_t_x(:,:,:)=0. ; flux_t_w(:,:,:)=0. - flux_q_x(:,:,:)=0. ; flux_q_w(:,:,:)=0. -! -!jyg< - flux_u_x(:,:,:)=0. ; flux_u_w(:,:,:)=0. - flux_v_x(:,:,:)=0. ; flux_v_w(:,:,:)=0. - fluxlat_x(:,:)=0. ; fluxlat_w(:,:)=0. -!>jyg -#ifdef ISO - flux_xt_x(:,:,:,:)=0. ; flux_xt_w(:,:,:,:)=0. -#endif -! -!jyg< -! pblh,plcl,capCL,cteiCL ... are meaningfull only over sub-surfaces -! actually present in the grid cell ==> value set to 999999. -! -!jyg< - ustar(:,:) = 999999. - wstar(:,:) = 999999. - windsp(:,:) = SQRT(u10m(:,:)**2 + v10m(:,:)**2 ) - u10m(:,:) = 999999. - v10m(:,:) = 999999. -!>jyg -! - pblh(:,:) = 999999. ! Hauteur de couche limite - plcl(:,:) = 999999. ! Niveau de condensation de la CLA - capCL(:,:) = 999999. ! CAPE de couche limite - oliqCL(:,:) = 999999. ! eau_liqu integree de couche limite - cteiCL(:,:) = 999999. ! cloud top instab. crit. couche limite - pblt(:,:) = 999999. ! T a la Hauteur de couche limite - therm(:,:) = 999999. - trmb1(:,:) = 999999. ! deep_cape - trmb2(:,:) = 999999. ! inhibition - trmb3(:,:) = 999999. ! Point Omega -! - t2m_x(:,:) = 999999. - q2m_x(:,:) = 999999. - ustar_x(:,:) = 999999. - wstar_x(:,:) = 999999. - u10m_x(:,:) = 999999. - v10m_x(:,:) = 999999. -! - pblh_x(:,:) = 999999. ! Hauteur de couche limite - plcl_x(:,:) = 999999. ! Niveau de condensation de la CLA - capCL_x(:,:) = 999999. ! CAPE de couche limite - oliqCL_x(:,:) = 999999. ! eau_liqu integree de couche limite - cteiCL_x(:,:) = 999999. ! cloud top instab. crit. couche limite - pblt_x(:,:) = 999999. ! T a la Hauteur de couche limite - therm_x(:,:) = 999999. - trmb1_x(:,:) = 999999. ! deep_cape - trmb2_x(:,:) = 999999. ! inhibition - trmb3_x(:,:) = 999999. ! Point Omega -! - t2m_w(:,:) = 999999. - q2m_w(:,:) = 999999. - ustar_w(:,:) = 999999. - wstar_w(:,:) = 999999. - u10m_w(:,:) = 999999. - v10m_w(:,:) = 999999. - - pblh_w(:,:) = 999999. ! Hauteur de couche limite - plcl_w(:,:) = 999999. ! Niveau de condensation de la CLA - capCL_w(:,:) = 999999. ! CAPE de couche limite - oliqCL_w(:,:) = 999999. ! eau_liqu integree de couche limite - cteiCL_w(:,:) = 999999. ! cloud top instab. crit. couche limite - pblt_w(:,:) = 999999. ! T a la Hauteur de couche limite - therm_w(:,:) = 999999. - trmb1_w(:,:) = 999999. ! deep_cape - trmb2_w(:,:) = 999999. ! inhibition - trmb3_w(:,:) = 999999. ! Point Omega -!!! -! -!!! -!**************************************************************************************** -! 3) - Calculate pressure thickness of each layer -! - Calculate the wind at first layer -! - Mean calculations of albedo -! - Calculate net radiance at sub-surface -!**************************************************************************************** - DO k = 1, klev - DO i = 1, klon - delp(i,k) = paprs(i,k)-paprs(i,k+1) - ENDDO - ENDDO - -!**************************************************************************************** -! Test for rugos........ from physiq.. A la fin plutot??? -! -!**************************************************************************************** - - DO nsrf = 1, nbsrf - DO i = 1, klon - z0m(i,nsrf) = MAX(z0m(i,nsrf),z0min) - z0h(i,nsrf) = MAX(z0h(i,nsrf),z0min) - ENDDO - ENDDO - - ! AM heterogeneous continental subsurfaces - ! compute time-independent effective surface parameters - IF (iflag_hetero_surf .GT. 0) THEN - CALL eff_surf_param(klon, nbtersrf, albedo_tersrf, frac_tersrf, 'ARI', albedo_eff) - ENDIF - -! Mean calculations of albedo -! -! * alb : mean albedo for whole SW interval -! -! Mean albedo for grid point -! * alb_m : mean albedo at whole SW interval - - alb_dir_m(:,:) = 0.0 - alb_dif_m(:,:) = 0.0 - DO k = 1, nsw - DO nsrf = 1, nbsrf - DO i = 1, klon - ! AM heterogeneous continental sub-surfaces - IF (nsrf .EQ. is_ter .AND. iflag_hetero_surf .GT. 0) THEN - alb_dir(i,k,nsrf) = albedo_eff(i) - alb_dif(i,k,nsrf) = albedo_eff(i) - ENDIF - ! - alb_dir_m(i,k) = alb_dir_m(i,k) + alb_dir(i,k,nsrf) * pctsrf(i,nsrf) - alb_dif_m(i,k) = alb_dif_m(i,k) + alb_dif(i,k,nsrf) * pctsrf(i,nsrf) - ENDDO - ENDDO - ENDDO - -! We here suppose the fraction f1 of incoming radiance of visible radiance -! as a fraction of all shortwave radiance - f1 = 0.5 -! f1 = 1 ! put f1=1 to recreate old calculations - -!f1 is already included with SFRWL values in each surf files - alb=0.0 - DO k=1,nsw - DO nsrf = 1, nbsrf - DO i = 1, klon - alb(i,nsrf) = alb(i,nsrf) + alb_dir(i,k,nsrf)*SFRWL(k) - ENDDO - ENDDO - ENDDO - - alb_m=0.0 - DO k = 1,nsw - DO i = 1, klon - alb_m(i) = alb_m(i) + alb_dir_m(i,k)*SFRWL(k) - ENDDO - ENDDO -!albedo SB <<< - - - -! Calculation of mean temperature at surface grid points - ztsol(:) = 0.0 - DO nsrf = 1, nbsrf - DO i = 1, klon - ztsol(i) = ztsol(i) + ts(i,nsrf)*pctsrf(i,nsrf) - ENDDO - ENDDO - -! Linear distrubution on sub-surface of long- and shortwave net radiance - DO nsrf = 1, nbsrf - DO i = 1, klon - sollw(i,nsrf) = sollw_m(i) + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ts(i,nsrf)) -!--OB this line is not satisfactory because alb is the direct albedo not total albedo - solsw(i,nsrf) = solsw_m(i) * (1.-alb(i,nsrf)) / (1.-alb_m(i)) - ENDDO - ENDDO -! -!al1 - -!--OB add diffuse fraction of SW down - DO n=1,nbcf_out - IF (cfname_out(n) == "swdownfdiff" ) fields_out(:,n) = solswfdiff_m(:) - ENDDO -! >> PC - IF (carbon_cycle_cpl .AND. carbon_cycle_tr .AND. nbcf_out.GT.0 ) THEN - r_co2_ppm(:) = co2_send(:) - DO n=1,nbcf_out - IF (cfname_out(n) == "atmco2" ) fields_out(:,n) = co2_send(:) - ENDDO - ENDIF - IF ( .NOT. carbon_cycle_tr .AND. nbcf_out.GT.0 ) THEN - r_co2_ppm(:) = co2_ppm ! Constant field - DO n=1,nbcf_out - IF (cfname_out(n) == "atmco2" ) fields_out(:,n) = co2_ppm - ENDDO - ENDIF -! << PC - -!**************************************************************************************** -! 4) Loop over different surfaces -! -! Only points containing a fraction of the sub surface will be treated. -! -!**************************************************************************************** - !<<<<<<<<<<<<< - loop_nbsrf: DO nsrf = 1, nbsrf !<<<<<<<<<<<<< - !<<<<<<<<<<<<< - IF (prt_level >=10) print *,' Loop nsrf ',nsrf -! - IF (iflag_split_ref == 3) THEN - IF (nsrf == is_oce) THEN - iflag_split = 1 - ELSE - iflag_split=0 - ENDIF !! (nsrf == is_oce) - ELSE - iflag_split = iflag_split_ref - ENDIF !! (iflag_split_ref == 3) - -! Search for index(ni) and size(knon) of domaine to treat - ni(:) = 0 - knon = 0 - DO i = 1, klon - IF (pctsrf(i,nsrf) > 0.) THEN - knon = knon + 1 - ni(knon) = i - ENDIF - ENDDO - -!!! jyg le 19/08/2012 -! IF (knon <= 0) THEN -! IF (prt_level >= 10) print *,' no grid point for nsrf= ',nsrf -! cycle loop_nbsrf -! ENDIF -!!! - - -!**************************************************************************************** -! 5) Compress variables -! -!**************************************************************************************** - -! -!jyg< (20190926) -! Provisional : set ybeta to standard values - IF (nsrf .NE. is_ter) THEN - ybeta(1:knon) = 1. - ELSE - IF (iflag_split .EQ. 0) THEN - ybeta(1:knon) = 1. - ELSE - DO j = 1, knon - i = ni(j) - ybeta(j) = beta(i,nsrf) - ENDDO - ENDIF ! (iflag_split .LE.1) - ENDIF ! (nsrf .NE. is_ter) -!>jyg -! - DO j = 1, knon - i = ni(j) - ypct(j) = pctsrf(i,nsrf) - yts(j) = ts(i,nsrf) - ysnow(j) = snow(i,nsrf) - yqsurf(j) = qsurf(i,nsrf) - yalb(j) = alb(i,nsrf) -!albedo SB >>> - yalb_vis(j) = alb_dir(i,1,nsrf) - IF (nsw==6) THEN - yalb_vis(j)=(alb_dir(i,1,nsrf)*SFRWL(1)+alb_dir(i,2,nsrf)*SFRWL(2) & - +alb_dir(i,3,nsrf)*SFRWL(3))/(SFRWL(1)+SFRWL(2)+SFRWL(3)) - ENDIF -!albedo SB <<< - yrain_f(j) = rain_f(i) - ysnow_f(j) = snow_f(i) - ybs_f(j) = bs_f(i) - yagesno(j) = agesno(i,nsrf) - yfder(j) = fder(i) - ylwdown(j) = lwdown_m(i) - ygustiness(j) = gustiness(i) - ysolsw(j) = solsw(i,nsrf) - ysollw(j) = sollw(i,nsrf) - yz0m(j) = z0m(i,nsrf) - yz0h(j) = z0h(i,nsrf) - yrugoro(j) = rugoro(i) - yu1(j) = u(i,1) - yv1(j) = v(i,1) - yqbs1(j) = qbs(i,1) - ypaprs(j,klev+1) = paprs(i,klev+1) -!jyg< -!! ywindsp(j) = SQRT(u10m(i,nsrf)**2 + v10m(i,nsrf)**2 ) - ywindsp(j) = windsp(i,nsrf) -!>jyg - ! Martin and Etienne - yzmea(j) = zmea(i) - yzsig(j) = zsig(i) - ycldt(j) = cldt(i) - yrmu0(j) = rmu0(i) - ! Martin -!!! nrlmd le 13/06/2011 - y_delta_tsurf(j)=delta_tsurf(i,nsrf) - yfluxbs(j)=0.0 - y_flux_bs(j) = 0.0 -!!! -#ifdef ISO - DO ixt=1,ntraciso - yxtrain_f(ixt,j) = xtrain_f(ixt,i) - yxtsnow_f(ixt,j) = xtsnow_f(ixt,i) - ENDDO - DO ixt=1,niso - yxtsnow(ixt,j) = xtsnow(ixt,i,nsrf) - ENDDO - !IF (nsrf == is_lic) THEN - DO ixt=1,niso - yRland_ice(ixt,j)= Rland_ice(ixt,i) - ENDDO - !endif !IF (nsrf == is_lic) THEN -#ifdef ISOVERIF - IF (iso_eau >= 0) THEN - call iso_verif_egalite_choix(ysnow_f(j), & - & yxtsnow_f(iso_eau,j),'pbl_surf_mod 862', & - & errmax,errmaxrel) - call iso_verif_egalite_choix(ysnow(j), & - & yxtsnow(iso_eau,j),'pbl_surf_mod 872', & - & errmax,errmaxrel) - ENDIF -#endif -#ifdef ISOVERIF - DO ixt=1,ntraciso - call iso_verif_noNaN(yxtsnow_f(ixt,j),'pbl_surf_mod 921') - ENDDO -#endif -#endif - ENDDO -! >> PC -!--compressing fields_out onto ORCHIDEE grid -!--these fields are shared and used directly surf_land_orchidee_mod - DO n = 1, nbcf_out - DO j = 1, knon - i = ni(j) - yfields_out(j,n) = fields_out(i,n) - ENDDO - ENDDO -! << PC - DO k = 1, klev - DO j = 1, knon - i = ni(j) - ypaprs(j,k) = paprs(i,k) - ypplay(j,k) = pplay(i,k) - ydelp(j,k) = delp(i,k) - ENDDO - ENDDO -! -!!! jyg le 07/02/2012 et le 10/04/2013 - DO k = 1, klev+1 - DO j = 1, knon - i = ni(j) -!jyg< -!! ytke(j,k) = tke(i,k,nsrf) - ytke(j,k) = tke_x(i,k,nsrf) - ENDDO - ENDDO -!>jyg - DO k = 1, klev - DO j = 1, knon - i = ni(j) - y_treedrg(j,k) = treedrg(i,k,nsrf) - yu(j,k) = u(i,k) - yv(j,k) = v(i,k) - yt(j,k) = t(i,k) - yq(j,k) = q(i,k) - yqbs(j,k)=qbs(i,k) -!! GG - ypphi(j,k) = pphi(i,k) -!! - -#ifdef ISO - DO ixt=1,ntraciso - yxt(ixt,j,k) = xt(ixt,i,k) - ENDDO !DO ixt=1,ntraciso -#endif - ENDDO - ENDDO -! - IF (iflag_split.GE.1) THEN -!!! nrlmd le 02/05/2011 - DO k = 1, klev - DO j = 1, knon - i = ni(j) - yu_x(j,k) = u(i,k) - yv_x(j,k) = v(i,k) - yt_x(j,k) = t(i,k)-wake_s(i)*wake_dlt(i,k) - yq_x(j,k) = q(i,k)-wake_s(i)*wake_dlq(i,k) - yu_w(j,k) = u(i,k) - yv_w(j,k) = v(i,k) - yt_w(j,k) = t(i,k)+(1.-wake_s(i))*wake_dlt(i,k) - yq_w(j,k) = q(i,k)+(1.-wake_s(i))*wake_dlq(i,k) -!!! -#ifdef ISO - DO ixt=1,ntraciso - yxt_x(ixt,j,k) = xt(ixt,i,k)-wake_s(i)*wake_dlxt(ixt,i,k) - yxt_w(ixt,j,k) = xt(ixt,i,k)+(1.-wake_s(i))*wake_dlxt(ixt,i,k) - ENDDO -#endif - ENDDO - ENDDO - - IF (prt_level .ge. 10) THEN - print *,'pbl_surface, wake_s(1), wake_dlt(1,:) ', wake_s(1), wake_dlt(1,:) - print *,'pbl_surface, wake_s(1), wake_dlq(1,:) ', wake_s(1), wake_dlq(1,:) - ENDIF - -!!! nrlmd le 02/05/2011 - DO k = 1, klev+1 - DO j = 1, knon - i = ni(j) -!jyg< -!! ytke_x(j,k) = tke(i,k,nsrf)-wake_s(i)*wake_dltke(i,k,nsrf) -!! ytke_w(j,k) = tke(i,k,nsrf)+(1.-wake_s(i))*wake_dltke(i,k,nsrf) -!! ywake_dltke(j,k) = wake_dltke(i,k,nsrf) -!! ytke(j,k) = tke(i,k,nsrf) -! - ytke_x(j,k) = tke_x(i,k,nsrf) - ytke(j,k) = tke_x(i,k,nsrf)+wake_s(i)*wake_dltke(i,k,nsrf) - ytke_w(j,k) = tke_x(i,k,nsrf)+wake_dltke(i,k,nsrf) - ywake_dltke(j,k) = wake_dltke(i,k,nsrf) - -!>jyg - ENDDO - ENDDO -!!! -!!! jyg le 07/02/2012 - DO j = 1, knon - i = ni(j) - ywake_s(j)=wake_s(i) - ywake_cstar(j)=wake_cstar(i) - ywake_dens(j)=wake_dens(i) - ENDDO -!!! -!!! nrlmd le 13/06/2011 - DO j=1,knon - yts_x(j)=yts(j)-ywake_s(j)*y_delta_tsurf(j) - yts_w(j)=yts(j)+(1.-ywake_s(j))*y_delta_tsurf(j) - ENDDO -!!! - ENDIF ! (iflag_split .ge.1) -!!! - DO k = 1, nsoilmx - DO j = 1, knon - i = ni(j) - ytsoil(j,k) = ftsoil(i,k,nsrf) - ENDDO - ENDDO - - ! qsol(water height in soil) only for bucket continental model - IF ( nsrf .EQ. is_ter .AND. .NOT. ok_veget ) THEN - DO j = 1, knon - i = ni(j) - yqsol(j) = qsol(i) -#ifdef ISO - DO ixt=1,niso - yxtsol(ixt,j) = xtsol(ixt,i) - ENDDO -#endif - ENDDO - ENDIF - - if (nsrf == is_oce .and. activate_ocean_skin >= 1) then - if (activate_ocean_skin == 2 .and. type_ocean == "couple") then - ydelta_sal(:knon) = delta_sal(ni(:knon)) - ydelta_sst(:knon) = delta_sst(ni(:knon)) - ydter(:knon) = dter(ni(:knon)) - ydser(:knon) = dser(ni(:knon)) - ydt_ds(:knon) = dt_ds(ni(:knon)) - end if - - yds_ns(:knon) = ds_ns(ni(:knon)) - ydt_ns(:knon) = dt_ns(ni(:knon)) - end if - -!**************************************************************************************** -! 6a) Calculate coefficients for turbulent diffusion at surface, cdragh et cdragm. -! -!**************************************************************************************** - - -!!! jyg le 07/02/2012 - IF (iflag_split .eq.0) THEN -!!! -!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012 -! Faire disparaitre les lignes commentees fin 2015 (le temps des tests) -! CALL clcdrag( knon, nsrf, ypaprs, ypplay, & -! yu(:,1), yv(:,1), yt(:,1), yq(:,1), & -! yts, yqsurf, yrugos, & -! ycdragm, ycdragh ) -! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag - DO i = 1, knon -! print*,'PBL ',i,RD -! print*,'PBL ',yt(i,1),ypaprs(i,1),ypplay(i,1) - zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) & - * (ypaprs(i,1)-ypplay(i,1)) - speed(i) = SQRT(yu(i,1)**2+yv(i,1)**2) - ENDDO - - !!! AM heterogeneous continental subsurfaces - IF (nsrf .EQ. is_ter) THEN - ! compute time-dependent effective surface parameters (function of zgeo1) !! AM - IF (iflag_hetero_surf .GT. 0) THEN - DO n=1,nbtersrf - DO j=1,knon - i = ni(j) - yfrac_tersrf(j,n) = frac_tersrf(i,n) - yz0m_tersrf(j,n) = z0m_tersrf(i,n) - IF (ratio_z0m_z0h_tersrf(i,n) .NE. 0.) THEN - yz0h_tersrf(j,n) = z0m_tersrf(i,n) / ratio_z0m_z0h_tersrf(i,n) - ELSE - yz0h_tersrf(j,n) = 0. - ENDIF - ENDDO - ENDDO - ! - CALL eff_surf_param(knon, nbtersrf, yz0m_tersrf, yfrac_tersrf, 'CDN', yz0m, zgeo1/RG) - CALL eff_surf_param(knon, nbtersrf, yz0h_tersrf, yfrac_tersrf, 'CDN', yz0h, zgeo1/RG) - ! - ENDIF - ENDIF - -! - CALL cdrag(knon, nsrf, & - speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1), s_pblh, & - yts, yqsurf, yz0m, yz0h, yri0, 0, & - ycdragm, ycdragh, zri1, pref, rain_f, zxtsol, ypplay(:,1)) - -! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013 - IF (ok_prescr_ust) THEN - DO i = 1, knon - print *,'ycdragm avant=',ycdragm(i) - vent= sqrt(yu(i,1)*yu(i,1)+yv(i,1)*yv(i,1)) -! ycdragm(i) = ust*ust/(1.+(yu(i,1)*yu(i,1)+yv(i,1)*yv(i,1))) -! ycdragm(i) = ust*ust/((1.+sqrt(yu(i,1)*yu(i,1)+yv(i,1)*yv(i,1))) & -! *sqrt(yu(i,1)*yu(i,1)+yv(i,1)*yv(i,1))) - ycdragm(i) = ust*ust/(1.+vent)/vent -! print *,'ycdragm ust yu yv apres=',ycdragm(i),ust,yu(i,1),yv(i,1) - ENDDO - ENDIF - - IF (prt_level >=10) print *,'cdrag -> ycdragh ', ycdragh(1:knon) - ELSE !(iflag_split .eq.0) - -! Faire disparaitre les lignes commentees fin 2015 (le temps des tests) -! CALL clcdrag( knon, nsrf, ypaprs, ypplay, & -! yu_x(:,1), yv_x(:,1), yt_x(:,1), yq_x(:,1), & -! yts_x, yqsurf, yrugos, & -! ycdragm_x, ycdragh_x ) -! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag - DO i = 1, knon - zgeo1_x(i) = RD * yt_x(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) & - * (ypaprs(i,1)-ypplay(i,1)) - speed_x(i) = SQRT(yu_x(i,1)**2+yv_x(i,1)**2) - ENDDO - - - CALL cdrag(knon, nsrf, & - speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),s_pblh_x,& - yts_x, yqsurf_x, yz0m, yz0h, yri0, 0, & - ycdragm_x, ycdragh_x, zri1_x, pref_x, rain_f, zxtsol, ypplay(:,1) ) - -! --- special Dice. JYG+MPL 25112013 - IF (ok_prescr_ust) THEN - DO i = 1, knon -! print *,'ycdragm_x avant=',ycdragm_x(i) - vent= sqrt(yu_x(i,1)*yu_x(i,1)+yv_x(i,1)*yv_x(i,1)) - ycdragm_x(i) = ust*ust/(1.+vent)/vent -! print *,'ycdragm_x ust yu yv apres=',ycdragm_x(i),ust,yu_x(i,1),yv_x(i,1) - ENDDO - ENDIF - IF (prt_level >=10) print *,'clcdrag -> ycdragh_x ', ycdragh_x(1:knon) -! -! Faire disparaitre les lignes commentees fin 2015 (le temps des tests) -! CALL clcdrag( knon, nsrf, ypaprs, ypplay, & -! yu_w(:,1), yv_w(:,1), yt_w(:,1), yq_w(:,1), & -! yts_w, yqsurf, yz0m, & -! ycdragm_w, ycdragh_w ) -! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag - DO i = 1, knon - zgeo1_w(i) = RD * yt_w(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) & - * (ypaprs(i,1)-ypplay(i,1)) - speed_w(i) = SQRT(yu_w(i,1)**2+yv_w(i,1)**2) - ENDDO - CALL cdrag(knon, nsrf, & - speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),s_pblh_w,& - yts_w, yqsurf_w, yz0m, yz0h, yri0, 0, & - ycdragm_w, ycdragh_w, zri1_w, pref_w, rain_f, zxtsol, ypplay(:,1) ) -! - IF(ok_bug_zg_wk_pbl) THEN - zgeo1(1:knon) = wake_s(1:knon)*zgeo1_w(1:knon) + (1.-wake_s(1:knon))*zgeo1_x(1:knon) - ELSE - zgeo1(1:knon) = ywake_s(1:knon)*zgeo1_w(1:knon) + (1.-ywake_s(1:knon))*zgeo1_x(1:knon) - ENDIF - -! --- special Dice. JYG+MPL 25112013 puis BOMEX - IF (ok_prescr_ust) THEN - DO i = 1, knon -! print *,'ycdragm_w avant=',ycdragm_w(i) - vent= sqrt(yu_w(i,1)*yu_w(i,1)+yv_w(i,1)*yv_w(i,1)) - ycdragm_w(i) = ust*ust/(1.+vent)/vent -! print *,'ycdragm_w ust yu yv apres=',ycdragm_w(i),ust,yu_w(i,1),yv_w(i,1) - ENDDO - ENDIF - IF (prt_level >=10) print *,'clcdrag -> ycdragh_w ', ycdragh_w(1:knon) -!!! - ENDIF ! (iflag_split .eq.0) -!!! - - -!**************************************************************************************** -! 6b) Calculate coefficients for turbulent diffusion in the atmosphere, ycoefh et ycoefm. -! -!**************************************************************************************** - -!!! jyg le 07/02/2012 - IF (iflag_split .eq.0) THEN -!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012 - IF (prt_level >=10) THEN - print *,' args coef_diff_turb: yu ', yu(1:knon,:) - print *,' args coef_diff_turb: yv ', yv(1:knon,:) - print *,' args coef_diff_turb: yq ', yq(1:knon,:) - print *,' args coef_diff_turb: yt ', yt(1:knon,:) - print *,' args coef_diff_turb: yts ', yts(1:knon) - print *,' args coef_diff_turb: yz0m ', yz0m(1:knon) - print *,' args coef_diff_turb: yqsurf ', yqsurf(1:knon) - print *,' args coef_diff_turb: ycdragm ', ycdragm(1:knon) - print *,' args coef_diff_turb: ycdragh ', ycdragh(1:knon) - print *,' args coef_diff_turb: ytke ', ytke(1:knon,:) - ENDIF - - IF (iflag_pbl>=50) THEN - CALL call_atke(dtime,knon,klev,nsrf,ni,ycdragm(1:knon), ycdragh(1:knon),yus0(1:knon),yvs0(1:knon),yts(1:knon), & - yu(1:knon,:),yv(1:knon,:),yt(1:knon,:),yq(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:), & - ytke(1:knon,:),yeps(1:knon,:), ycoefm(1:knon,:), ycoefh(1:knon,:)) - - ELSE - - CALL coef_diff_turb(dtime, nsrf, knon, ni, & - ypaprs, ypplay, yu, yv, yq, yt, yts, yqsurf, ycdragm, & - ycoefm, ycoefh, ytke, yeps, y_treedrg) -! ycoefm, ycoefh, ytke) -!FC y_treedrg ajoute - IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN -! In this case, coef_diff_turb is called for the Cd only - DO k = 2, klev - DO j = 1, knon - i = ni(j) - ycoefh(j,k) = zcoefh(i,k,nsrf) - ycoefm(j,k) = zcoefm(i,k,nsrf) - ENDDO - ENDDO - ENDIF - - ENDIF ! iflag_pbl >= 50 - - IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh ',ycoefh(1:knon,:) - - - ELSE !(iflag_split .eq.0) - - - IF (prt_level >=10) THEN - print *,' args coef_diff_turb: yu_x ', yu_x(1:knon,:) - print *,' args coef_diff_turb: yv_x ', yv_x(1:knon,:) - print *,' args coef_diff_turb: yq_x ', yq_x(1:knon,:) - print *,' args coef_diff_turb: yt_x ', yt_x(1:knon,:) - print *,' args coef_diff_turb: yts_x ', yts_x(1:knon) - print *,' args coef_diff_turb: yqsurf ', yqsurf(1:knon) - print *,' args coef_diff_turb: ycdragm_x ', ycdragm_x(1:knon) - print *,' args coef_diff_turb: ycdragh_x ', ycdragh_x(1:knon) - print *,' args coef_diff_turb: ytke_x ', ytke_x(1:knon,:) - ENDIF - - - IF (iflag_pbl>=50) THEN - - 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), & - yu_x(1:knon,:),yv_x(1:knon,:),yt_x(1:knon,:),yq_x(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:), & - ytke_x(1:knon,:),yeps_x(1:knon,:),ycoefm_x(1:knon,:), ycoefh_x(1:knon,:)) - - ELSE - - CALL coef_diff_turb(dtime, nsrf, knon, ni, & - ypaprs, ypplay, yu_x, yv_x, yq_x, yt_x, yts_x, yqsurf_x, ycdragm_x, & - ycoefm_x, ycoefh_x, ytke_x,yeps_x,y_treedrg) -! ycoefm_x, ycoefh_x, ytke_x) -!FC doit on le mettre ( on ne l utilise pas si il y a du spliting) - IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN -! In this case, coef_diff_turb is called for the Cd only - DO k = 2, klev - DO j = 1, knon - i = ni(j) - ycoefh_x(j,k) = zcoefh(i,k,nsrf) - ycoefm_x(j,k) = zcoefm(i,k,nsrf) - ENDDO - ENDDO - ENDIF - - ENDIF ! iflag_pbl >= 50 - - IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_x ',ycoefh_x(1:knon,:) -! - IF (prt_level >=10) THEN - print *,' args coef_diff_turb: yu_w ', yu_w(1:knon,:) - print *,' args coef_diff_turb: yv_w ', yv_w(1:knon,:) - print *,' args coef_diff_turb: yq_w ', yq_w(1:knon,:) - print *,' args coef_diff_turb: yt_w ', yt_w(1:knon,:) - print *,' args coef_diff_turb: yts_w ', yts_w(1:knon) - print *,' args coef_diff_turb: yqsurf ', yqsurf(1:knon) - print *,' args coef_diff_turb: ycdragm_w ', ycdragm_w(1:knon) - print *,' args coef_diff_turb: ycdragh_w ', ycdragh_w(1:knon) - print *,' args coef_diff_turb: ytke_w ', ytke_w(1:knon,:) - ENDIF - - IF (iflag_pbl>=50) THEN - - 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), & - yu_w(1:knon,:),yv_w(1:knon,:),yt_w(1:knon,:),yq_w(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:), & - ytke_w(1:knon,:),yeps_w(1:knon,:),ycoefm_w(1:knon,:),ycoefh_w(1:knon,:)) - - ELSE - - CALL coef_diff_turb(dtime, nsrf, knon, ni, & - ypaprs, ypplay, yu_w, yv_w, yq_w, yt_w, yts_w, yqsurf_w, ycdragm_w, & - ycoefm_w, ycoefh_w, ytke_w,yeps_w,y_treedrg) -! ycoefm_w, ycoefh_w, ytke_w) - IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN -! In this case, coef_diff_turb is called for the Cd only - DO k = 2, klev - DO j = 1, knon - i = ni(j) - ycoefh_w(j,k) = zcoefh(i,k,nsrf) - ycoefm_w(j,k) = zcoefm(i,k,nsrf) - ENDDO - ENDDO - ENDIF - - ENDIF ! iflag_pbl >= 50 - - - IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_w ',ycoefh_w(1:knon,:) - -!!!jyg le 10/04/2013 -!! En attendant de traiter le transport des traceurs dans les poches froides, formule -!! arbitraire pour ycoefh et ycoefm - DO k = 2,klev - DO j = 1,knon - ycoefh(j,k) = ycoefh_x(j,k) + ywake_s(j)*(ycoefh_w(j,k) - ycoefh_x(j,k)) - ycoefm(j,k) = ycoefm_x(j,k) + ywake_s(j)*(ycoefm_w(j,k) - ycoefm_x(j,k)) - ENDDO - ENDDO - - - ENDIF ! (iflag_split .eq.0) - - -!**************************************************************************************** -! -! 8) "La descente" - "The downhill" -! -! climb_hq_down and climb_wind_down calculate the coefficients -! Ccoef_X et Dcoef_X for X=[H, Q, U, V]. -! Only the coefficients at surface for H and Q are returned. -! -!**************************************************************************************** - -! - Calculate the coefficients Ccoef_H, Ccoef_Q, Dcoef_H and Dcoef_Q -!!! jyg le 07/02/2012 - IF (iflag_split .eq.0) THEN -!!! -!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012 - CALL climb_hq_down(knon, ni, ycoefh, ypaprs, ypplay, & - ydelp, yt, yq, dtime, & -!!! jyg le 09/05/2011 - CcoefH, CcoefQ, DcoefH, DcoefQ, & - Kcoef_hq, gama_q, gama_h, & -!!! - AcoefH, AcoefQ, BcoefH, BcoefQ & -#ifdef ISO - & ,yxt, CcoefXT, DcoefXT, gama_xt, AcoefXT, BcoefXT & -#endif - & ) - ELSE !(iflag_split .eq.0) - CALL climb_hq_down(knon, ni, ycoefh_x, ypaprs, ypplay, & - ydelp, yt_x, yq_x, dtime, & -!!! nrlmd le 02/05/2011 - CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x, & - Kcoef_hq_x, gama_q_x, gama_h_x, & -!!! - AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x & -#ifdef ISO - & ,yxt_x, CcoefXT_x, DcoefXT_x, gama_xt_x, AcoefXT_x, BcoefXT_x & -#endif - & ) -!!! - IF (prt_level >=10) THEN - PRINT *,'pbl_surface (climb_hq_down.x->) AcoefH_x ',AcoefH_x - PRINT *,'pbl_surface (climb_hq_down.x->) AcoefQ_x ',AcoefQ_x - PRINT *,'pbl_surface (climb_hq_down.x->) BcoefH_x ',BcoefH_x - PRINT *,'pbl_surface (climb_hq_down.x->) BcoefQ_x ',BcoefQ_x - ENDIF -! - CALL climb_hq_down(knon, ni, ycoefh_w, ypaprs, ypplay, & - ydelp, yt_w, yq_w, dtime, & -!!! nrlmd le 02/05/2011 - CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w, & - Kcoef_hq_w, gama_q_w, gama_h_w, & -!!! - AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w & -#ifdef ISO - & ,yxt_w, CcoefXT_w, DcoefXT_w, gama_xt_w, AcoefXT_w, BcoefXT_w & -#endif - & ) -!!! - IF (prt_level >=10) THEN - PRINT *,'pbl_surface (climb_hq_down.w->) AcoefH_w ',AcoefH_w - PRINT *,'pbl_surface (climb_hq_down.w->) AcoefQ_w ',AcoefQ_w - PRINT *,'pbl_surface (climb_hq_down.w->) BcoefH_w ',BcoefH_w - PRINT *,'pbl_surface (climb_hq_down.w->) BcoefQ_w ',BcoefQ_w - ENDIF -!!! - ENDIF ! (iflag_split .eq.0) -!!! - -! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V -!!! jyg le 07/02/2012 - IF (iflag_split .eq.0) THEN -!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012 - CALL climb_wind_down(knon, ni, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv, & -!!! jyg le 09/05/2011 - CcoefU, CcoefV, DcoefU, DcoefV, & - Kcoef_m, alf_1, alf_2, & -!!! - AcoefU, AcoefV, BcoefU, BcoefV) - ELSE ! (iflag_split .eq.0) - CALL climb_wind_down(knon, ni, dtime, ycoefm_x, ypplay, ypaprs, yt_x, ydelp, yu_x, yv_x, & -!!! nrlmd le 02/05/2011 - CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x, & - Kcoef_m_x, alf_1_x, alf_2_x, & -!!! - AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x) -! - CALL climb_wind_down(knon, ni, dtime, ycoefm_w, ypplay, ypaprs, yt_w, ydelp, yu_w, yv_w, & -!!! nrlmd le 02/05/2011 - CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w, & - Kcoef_m_w, alf_1_w, alf_2_w, & -!!! - AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w) -!!! - ENDIF ! (iflag_split .eq.0) -!!! - -! For blowing snow: - IF (ok_bs) THEN - ! following Bintanja et al 2000, part II and Vionnet V PhD thesis - ! we assume that the eddy diffsivity coefficient for - ! suspended particles is a fraction of Kh - do k=1,klev - do j=1,knon - ycoefqbs(j,k)=ycoefh(j,k)*zeta_bs - enddo - enddo - CALL climb_qbs_down(knon, ni, ycoefqbs, ypaprs, ypplay, & - ydelp, yt, yqbs, dtime, & - CcoefQBS, DcoefQBS, & - Kcoef_qbs, gama_qbs, & - AcoefQBS, BcoefQBS) - ENDIF - -!**************************************************************************************** -! 9) Small calculations -! -!**************************************************************************************** - -! - Reference pressure is given the values at surface level - ypsref(:) = ypaprs(:,1) - -! - CO2 field on 2D grid to be sent to ORCHIDEE -! Transform to compressed field - IF (carbon_cycle_cpl) THEN - DO i=1,knon - r_co2_ppm(i) = co2_send(ni(i)) - ENDDO - ELSE - r_co2_ppm(:) = co2_ppm ! Constant field - ENDIF - -!!! nrlmd le 02/05/2011 -----------------------On raccorde les 2 colonnes dans la couche 1 -!---------------------------------------------------------------------------------------- -!!! jyg le 07/02/2012 -!!! jyg le 01/02/2017 - IF (iflag_split .eq. 0) THEN - yt1(:) = yt(:,1) - yq1(:) = yq(:,1) -#ifdef ISO - yxt1(:,:) = yxt(:,:,1) -#endif - - ELSE IF (iflag_split .ge. 1) THEN -#ifdef ISO - call abort_physic('pbl_surface_mod 2149','isos pas encore dans iflag_split=1',1) -#endif - -! -! Cdragq computation -! ------------------ - !****************************************************************************** - ! Cdragq computed from cdrag - ! The difference comes only from a factor (f_z0qh_oce) on z0, so that - ! it can be computed inside wx_pbl0_merge - ! More complicated appraches may require the propagation through - ! pbl_surface of an independant cdragq variable. - !****************************************************************************** -! - IF ( f_z0qh_oce .ne. 1. .and. nsrf .eq.is_oce) THEN - ! Si on suit les formulations par exemple de Tessel, on - ! a z0h=0.4*nu/u*, z0q=0.62*nu/u*, d'ou f_z0qh_oce=0.62/0.4=1.55 -!! ycdragq_x(1:knon)=ycdragh_x(1:knon)* & -!! log(z1lay(1:knon)/yz0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*yz0h(1:knon))) -!! ycdragq_w(1:knon)=ycdragh_w(1:knon)* & -!! log(z1lay(1:knon)/yz0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*yz0h(1:knon))) -! - DO j = 1,knon - z1lay = zgeo1(j)/RG - fact_cdrag = log(z1lay/yz0h(j))/log(z1lay/(f_z0qh_oce*yz0h(j))) - ycdragq_x(j)=ycdragh_x(j)*fact_cdrag - ycdragq_w(j)=ycdragh_w(j)*fact_cdrag -!! Print *,'YYYYpbl0: fact_cdrag ', fact_cdrag - ENDDO ! j = 1,knon -! -!! Print *,'YYYYpbl0: z1lay, yz0h, f_z0qh_oce, ycdragh_w, ycdragq_w ', & -!! z1lay, yz0h(1:knon), f_z0qh_oce, ycdragh_w(1:knon), ycdragq_w(1:knon) - ELSE - ycdragq_x(1:knon)=ycdragh_x(1:knon) - ycdragq_w(1:knon)=ycdragh_w(1:knon) - ENDIF ! ( f_z0qh_oce .ne. 1. .and. nsrf .eq.is_oce) -! - CALL wx_pbl_prelim_0(knon, nsrf, dtime, ypplay, ypaprs, ywake_s, & - yts, y_delta_tsurf, ygustiness, & - yt_x, yt_w, yq_x, yq_w, & - yu_x, yu_w, yv_x, yv_w, & - ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, & - ycdragm_x, ycdragm_w, & - AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, & - AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, & - BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, & - BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, & - Kech_h_x, Kech_h_w, Kech_h & - ) - CALL wx_pbl_prelim_beta(knon, dtime, ywake_s, ybeta, & - BcoefQ_x, BcoefQ_w & - ) - CALL wx_pbl0_merge(knon, ypplay, ypaprs, & - ywake_s, ydTs0, ydqs0, & - yt_x, yt_w, yq_x, yq_w, & - yu_x, yu_w, yv_x, yv_w, & - ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, & - ycdragm_x, ycdragm_w, & - AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, & - AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, & - BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, & - BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, & - AcoefH_0, AcoefQ_0, AcoefU, AcoefV, & - BcoefH_0, BcoefQ_0, BcoefU, BcoefV, & - ycdragh, ycdragq, ycdragm, & - yt1, yq1, yu1, yv1 & - ) - IF (iflag_split .eq. 2 .AND. nsrf .ne. is_oce) THEN - CALL wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, & - ywake_s, ybeta, ywake_cstar, ywake_dens, & - AcoefH_x, AcoefH_w, & - BcoefH_x, BcoefH_w, & - AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, & - AcoefH, AcoefQ, BcoefH, BcoefQ, & - HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, & - phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, & - yg_T, yg_Q, & - yGamma_dTs_phiT, yGamma_dQs_phiQ, & - ydTs_ins, ydqs_ins & - ) - ELSE ! - AcoefH(:) = AcoefH_0(:) - AcoefQ(:) = AcoefQ_0(:) - BcoefH(:) = BcoefH_0(:) - BcoefQ(:) = BcoefQ_0(:) - yg_T(:) = 0. - yg_Q(:) = 0. - yGamma_dTs_phiT(:) = 0. - yGamma_dQs_phiQ(:) = 0. - ydTs_ins(:) = 0. - ydqs_ins(:) = 0. - ENDIF ! (iflag_split .eq. 2) - ENDIF ! (iflag_split .eq.0) -!!! - IF (prt_level >=10) THEN - DO i = 1, min(1,knon) - PRINT *,'pbl_surface (merge->): yt(1,:) ',yt(i,:) - PRINT *,'pbl_surface (merge->): yq(1,:) ',yq(i,:) - PRINT *,'pbl_surface (merge->): yu(1,:) ',yu(i,:) - PRINT *,'pbl_surface (merge->): yv(1,:) ',yv(i,:) - PRINT *,'pbl_surface (merge->): AcoefH(1), AcoefQ(1), AcoefU(1), AcoefV(1) ', & - AcoefH(i), AcoefQ(i), AcoefU(i), AcoefV(i) - PRINT *,'pbl_surface (merge->): BcoefH(1), BcoefQ(1), BcoefU(1), BcoefV(1) ', & - BcoefH(i), BcoefQ(i), BcoefU(i), BcoefV(i) - ENDDO - - ENDIF - -! Save initial value of z0h for use in evappot (z0h wiil be computed again in the surface models) - yz0h_old(1:knon) = yz0h(1:knon) -! -!**************************************************************************************** -! -! Calulate t2m and q2m for the case of calculation at land grid points -! t2m and q2m are needed as input to ORCHIDEE -! -!**************************************************************************************** - IF (nsrf == is_ter) THEN - - DO i = 1, knon - zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) & - * (ypaprs(i,1)-ypplay(i,1)) - ENDDO - - ! Calculate the temperature et relative humidity at 2m and the wind at 10m - IF (iflag_new_t2mq2m==1) THEN - CALL stdlevvarn(klon, knon, is_ter, zxli, & - yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, & - yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), & - yt2m, yq2m, yt10m, yq10m, yu10m, yustar, & - yn2mout(:, nsrf, :)) - ELSE - CALL stdlevvar(klon, knon, is_ter, zxli, & - yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, & - yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), & - yt2m, yq2m, yt10m, yq10m, yu10m, yustar, ypblh, rain_f, zxtsol) - ENDIF - - ENDIF - -!**************************************************************************************** -! -! 10) Switch according to current surface -! It is necessary to start with the continental surfaces because the ocean -! needs their run-off. -! -!**************************************************************************************** - SELECT CASE(nsrf) - - CASE(is_ter) -! print*,"DEBUGTS",yts(knon/2),ylwdown(knon/2) - CALL surf_land(itap, dtime, date0, jour, knon, ni,& - rlon, rlat, yrmu0, & - debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, & -!!jyg yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),& - yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt1, yq1,& - AcoefH, AcoefQ, BcoefH, BcoefQ, & - AcoefU, AcoefV, BcoefU, BcoefV, & - ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, & - ylwdown, yq2m, yt2m, & - ysnow, yqsol, yagesno, ytsoil, & - yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,yfluxbs,& - yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, & - y_flux_u1, y_flux_v1, & - yveget,ylai,yheight, tsurf_tersrf, tsoil_tersrf, qsurf_tersrf, tsurf_new_tersrf, & - cdragm_tersrf, cdragh_tersrf, & - swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf & -!GG -! yveget,ylai,yheight,hice,tice,bilg_cumul, & -! fcds, fcdi, dh_basal_growth, dh_basal_melt, dh_top_melt, dh_snow2sic, & -! dtice_melt, dtice_snow2sic) - !GG -#ifdef ISO - & ,yxtrain_f, yxtsnow_f,yxt1, & - & yxtsnow,yxtsol,yxtevap,h1, & - & yrunoff_diag,yxtrunoff_diag,yRland_ice & -#endif - & ) - - tsurf_tersrf(:,:) = tsurf_new_tersrf(:,:) ! for next time step - -!FC quid qd yveget ylai yheight ne sont pas definit -!FC yveget,ylai,yheight, & - IF (ifl_pbltree .ge. 1) THEN - CALL freinage(knon, yu, yv, yt, & -! yveget,ylai, yheight,ypaprs,ypplay,y_d_u_frein,y_d_v_frein) - yveget,ylai, yheight,ypaprs,ypplay,y_treedrg, y_d_u_frein,y_d_v_frein) - ENDIF - - -! Special DICE MPL 05082013 puis BOMEX - IF (ok_prescr_ust) THEN - DO j=1,knon -! ysnow(:)=0. -! yqsol(:)=0. -! yagesno(:)=50. -! ytsoil(:,:)=300. -! yz0_new(:)=0.001 -! yevap(:)=flat/RLVTT -! yfluxlat(:)=-flat -! yfluxsens(:)=-fsens -! yqsurf(:)=0. -! ytsurf_new(:)=tg -! y_dflux_t(:)=0. -! y_dflux_q(:)=0. - 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) - 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) - ENDDO - ENDIF - -#ifdef ISOVERIF - DO j=1,knon - DO ixt=1,ntraciso - CALL iso_verif_noNaN(yxtevap(ixt,j), & - & 'pbl_surface 1056a: apres surf_land') - ENDDO - DO ixt=1,niso - CALL iso_verif_noNaN(yxtsol(ixt,j), & - & 'pbl_surface 1056b: apres surf_land') - ENDDO - ENDDO -#endif -#ifdef ISOVERIF -! write(*,*) 'pbl_surface_mod 1038: sortie surf_land' - DO j=1,knon - IF (iso_eau >= 0) THEN - CALL iso_verif_egalite(yxtsnow(iso_eau,j), & - & ysnow(j),'pbl_surf_mod 1043') - ENDIF !if (iso_eau.gt.0) then - ENDDO !DO i=1,klon -#endif - - CASE(is_lic) - ! Martin - - IF (landice_opt .LT. 2) THEN - ! Land ice is treated by LMDZ and not by ORCHIDEE - CALL surf_landice(itap, dtime, knon, ni, & - rlon, rlat, debut, lafin, & - yrmu0, ylwdown, yalb, zgeo1, & - ysolsw, ysollw, yts, ypplay(:,1), & - ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt1, yq1,& - AcoefH, AcoefQ, BcoefH, BcoefQ, & - AcoefU, AcoefV, BcoefU, BcoefV, & - AcoefQBS, BcoefQBS, & - ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, & - ysnow, yqsurf, yqsol,yqbs1, yagesno, & - ytsoil, yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yicesub_lic, yfluxsens,yfluxlat, & - yfluxbs, ytsurf_new, y_dflux_t, y_dflux_q, & - yzmea, yzsig, ycldt, & - ysnowhgt, yqsnow, ytoice, ysissnow, & - yalb3_new, yrunoff, & - y_flux_u1, y_flux_v1 & -#ifdef ISO - & ,yxtrain_f, yxtsnow_f,yxt1,yRland_ice & - & ,yxtsnow,yxtsol,yxtevap & -#endif - & ) - - !jyg< - !! alb3_lic(:)=0. - !>jyg - DO j = 1, knon - i = ni(j) - alb3_lic(i) = yalb3_new(j) - snowhgt(i) = ysnowhgt(j) - qsnow(i) = yqsnow(j) - to_ice(i) = ytoice(j) - sissnow(i) = ysissnow(j) - runoff(i) = yrunoff(j) - icesub_lic(i) = yicesub_lic(j)*ypct(j) - ENDDO - ! Martin - ! Special DICE MPL 05082013 puis BOMEX MPL 20150410 - IF (ok_prescr_ust) THEN - DO j=1,knon - 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) - 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) - ENDDO - ENDIF - -#ifdef ISOVERIF - DO j=1,knon - DO ixt=1,ntraciso - CALL iso_verif_noNaN(yxtevap(ixt,j), & - & 'pbl_surface 1095a: apres surf_landice') - ENDDO - do ixt=1,niso - call iso_verif_noNaN(yxtsol(ixt,j), & - & 'pbl_surface 1095b: apres surf_landice') - enddo - enddo -#endif -#ifdef ISOVERIF - !write(*,*) 'pbl_surface_mod 1060: sortie surf_landice' - do j=1,knon - IF (iso_eau >= 0) THEN - CALL iso_verif_egalite(yxtsnow(iso_eau,j), & - & ysnow(j),'pbl_surf_mod 1064') - ENDIF !if (iso_eau >= 0) THEN - ENDDO !DO i=1,klon -#endif - - END IF - - CASE(is_oce) - -!GG -! calculate length scale PBL - - if (iflag_leads == 1) then - ydthetadz = 999999. - ypphii = 999999. - ytheta = 999999. - - DO k = 1, klev - DO j = 1, knon - ytheta(j,k) = yt(j,k)*(ypplay(j,k)/1.e5)**(RD/RCPD) - ENDDO - ENDDO - - DO k = 2, klev - DO j = 1, knon - ydthetadz(j,k) = RG*( ytheta(j,k) - ytheta(j,k-1) ) / ( ypphi(j,k) - ypphi(j,k-1) ) - ypphii(j,k) = (ypphi(j,k)+ypphi(j,k-1))/(RG*2.) - ENDDO - ENDDO - - DO j = 1, knon - ! print *, "ypphii(j,:)=", ypphii(j,:) - ! print *, "ypplay(j,:)=", ypplay(j,:) - ! print *, "ytheta(j,:)=", ytheta(j,:) - ! print *, "minloc(abs(ypphii(j,:)-300))=", - ! minloc(abs(ypphii(j,:)-300),1) - k= minloc(abs(ypphii(j,:)-300),1) - ydthetadz300(j)=ydthetadz(j,k) - ENDDO - end if -!GG - CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb_vis, & - ywindsp, yrmu0, yfder, yts, & - itap, dtime, jour, knon, ni, & -!!jyg ypplay(:,1), zgeo1/RG, ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),& - ypplay(:,1), zgeo1(1:knon)/RG, ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt(:,1), yq(:,1),& ! ym missing init - AcoefH, AcoefQ, BcoefH, BcoefQ, & - AcoefU, AcoefV, BcoefU, BcoefV, & - ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, & - ysnow, yqsurf, yagesno, & - yz0m, yz0h, SFRWL,yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,& - ytsurf_new, y_dflux_t, y_dflux_q, slab_wfbils, & - y_flux_u1, y_flux_v1, ydelta_sst(:knon), ydelta_sal(:knon), & - yds_ns(:knon), ydt_ns(:knon), ydter(:knon), ydser(:knon), & - !GG ydt_ds(:knon), ytkt(:knon), ytks(:knon), ytaur(:knon), ysss) - ydt_ds(:knon), ytkt(:knon), ytks(:knon), ytaur(:knon), ysss, & - ydthetadz300,Ampl & - !GG -#ifdef ISO - & ,yxtrain_f, yxtsnow_f,yxt1,Roce, & - & yxtsnow,yxtevap,h1 & -#endif - & ) - IF (prt_level >=10) THEN - print *,'arg de surf_ocean: ycdragh ',ycdragh(1:knon) - print *,'arg de surf_ocean: ycdragm ',ycdragm(1:knon) - print *,'arg de surf_ocean: yt ', yt(1:knon,:) - print *,'arg de surf_ocean: yq ', yq(1:knon,:) - print *,'arg de surf_ocean: yts ', yts(1:knon) - print *,'arg de surf_ocean: AcoefH ',AcoefH(1:knon) - print *,'arg de surf_ocean: AcoefQ ',AcoefQ(1:knon) - print *,'arg de surf_ocean: BcoefH ',BcoefH(1:knon) - print *,'arg de surf_ocean: BcoefQ ',BcoefQ(1:knon) - print *,'arg de surf_ocean: yevap ',yevap(1:knon) - print *,'arg de surf_ocean: yfluxsens ',yfluxsens(1:knon) - print *,'arg de surf_ocean: yfluxlat ',yfluxlat(1:knon) - print *,'arg de surf_ocean: ytsurf_new ',ytsurf_new(1:knon) - ENDIF -! Special DICE MPL 05082013 puis BOMEX MPL 20150410 - IF (ok_prescr_ust) THEN - DO j=1,knon - 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) - 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) - ENDDO - ENDIF - - CASE(is_sic) - CALL surf_seaice( & -!albedo SB >>> - rlon, rlat, ysolsw, ysollw, yalb_vis, yfder, & -!albedo SB <<< - itap, dtime, jour, knon, ni, & - lafin, & -!!jyg yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),& - yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,& - AcoefH, AcoefQ, BcoefH, BcoefQ, & - AcoefU, AcoefV, BcoefU, BcoefV, & - ypsref, yu1, yv1, ygustiness, pctsrf, & - ysnow, yqsurf, yqsol, yagesno, ytsoil, & -!albedo SB >>> - yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,& -!albedo SB <<< - ytsurf_new, y_dflux_t, y_dflux_q, & -!GG y_flux_u1, y_flux_v1) - y_flux_u1, y_flux_v1, & - hice,tice,bilg_cumul, & - fcds, fcdi, dh_basal_growth, dh_basal_melt, dh_top_melt, dh_snow2sic, & - dtice_melt, dtice_snow2sic & -!GG -#ifdef ISO - & ,yxtrain_f, yxtsnow_f,yxt1,Roce, & - & yxtsnow,yxtsol,yxtevap,Rland_ice & -#endif - & ) - -! Special DICE MPL 05082013 puis BOMEX MPL 20150410 - IF (ok_prescr_ust) THEN - DO j=1,knon - 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) - 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) - ENDDO - ENDIF - -#ifdef ISOVERIF - DO j=1,knon - DO ixt=1,ntraciso - CALL iso_verif_noNaN(yxtevap(ixt,j), & - & 'pbl_surface 1165a: apres surf_seaice') - ENDDO - DO ixt=1,niso - CALL iso_verif_noNaN(yxtsol(ixt,j), & - & 'pbl_surface 1165b: apres surf_seaice') - ENDDO - ENDDO -#endif -#ifdef ISOVERIF - !write(*,*) 'pbl_surface_mod 1077: sortie surf_seaice' - DO j=1,knon - IF (iso_eau >= 0) THEN - CALL iso_verif_egalite(yxtsnow(iso_eau,j), & - & ysnow(j),'pbl_surf_mod 1106') - ENDIF !IF (iso_eau >= 0) THEN - ENDDO !DO i=1,klon -#endif - - CASE DEFAULT - WRITE(lunout,*) 'Surface index = ', nsrf - abort_message = 'Surface index not valid' - CALL abort_physic(modname,abort_message,1) - END SELECT - - -!**************************************************************************************** -! 11) - Calcul the increment of surface temperature -! -!**************************************************************************************** - - IF (evap0>=0.) THEN - yevap(1:knon)=evap0 - yevap(1:knon)=RLVTT*evap0 - ENDIF - - y_d_ts(1:knon) = ytsurf_new(1:knon) - yts(1:knon) - -!**************************************************************************************** -! -! 12) "La remontee" - "The uphill" -! -! The fluxes (y_flux_X) and tendancy (y_d_X) are calculated -! for X=H, Q, U and V, for all vertical levels. -! -!**************************************************************************************** -!! -!!! -!!! jyg le 10/04/2013 et EV 10/2020 - - IF (ok_forc_tsurf) THEN - DO j=1,knon - ytsurf_new(j)=tg - y_d_ts(j) = ytsurf_new(j) - yts(j) - ENDDO - ENDIF ! ok_forc_tsurf - -!!! - IF (ok_flux_surf) THEN - IF (prt_level >=10) THEN - PRINT *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT - ENDIF - y_flux_t1(:) = fsens - y_flux_q1(:) = flat/RLVTT - yfluxlat(:) = flat -! -!! Test sur iflag_split retire le 2/02/2018, sans vraiment comprendre la raison de ce test. (jyg) -!! IF (iflag_split .eq.0) THEN - DO j=1,knon - Kech_h(j) = ycdragh(j) * (1.0+SQRT(yu(j,1)**2+yv(j,1)**2)) * & - ypplay(j,1)/(RD*yt(j,1)) - ENDDO -!! ENDIF ! (iflag_split .eq.0) - - DO j = 1, knon - yt1_new=(1./RCPD)*(AcoefH(j)+BcoefH(j)*y_flux_t1(j)*dtime) - ytsurf_new(j)=yt1_new-y_flux_t1(j)/(Kech_h(j)*RCPD) - ! for cases forced in flux and for which forcing in Ts is needed - ! to prevent the latter to reach unrealistic value (even if not used, - ! Ts is calculated and hgardfou can appear during the calculation - ! of surface saturation humidity for example - if (ok_forc_tsurf) ytsurf_new(j)=tg - ENDDO - - DO j=1,knon - y_d_ts(j) = ytsurf_new(j) - yts(j) - ENDDO - - ELSE ! (ok_flux_surf) - DO j=1,knon - y_flux_t1(j) = yfluxsens(j) - y_flux_q1(j) = -yevap(j) -#ifdef ISO - y_flux_xt1(:,:) = -yxtevap(:,:) -#endif - ENDDO - ENDIF ! (ok_flux_surf) - - ! flux of blowing snow at the first level - IF (ok_bs) THEN - DO j=1,knon - y_flux_bs(j)=yfluxbs(j) - ENDDO - ENDIF -! -! ------------------------------------------------------------------------------ -! 12a) Splitting -! ------------------------------------------------------------------------------ - - IF (iflag_split .GE. 1) THEN -#ifdef ISO - call abort_physic('pbl_surface_mod 2607','isos pas encore dans iflag_split=1',1) -#endif -! -! - IF (nsrf .ne. is_oce) THEN -! -! Compute potential evaporation and aridity factor (jyg, 20200328) - ybeta_prev(:) = ybeta(:) - DO j = 1, knon - yqa(j) = AcoefQ(j) - BcoefQ(j)*yevap(j)*dtime - ENDDO -! - CALL wx_evappot(knon, yqa, yTsurf_new, yevap_pot) -! - ybeta(1:knon) = min(yevap(1:knon)/yevap_pot(1:knon), 1.) - - IF (prt_level >=10) THEN - DO j=1,knon - print*,'y_flux_t1,yfluxlat,wakes' & - & , y_flux_t1(j), yfluxlat(j), ywake_s(j) - print*,'beta_prev, beta, ytsurf_new', ybeta_prev(j), ybeta(j), ytsurf_new(j) - print*,'inertia,facteur,cstar', inertia, facteur,wake_cstar(j) - ENDDO - ENDIF ! (prt_level >=10) -! -! Second call to wx_pbl0_merge and wx_pbl_dts_merge in order to take into account -! the update of the aridity coeficient beta. -! - CALL wx_pbl_prelim_beta(knon, dtime, ywake_s, ybeta, & - BcoefQ_x, BcoefQ_w & - ) - CALL wx_pbl0_merge(knon, ypplay, ypaprs, & - ywake_s, ydTs0, ydqs0, & - yt_x, yt_w, yq_x, yq_w, & - yu_x, yu_w, yv_x, yv_w, & - ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, & - ycdragm_x, ycdragm_w, & - AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, & - AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, & - BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, & - BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, & - AcoefH_0, AcoefQ_0, AcoefU, AcoefV, & - BcoefH_0, BcoefQ_0, BcoefU, BcoefV, & - ycdragh, ycdragq, ycdragm, & - yt1, yq1, yu1, yv1 & - ) - IF (iflag_split .eq. 2) THEN - CALL wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, & - ywake_s, ybeta, ywake_cstar, ywake_dens, & - AcoefH_x, AcoefH_w, & - BcoefH_x, BcoefH_w, & - AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, & - AcoefH, AcoefQ, BcoefH, BcoefQ, & - HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, & - phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, & - yg_T, yg_Q, & - yGamma_dTs_phiT, yGamma_dQs_phiQ, & - ydTs_ins, ydqs_ins & - ) - ELSE ! - AcoefH(:) = AcoefH_0(:) - AcoefQ(:) = AcoefQ_0(:) - BcoefH(:) = BcoefH_0(:) - BcoefQ(:) = BcoefQ_0(:) - yg_T(:) = 0. - yg_Q(:) = 0. - yGamma_dTs_phiT(:) = 0. - yGamma_dQs_phiQ(:) = 0. - ydTs_ins(:) = 0. - ydqs_ins(:) = 0. - ENDIF ! (iflag_split .eq. 2) -! - ELSE ! (nsrf .ne. is_oce) - ybeta(1:knon) = 1. - yevap_pot(1:knon) = yevap(1:knon) - AcoefH(:) = AcoefH_0(:) - AcoefQ(:) = AcoefQ_0(:) - BcoefH(:) = BcoefH_0(:) - BcoefQ(:) = BcoefQ_0(:) - yg_T(:) = 0. - yg_Q(:) = 0. - yGamma_dTs_phiT(:) = 0. - yGamma_dQs_phiQ(:) = 0. - ydTs_ins(:) = 0. - ydqs_ins(:) = 0. - ENDIF ! (nsrf .ne. is_oce) -! - CALL wx_pbl_split(knon, nsrf, dtime, ywake_s, ybeta, iflag_split, & - yg_T, yg_Q, & - yGamma_dTs_phiT, yGamma_dQs_phiQ, & - ydTs_ins, ydqs_ins, & - y_flux_t1, y_flux_q1, y_flux_u1, y_flux_v1, & -!!!! HTRn_b, dd_HTRn, HTphiT_b, dd_HTphiT, & - phiQ0_b, phiT0_b, & - y_flux_t1_x, y_flux_t1_w, & - y_flux_q1_x, y_flux_q1_w, & - y_flux_u1_x, y_flux_u1_w, & - y_flux_v1_x, y_flux_v1_w, & - yfluxlat_x, yfluxlat_w, & - y_delta_qsats, & - y_delta_tsurf_new, y_delta_qsurf & - ) -! - CALL wx_pbl_check(knon, dtime, ypplay, ypaprs, ywake_s, ybeta, iflag_split, & - yTs, y_delta_tsurf, & - yqsurf, yTsurf_new, & - y_delta_tsurf_new, y_delta_qsats, & - AcoefH_x, AcoefH_w, & - BcoefH_x, BcoefH_w, & - AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, & - AcoefH, AcoefQ, BcoefH, BcoefQ, & - y_flux_t1, y_flux_q1, & - y_flux_t1_x, y_flux_t1_w, & - y_flux_q1_x, y_flux_q1_w) -! - IF (nsrf .ne. is_oce) THEN - CALL wx_pbl_dts_check(knon, dtime, ypplay, ypaprs, ywake_s, ybeta, iflag_split, & - yTs, y_delta_tsurf, & - yqsurf, yTsurf_new, & - y_delta_qsats, y_delta_tsurf_new, y_delta_qsurf, & - AcoefH_x, AcoefH_w, & - BcoefH_x, BcoefH_w, & - AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, & - AcoefH, AcoefQ, BcoefH, BcoefQ, & - HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, & - phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, & - yg_T, yg_Q, & - yGamma_dTs_phiT, yGamma_dQs_phiQ, & - ydTs_ins, ydqs_ins, & - y_flux_t1, y_flux_q1, & - y_flux_t1_x, y_flux_t1_w, & - y_flux_q1_x, y_flux_q1_w ) - ENDIF ! (nsrf .ne. is_oce) -! - ELSE ! (iflag_split .ge. 1) - ybeta(1:knon) = 1. - yevap_pot(1:knon) = yevap(1:knon) - ENDIF ! (iflag_split .ge. 1) -! - IF (prt_level >= 10) THEN - print *,'pbl_surface, ybeta , yevap, yevap_pot ', & - ybeta(1:knon) , yevap(1:knon), yevap_pot(1:knon) - ENDIF ! (prt_level >= 10) -! -!>jyg -! - -!!jyg!! A reprendre apres reflexion =============================================== -!!jyg!! -!!jyg!! DO j=1,knon -!!jyg!!!!! nrlmd le 13/06/2011 -!!jyg!! -!!jyg!!!----Diffusion dans le sol dans le cas continental seulement -!!jyg!! IF (nsrf.eq.is_ter) THEN -!!jyg!!!----Calcul du coefficient delta_coeff -!!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))) -!!jyg!! -!!jyg!!! delta_coef(j)=dtime/(inertia*sqrt(tau_eq(j))) -!!jyg!! delta_coef(j)=facteur*sqrt(tau_eq(j))/inertia -!!jyg!!! delta_coef(j)=0. -!!jyg!! ELSE -!!jyg!! delta_coef(j)=0. -!!jyg!! ENDIF -!!jyg!! -!!jyg!!!----Calcul de delta_tsurf -!!jyg!! y_delta_tsurf(j)=delta_coef(j)*y_delta_flux_t1(j) -!!jyg!! -!!jyg!!!----Si il n'y a pas des poches... -!!jyg!! IF (wake_cstar(j).le.0.01) THEN -!!jyg!! y_delta_tsurf(j)=0. -!!jyg!! y_delta_flux_t1(j)=0. -!!jyg!! ENDIF -!!jyg!! -!!jyg!!!-----Calcul de ybeta (evap_r\'eelle/evap_potentielle) -!!jyg!!!!!!! jyg le 23/02/2012 -!!jyg!!!!!!! -!!jyg!!!! ybeta(j)=y_flux_q1(j) / & -!!jyg!!!! & (Kech_h(j)*(yq(j,1)-yqsatsurf(j))) -!!jyg!!!!!! ybeta(j)=-1.*yevap(j) / & -!!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))) -!!jyg!!!!!!! fin jyg -!!jyg!!!!! -!!jyg!! -!!jyg!! ENDDO -!!jyg!! -!!jyg!!!!! fin nrlmd le 13/06/2011 -!!jyg!! - IF (iflag_split .ge. 1) THEN - IF (prt_level >=10) THEN - DO j = 1, knon - print*,'Chx,Chw,Ch', ycdragh_x(j), ycdragh_w(j), ycdragh(j) - print*,'Khx,Khw,Kh', Kech_h_x(j), Kech_h_w(j), Kech_h(j) - print*,'t1x, t1w, t1, t1_ancien', & - & yt_x(j,1), yt_w(j,1), yt(j,1), t(j,1) - print*,'delta_coef,delta_flux,delta_tsurf,tau', delta_coef(j), y_delta_flux_t1(j), y_delta_tsurf(j), tau_eq(j) - ENDDO - - DO j=1,knon - print*,'fluxT_x, fluxT_w, y_flux_t1, fluxQ_x, fluxQ_w, yfluxlat, wakes' & - & , 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) - print*,'beta, ytsurf_new ', ybeta(j), ytsurf_new(j) - print*,'inertia, facteur, cstar', inertia, facteur,wake_cstar(j) - ENDDO - ENDIF ! (prt_level >=10) - -!!! jyg le 07/02/2012 - ENDIF ! (iflag_split .ge.1) -!!! - -!!! jyg le 07/02/2012 - IF (iflag_split .eq.0) THEN -!!! -!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012 - CALL climb_hq_up(knon, ni, dtime, yt, yq, & - y_flux_q1, y_flux_t1, ypaprs, ypplay, & -!!! jyg le 07/02/2012 - AcoefH, AcoefQ, BcoefH, BcoefQ, & - CcoefH, CcoefQ, DcoefH, DcoefQ, & - Kcoef_hq, gama_q, gama_h, & -!!! - y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:) & -#ifdef ISO - & ,yxt,y_flux_xt1 & - & ,AcoefXT,BcoefXT,CcoefXT,DcoefXT,gama_xt & - & ,y_flux_xt(:,:,:),y_d_xt(:,:,:) & -#endif - & ) - ELSE !(iflag_split .eq.0) - CALL climb_hq_up(knon, ni, dtime, yt_x, yq_x, & - y_flux_q1_x, y_flux_t1_x, ypaprs, ypplay, & -!!! nrlmd le 02/05/2011 - AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x, & - CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x, & - Kcoef_hq_x, gama_q_x, gama_h_x, & -!!! - y_flux_q_x(:,:), y_flux_t_x(:,:), y_d_q_x(:,:), y_d_t_x(:,:) & -#ifdef ISO - & ,yxt_x,y_flux_xt1_x & - & ,AcoefXT_x,BcoefXT_x,CcoefXT_x,DcoefXT_x,gama_xt_x & - & ,y_flux_xt_x(:,:,:),y_d_xt_x(:,:,:) & -#endif - & ) -! - CALL climb_hq_up(knon, ni, dtime, yt_w, yq_w, & - y_flux_q1_w, y_flux_t1_w, ypaprs, ypplay, & -!!! nrlmd le 02/05/2011 - AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w, & - CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w, & - Kcoef_hq_w, gama_q_w, gama_h_w, & -!!! - y_flux_q_w(:,:), y_flux_t_w(:,:), y_d_q_w(:,:), y_d_t_w(:,:) & -#ifdef ISO - & ,yxt_w,y_flux_xt1_w & - & ,AcoefXT_w,BcoefXT_w,CcoefXT_w,DcoefXT_w,gama_xt_w & - & ,y_flux_xt_w(:,:,:),y_d_xt_w(:,:,:) & -#endif - & ) -!!! - ENDIF ! (iflag_split .eq.0) -!!! - -!!! jyg le 07/02/2012 - IF (iflag_split .eq.0) THEN -!!! -!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012 - CALL climb_wind_up(knon, ni, dtime, yu, yv, y_flux_u1, y_flux_v1, & -!!! jyg le 07/02/2012 - AcoefU, AcoefV, BcoefU, BcoefV, & - CcoefU, CcoefV, DcoefU, DcoefV, & - Kcoef_m, & -!!! - y_flux_u, y_flux_v, y_d_u, y_d_v) - y_d_t_diss(:,:)=0. - IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN - CALL yamada_c(knon, knon,dtime,ypaprs,ypplay & - & ,yu,yv,yt,y_d_u,y_d_v,y_d_t,ycdragm,ytke,ycoefm,ycoefh,ycoefq,y_d_t_diss,yustar & - & ,iflag_pbl) - ENDIF -! print*,'yamada_c OK' - - ELSE !(iflag_split .eq.0) - CALL climb_wind_up(knon, ni, dtime, yu_x, yv_x, y_flux_u1_x, y_flux_v1_x, & -!!! nrlmd le 02/05/2011 - AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x, & - CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x, & - Kcoef_m_x, & -!!! - y_flux_u_x, y_flux_v_x, y_d_u_x, y_d_v_x) -! - y_d_t_diss_x(:,:)=0. - IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN - CALL yamada_c(knon, knon,dtime,ypaprs,ypplay & - & ,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 & - ,ycoefq_x,y_d_t_diss_x,yustar_x & - & ,iflag_pbl) - ENDIF -! print*,'yamada_c OK' - - CALL climb_wind_up(knon, ni, dtime, yu_w, yv_w, y_flux_u1_w, y_flux_v1_w, & -!!! nrlmd le 02/05/2011 - AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w, & - CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w, & - Kcoef_m_w, & -!!! - y_flux_u_w, y_flux_v_w, y_d_u_w, y_d_v_w) -!!! - y_d_t_diss_w(:,:)=0. - IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN - CALL yamada_c(knon, knon,dtime,ypaprs,ypplay & - & ,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 & - ,ycoefq_w,y_d_t_diss_w,yustar_w & - & ,iflag_pbl) - ENDIF -! print*,'yamada_c OK' -! - IF (prt_level >=10) THEN - print *, 'After climbing up, lfuxlat_x, fluxlat_w ', & - yfluxlat_x(1:knon), yfluxlat_w(1:knon) - ENDIF -! - ENDIF ! (iflag_split .eq.0) - - IF (ok_bs) THEN - CALL climb_qbs_up(knon, ni, dtime, yqbs, & - y_flux_bs, ypaprs, ypplay, & - AcoefQBS, BcoefQBS, & - CcoefQBS, DcoefQBS, & - Kcoef_qbs, gama_qbs, & - y_flux_qbs(:,:), y_d_qbs(:,:)) - ENDIF - -!!! -!! -!! DO j = 1, knon -!! y_dflux_t(j) = y_dflux_t(j) * ypct(j) -!! y_dflux_q(j) = y_dflux_q(j) * ypct(j) -!! ENDDO -!! -!**************************************************************************************** -! 13) Transform variables for output format : -! - Decompress -! - Multiply with pourcentage of current surface -! - Cumulate in global variable -! -!**************************************************************************************** - - -!!! jyg le 07/02/2012 - IF (iflag_split.EQ.0) THEN -!!! - DO k = 1, klev - DO j = 1, knon - i = ni(j) - y_d_t_diss(j,k) = y_d_t_diss(j,k) * ypct(j) - y_d_t(j,k) = y_d_t(j,k) * ypct(j) - y_d_q(j,k) = y_d_q(j,k) * ypct(j) - y_d_u(j,k) = y_d_u(j,k) * ypct(j) - y_d_v(j,k) = y_d_v(j,k) * ypct(j) -!FC - IF (nsrf .EQ. is_ter .and. ifl_pbltree .GE. 1) THEN -! if (y_d_u_frein(j,k).ne.0. ) then -! print*, nsrf,'IS_TER ++', y_d_u_frein(j,k)*ypct(j),y_d_u(j,k),j,k -! ENDIF - y_d_u(j,k) =y_d_u(j,k) + y_d_u_frein(j,k)*ypct(j) - y_d_v(j,k) =y_d_v(j,k) + y_d_v_frein(j,k)*ypct(j) - treedrg(i,k,nsrf)=y_treedrg(j,k) - ELSE - treedrg(i,k,nsrf)=0. - ENDIF -!FC - flux_t(i,k,nsrf) = y_flux_t(j,k) - flux_q(i,k,nsrf) = y_flux_q(j,k) - flux_u(i,k,nsrf) = y_flux_u(j,k) - flux_v(i,k,nsrf) = y_flux_v(j,k) - -#ifdef ISO - DO ixt=1,ntraciso - y_d_xt(ixt,j,k) = y_d_xt(ixt,j,k) * ypct(j) - flux_xt(ixt,i,k,nsrf) = y_flux_xt(ixt,j,k) - ENDDO ! DO ixt=1,ntraciso - h1_diag(i)=h1(j) -#endif - - ENDDO - ENDDO - -#ifdef ISO -#ifdef ISOVERIF - if (iso_eau.gt.0) then - call iso_verif_egalite_vect2D( & - y_d_xt,y_d_q, & - 'pbl_surface_mod 2600',ntraciso,klon,klev) - endif -#endif -#endif - - ELSE !(iflag_split .eq.0) - -! Tendances hors poches - DO k = 1, klev - DO j = 1, knon - i = ni(j) - y_d_t_diss_x(j,k) = y_d_t_diss_x(j,k) * ypct(j) - y_d_t_x(j,k) = y_d_t_x(j,k) * ypct(j) - y_d_q_x(j,k) = y_d_q_x(j,k) * ypct(j) - y_d_u_x(j,k) = y_d_u_x(j,k) * ypct(j) - y_d_v_x(j,k) = y_d_v_x(j,k) * ypct(j) - - flux_t_x(i,k,nsrf) = y_flux_t_x(j,k) - flux_q_x(i,k,nsrf) = y_flux_q_x(j,k) - flux_u_x(i,k,nsrf) = y_flux_u_x(j,k) - flux_v_x(i,k,nsrf) = y_flux_v_x(j,k) - -#ifdef ISO - DO ixt=1,ntraciso - y_d_xt_x(ixt,j,k) = y_d_xt_x(ixt,j,k) * ypct(j) - flux_xt_x(ixt,i,k,nsrf) = y_flux_xt_x(ixt,j,k) - ENDDO ! DO ixt=1,ntraciso -#endif - ENDDO - ENDDO - -! Tendances dans les poches - DO k = 1, klev - DO j = 1, knon - i = ni(j) - y_d_t_diss_w(j,k) = y_d_t_diss_w(j,k) * ypct(j) - y_d_t_w(j,k) = y_d_t_w(j,k) * ypct(j) - y_d_q_w(j,k) = y_d_q_w(j,k) * ypct(j) - y_d_u_w(j,k) = y_d_u_w(j,k) * ypct(j) - y_d_v_w(j,k) = y_d_v_w(j,k) * ypct(j) - - flux_t_w(i,k,nsrf) = y_flux_t_w(j,k) - flux_q_w(i,k,nsrf) = y_flux_q_w(j,k) - flux_u_w(i,k,nsrf) = y_flux_u_w(j,k) - flux_v_w(i,k,nsrf) = y_flux_v_w(j,k) - -#ifdef ISO - DO ixt=1,ntraciso - y_d_xt_w(ixt,j,k) = y_d_xt_w(ixt,j,k) * ypct(j) - flux_xt_w(ixt,i,k,nsrf) = y_flux_xt_w(ixt,j,k) - ENDDO ! do ixt=1,ntraciso -#endif - - ENDDO - ENDDO - -! Flux, tendances et Tke moyenne dans la maille - DO k = 1, klev - DO j = 1, knon - i = ni(j) - 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)) - 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)) - 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)) - 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)) -#ifdef ISO - DO ixt=1,ntraciso - 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)) - ENDDO ! do ixt=1,ntraciso -#endif - ENDDO - ENDDO - DO j=1,knon - yfluxlat(j)=yfluxlat_x(j)+ywake_s(j)*(yfluxlat_w(j)-yfluxlat_x(j)) - ENDDO - IF (prt_level >=10) THEN - print *,' nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf) ', & - nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf) - ENDIF - - DO k = 1, klev - DO j = 1, knon - 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)) - 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)) - 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)) - 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)) - 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)) - ENDDO - ENDDO - - ENDIF ! (iflag_split .eq.0) - - - ! tendencies of blowing snow - IF (ok_bs) THEN - DO k = 1, klev - DO j = 1, knon - i = ni(j) - y_d_qbs(j,k)=y_d_qbs(j,k) * ypct(j) - flux_qbs(i,k,nsrf) = y_flux_qbs(j,k) - ENDDO - ENDDO - ENDIF - - - DO j = 1, knon - i = ni(j) - evap(i,nsrf) = - flux_q(i,1,nsrf) !jyg - if (ok_bs) then ; snowerosion(i,nsrf)=flux_qbs(i,1,nsrf); endif - beta(i,nsrf) = ybeta(j) !jyg - d_ts(i,nsrf) = y_d_ts(j) -!albedo SB >>> - DO k=1,nsw - alb_dir(i,k,nsrf) = yalb_dir_new(j,k) - alb_dif(i,k,nsrf) = yalb_dif_new(j,k) - ENDDO -!albedo SB <<< - snow(i,nsrf) = ysnow(j) - qsurf(i,nsrf) = yqsurf(j) - z0m(i,nsrf) = yz0m(j) - z0h(i,nsrf) = yz0h(j) - fluxlat(i,nsrf) = yfluxlat(j) - agesno(i,nsrf) = yagesno(j) - cdragh(i) = cdragh(i) + ycdragh(j)*ypct(j) - cdragm(i) = cdragm(i) + ycdragm(j)*ypct(j) - dflux_t(i) = dflux_t(i) + y_dflux_t(j)*ypct(j) - dflux_q(i) = dflux_q(i) + y_dflux_q(j)*ypct(j) -#ifdef ISO - DO ixt=1,niso - xtsnow(ixt,i,nsrf) = yxtsnow(ixt,j) - ENDDO - DO ixt=1,ntraciso - xtevap(ixt,i,nsrf) = - flux_xt(ixt,i,1,nsrf) - dflux_xt(ixt,i) = dflux_xt(ixt,i) + y_dflux_xt(ixt,j)*ypct(j) - ENDDO - IF (nsrf == is_lic) THEN - DO ixt=1,niso - Rland_ice(ixt,i) = yRland_ice(ixt,j) - ENDDO - ENDIF !IF (nsrf == is_lic) THEN -#ifdef ISOVERIF - IF (iso_eau.gt.0) THEN - call iso_verif_egalite_choix(Rland_ice(iso_eau,i),1.0, & - & 'pbl_surf_mod 1230',errmax,errmaxrel) - ENDIF !if (iso_eau.gt.0) then -#endif -#endif - ENDDO - -! print*,'Dans pbl OK2' - -!!! jyg le 07/02/2012 - IF (iflag_split .ge.1) THEN -!!! -!!! nrlmd le 02/05/2011 - DO j = 1, knon - i = ni(j) - fluxlat_x(i,nsrf) = yfluxlat_x(j) - fluxlat_w(i,nsrf) = yfluxlat_w(j) -!!! -!!! nrlmd le 13/06/2011 -!!jyg20170131 delta_tsurf(i,nsrf)=y_delta_tsurf(j)*ypct(j) -!!jyg20210118 delta_tsurf(i,nsrf)=y_delta_tsurf(j) - delta_tsurf(i,nsrf)=y_delta_tsurf_new(j) -! - delta_qsurf(i) = delta_qsurf(i) + y_delta_qsurf(j)*ypct(j) -! - cdragh_x(i) = cdragh_x(i) + ycdragh_x(j)*ypct(j) - cdragh_w(i) = cdragh_w(i) + ycdragh_w(j)*ypct(j) - cdragm_x(i) = cdragm_x(i) + ycdragm_x(j)*ypct(j) - cdragm_w(i) = cdragm_w(i) + ycdragm_w(j)*ypct(j) - kh(i) = kh(i) + Kech_h(j)*ypct(j) - kh_x(i) = kh_x(i) + Kech_h_x(j)*ypct(j) - kh_w(i) = kh_w(i) + Kech_h_w(j)*ypct(j) -!!! - ENDDO -!!! - ENDIF ! (iflag_split .ge.1) -!!! -!!! nrlmd le 02/05/2011 -!!jyg le 20/02/2011 -!! tke_x(:,:,nsrf)=0. -!! tke_w(:,:,nsrf)=0. -!!jyg le 20/02/2011 -!! DO k = 1, klev+1 -!! DO j = 1, knon -!! i = ni(j) -!! wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k) -!! tke(i,k,nsrf) = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf) -!! ENDDO -!! ENDDO -!!jyg le 20/02/2011 -!! DO k = 1, klev+1 -!! DO j = 1, knon -!! i = ni(j) -!! tke(i,k,nsrf)=(1.-ywake_s(j))*tke_x(i,k,nsrf)+ywake_s(j)*tke_w(i,k,nsrf) -!! ENDDO -!! ENDDO -!!! - IF (iflag_split .eq.0) THEN - wake_dltke(:,:,nsrf) = 0. - DO k = 1, klev+1 - DO j = 1, knon - i = ni(j) -!jyg< -!! tke(i,k,nsrf) = ytke(j,k) -!! tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)*ypct(j) - tke_x(i,k,nsrf) = ytke(j,k) - tke_x(i,k,is_ave) = tke_x(i,k,is_ave) + ytke(j,k)*ypct(j) - eps_x(i,k,nsrf) = yeps(j,k) - eps_x(i,k,is_ave) = eps_x(i,k,is_ave) + yeps(j,k)*ypct(j) -!>jyg - ENDDO - ENDDO - - ELSE ! (iflag_split .eq.0) - DO k = 1, klev+1 - DO j = 1, knon - i = ni(j) - wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k) -!jyg< -!! tke(i,k,nsrf) = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf) -!! tke(i,k,is_ave) = tke(i,k,is_ave) + tke(i,k,nsrf)*ypct(j) - tke_x(i,k,nsrf) = ytke_x(j,k) - tke_x(i,k,is_ave) = tke_x(i,k,is_ave) + tke_x(i,k,nsrf)*ypct(j) - eps_x(i,k,nsrf) = yeps_x(j,k) - eps_x(i,k,is_ave) = eps_x(i,k,is_ave) + eps_x(i,k,nsrf)*ypct(j) - wake_dltke(i,k,is_ave) = wake_dltke(i,k,is_ave) + wake_dltke(i,k,nsrf)*ypct(j) - - -!>jyg - ENDDO - ENDDO - ENDIF ! (iflag_split .eq.0) -!!! - DO k = 2, klev - DO j = 1, knon - i = ni(j) - zcoefh(i,k,nsrf) = ycoefh(j,k) - zcoefm(i,k,nsrf) = ycoefm(j,k) - zcoefh(i,k,is_ave) = zcoefh(i,k,is_ave) + ycoefh(j,k)*ypct(j) - zcoefm(i,k,is_ave) = zcoefm(i,k,is_ave) + ycoefm(j,k)*ypct(j) - ENDDO - ENDDO - -! print*,'Dans pbl OK3' - - IF ( nsrf .EQ. is_ter ) THEN - DO j = 1, knon - i = ni(j) - qsol(i) = yqsol(j) -#ifdef ISO - runoff_diag(i)=yrunoff_diag(j) - DO ixt=1,niso - xtsol(ixt,i) = yxtsol(ixt,j) - xtrunoff_diag(ixt,i)=yxtrunoff_diag(ixt,j) - ENDDO -#endif - ENDDO - ENDIF - -!jyg< -!! ftsoil(:,:,nsrf) = 0. -!>jyg - DO k = 1, nsoilmx - DO j = 1, knon - i = ni(j) - ftsoil(i, k, nsrf) = ytsoil(j,k) - ENDDO - ENDDO - -#ifdef ISO -#ifdef ISOVERIF - !write(*,*) 'pbl_surface 2858' - DO i = 1, klon - DO ixt=1,niso - call iso_verif_noNaN(xtsol(ixt,i),'pbl_surface 1405') - ENDDO - ENDDO -#endif -#ifdef ISOVERIF - IF (iso_eau.gt.0) THEN - call iso_verif_egalite_vect2D( & - y_d_xt,y_d_q, & - 'pbl_surface_mod 1261',ntraciso,klon,klev) - ENDIF !if (iso_eau.gt.0) then -#endif -#endif -!!! jyg le 07/02/2012 - IF (iflag_split .ge.1) THEN -!!! -!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012 - DO k = 1, klev - DO j = 1, knon - i = ni(j) - d_t_diss_x(i,k) = d_t_diss_x(i,k) + y_d_t_diss_x(j,k) - d_t_x(i,k) = d_t_x(i,k) + y_d_t_x(j,k) - d_q_x(i,k) = d_q_x(i,k) + y_d_q_x(j,k) - d_u_x(i,k) = d_u_x(i,k) + y_d_u_x(j,k) - d_v_x(i,k) = d_v_x(i,k) + y_d_v_x(j,k) -! - d_t_diss_w(i,k) = d_t_diss_w(i,k) + y_d_t_diss_w(j,k) - d_t_w(i,k) = d_t_w(i,k) + y_d_t_w(j,k) - d_q_w(i,k) = d_q_w(i,k) + y_d_q_w(j,k) - d_u_w(i,k) = d_u_w(i,k) + y_d_u_w(j,k) - d_v_w(i,k) = d_v_w(i,k) + y_d_v_w(j,k) -#ifdef ISO - DO ixt=1,ntraciso - d_xt_x(ixt,i,k) = d_xt_x(ixt,i,k) + y_d_xt_x(ixt,j,k) - d_xt_w(ixt,i,k) = d_xt_w(ixt,i,k) + y_d_xt_w(ixt,j,k) - ENDDO ! DO ixt=1,ntraciso -#endif - -! -!! d_wake_dlt(i,k) = d_wake_dlt(i,k) + y_d_t_w(i,k)-y_d_t_x(i,k) -!! d_wake_dlq(i,k) = d_wake_dlq(i,k) + y_d_q_w(i,k)-y_d_q_x(i,k) - ENDDO - ENDDO -!!! - ENDIF ! (iflag_split .ge.1) -!!! - - DO k = 1, klev - DO j = 1, knon - i = ni(j) - d_t_diss(i,k) = d_t_diss(i,k) + y_d_t_diss(j,k) - d_t(i,k) = d_t(i,k) + y_d_t(j,k) - d_q(i,k) = d_q(i,k) + y_d_q(j,k) -#ifdef ISO - DO ixt=1,ntraciso - d_xt(ixt,i,k) = d_xt(ixt,i,k) + y_d_xt(ixt,j,k) - ENDDO !DO ixt=1,ntraciso -#endif - d_u(i,k) = d_u(i,k) + y_d_u(j,k) - d_v(i,k) = d_v(i,k) + y_d_v(j,k) - ENDDO - ENDDO - - - IF (ok_bs) THEN - DO k = 1, klev - DO j = 1, knon - i = ni(j) - d_qbs(i,k) = d_qbs(i,k) + y_d_qbs(j,k) - ENDDO - ENDDO - ENDIF - -#ifdef ISO -#ifdef ISOVERIF -! write(*,*) 'd_q,d_xt(iso_eau,554,19)=',d_q(554,19),d_xt(iso_eau,554,19) -! write(*,*) 'pbl_surface 2929: d_q,d_xt(iso_eau,2,1)=',d_q(2,1),d_xt(iso_eau,2,1) -! write(*,*) 'y_d_q,y_d_xt(iso_eau,2,1)=',y_d_q(2,1),y_d_xt(iso_eau,2,1) -! write(*,*) 'iso_eau.gt.0=',iso_eau.gt.0 - call iso_verif_noNaN_vect2D( & - & d_xt, & - & 'pbl_surface 1385',ntraciso,klon,klev) - IF (iso_eau >= 0) THEN - call iso_verif_egalite_vect2D( & - y_d_xt,y_d_q, & - 'pbl_surface_mod 2945',ntraciso,klon,klev) - call iso_verif_egalite_vect2D( & - d_xt,d_q, & - 'pbl_surface_mod 1276',ntraciso,klon,klev) - ENDIF !IF (iso_eau >= 0) THEN -#endif -#endif - -! print*,'Dans pbl OK4' - - IF (prt_level >=10) THEN - print *, 'pbl_surface tendencies for w: d_t_w, d_t_x, d_t ', & - d_t_w(1:knon,1), d_t_x(1:knon,1), d_t(1:knon,1) - ENDIF - - if (nsrf == is_oce .and. activate_ocean_skin >= 1) then - delta_sal = missing_val - ds_ns = missing_val - dt_ns = missing_val - delta_sst = missing_val - dter = missing_val - dser = missing_val - tkt = missing_val - tks = missing_val - taur = missing_val - sss = missing_val - - delta_sal(ni(:knon)) = ydelta_sal(:knon) - ds_ns(ni(:knon)) = yds_ns(:knon) - dt_ns(ni(:knon)) = ydt_ns(:knon) - delta_sst(ni(:knon)) = ydelta_sst(:knon) - dter(ni(:knon)) = ydter(:knon) - dser(ni(:knon)) = ydser(:knon) - tkt(ni(:knon)) = ytkt(:knon) - tks(ni(:knon)) = ytks(:knon) - taur(ni(:knon)) = ytaur(:knon) - sss(ni(:knon)) = ysss(:knon) - - if (activate_ocean_skin == 2 .and. type_ocean == "couple") then - dt_ds = missing_val - dt_ds(ni(:knon)) = ydt_ds(:knon) - end if - end if - - -!**************************************************************************************** -! 14) Calculate the temperature and relative humidity at 2m and the wind at 10m -! Call HBTM -! -!**************************************************************************************** -!!! -! -#undef T2m -#define T2m -#ifdef T2m -! Calculations of diagnostic t,q at 2m and u, v at 10m - -! print*,'Dans pbl OK41' -! print*,'tair1,yt(:,1),y_d_t(:,1)' -! print*, tair1,yt(:,1),y_d_t(:,1) -!!! jyg le 07/02/2012 - IF (iflag_split .eq.0) THEN - DO j=1, knon - uzon(j) = yu(j,1) + y_d_u(j,1) - vmer(j) = yv(j,1) + y_d_v(j,1) - tair1(j) = yt(j,1) + y_d_t(j,1) + y_d_t_diss(j,1) - qair1(j) = yq(j,1) + y_d_q(j,1) - zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) & - * (ypaprs(j,1)-ypplay(j,1)) - tairsol(j) = yts(j) + y_d_ts(j) - qairsol(j) = yqsurf(j) - ENDDO - ELSE ! (iflag_split .eq.0) - DO j=1, knon - uzon_x(j) = yu_x(j,1) + y_d_u_x(j,1) - vmer_x(j) = yv_x(j,1) + y_d_v_x(j,1) - tair1_x(j) = yt_x(j,1) + y_d_t_x(j,1) + y_d_t_diss_x(j,1) - qair1_x(j) = yq_x(j,1) + y_d_q_x(j,1) - zgeo1_x(j) = RD * tair1_x(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) & - * (ypaprs(j,1)-ypplay(j,1)) - tairsol(j) = yts(j) + y_d_ts(j) -!! tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf(j) - tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf_new(j) - qairsol(j) = yqsurf(j) - ENDDO - DO j=1, knon - uzon_w(j) = yu_w(j,1) + y_d_u_w(j,1) - vmer_w(j) = yv_w(j,1) + y_d_v_w(j,1) - tair1_w(j) = yt_w(j,1) + y_d_t_w(j,1) + y_d_t_diss_w(j,1) - qair1_w(j) = yq_w(j,1) + y_d_q_w(j,1) - zgeo1_w(j) = RD * tair1_w(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) & - * (ypaprs(j,1)-ypplay(j,1)) - tairsol_w(j) = tairsol(j) + (1.- ywake_s(j))*y_delta_tsurf(j) - qairsol(j) = yqsurf(j) - ENDDO -!!! - ENDIF ! (iflag_split .eq.0) -!!! - DO j=1, knon -! i = ni(j) -! yz0h_oupas(j) = yz0m(j) -! IF(nsrf.EQ.is_oce) THEN -! yz0h_oupas(j) = z0m(i,nsrf) -! ENDIF - psfce(j)=ypaprs(j,1) - patm(j)=ypplay(j,1) - ENDDO - - IF (iflag_pbl_surface_t2m_bug==1) THEN - yz0h_oupas(1:knon)=yz0m(1:knon) - ELSE - yz0h_oupas(1:knon)=yz0h(1:knon) - ENDIF - -! print*,'Dans pbl OK42A' -! print*,'tair1,yt(:,1),y_d_t(:,1)' -! print*, tair1,yt(:,1),y_d_t(:,1) - -! Calculate the temperature and relative humidity at 2m and the wind at 10m -!!! jyg le 07/02/2012 - IF (iflag_split .eq.0) THEN - IF (iflag_new_t2mq2m==1) THEN - CALL stdlevvarn(klon, knon, nsrf, zxli, & - uzon, vmer, tair1, qair1, zgeo1, & - tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, & - yt2m, yq2m, yt10m, yq10m, yu10m, yustar, & - yn2mout(:, nsrf, :)) - ELSE - CALL stdlevvar(klon, knon, nsrf, zxli, & - uzon, vmer, tair1, qair1, zgeo1, & - tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, & - yt2m, yq2m, yt10m, yq10m, yu10m, yustar, ypblh, rain_f, zxtsol) - ENDIF - ELSE !(iflag_split .eq.0) - IF (iflag_new_t2mq2m==1) THEN - CALL stdlevvarn(klon, knon, nsrf, zxli, & - uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, & - tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, & - yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x, & - yn2mout_x(:, nsrf, :)) - CALL stdlevvarn(klon, knon, nsrf, zxli, & - uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, & - tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, & - yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w, & - yn2mout_w(:, nsrf, :)) - ELSE - CALL stdlevvar(klon, knon, nsrf, zxli, & - uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, & - tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, & - yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x, ypblh_x, rain_f, zxtsol) - CALL stdlevvar(klon, knon, nsrf, zxli, & - uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, & - tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, & - yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w, ypblh_w, rain_f, zxtsol) - ENDIF -!!! - ENDIF ! (iflag_split .eq.0) -!!! -!!! jyg le 07/02/2012 - IF (iflag_split .eq.0) THEN - DO j=1, knon - i = ni(j) - t2m(i,nsrf)=yt2m(j) - q2m(i,nsrf)=yq2m(j) - ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman - ustar(i,nsrf)=yustar(j) - u10m(i,nsrf)=(yu10m(j) * uzon(j))/max(SQRT(uzon(j)**2+vmer(j)**2), smallestreal) - v10m(i,nsrf)=(yu10m(j) * vmer(j))/max(SQRT(uzon(j)**2+vmer(j)**2), smallestreal) -! - DO k = 1, 6 - n2mout(i,nsrf,k) = yn2mout(j,nsrf,k) - END DO -! - ENDDO - ELSE !(iflag_split .eq.0) - DO j=1, knon - i = ni(j) - t2m_x(i,nsrf)=yt2m_x(j) - q2m_x(i,nsrf)=yq2m_x(j) - ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman - ustar_x(i,nsrf)=yustar_x(j) - u10m_x(i,nsrf)=(yu10m_x(j) * uzon_x(j))/max(SQRT(uzon_x(j)**2+vmer_x(j)**2), smallestreal) - v10m_x(i,nsrf)=(yu10m_x(j) * vmer_x(j))/max(SQRT(uzon_x(j)**2+vmer_x(j)**2), smallestreal) -! - DO k = 1, 6 - n2mout_x(i,nsrf,k) = yn2mout_x(j,nsrf,k) - END DO -! - ENDDO - DO j=1, knon - i = ni(j) - t2m_w(i,nsrf)=yt2m_w(j) - q2m_w(i,nsrf)=yq2m_w(j) - ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman - ustar_w(i,nsrf)=yustar_w(j) - u10m_w(i,nsrf)=(yu10m_w(j) * uzon_w(j))/max(SQRT(uzon_w(j)**2+vmer_w(j)**2), smallestreal) - v10m_w(i,nsrf)=(yu10m_w(j) * vmer_w(j))/max(SQRT(uzon_w(j)**2+vmer_w(j)**2), smallestreal) -! - ustar(i,nsrf) = ustar_x(i,nsrf) + wake_s(i)*(ustar_w(i,nsrf)-ustar_x(i,nsrf)) - u10m(i,nsrf) = u10m_x(i,nsrf) + wake_s(i)*(u10m_w(i,nsrf)-u10m_x(i,nsrf)) - v10m(i,nsrf) = v10m_x(i,nsrf) + wake_s(i)*(v10m_w(i,nsrf)-v10m_x(i,nsrf)) -! - DO k = 1, 6 - n2mout_w(i,nsrf,k) = yn2mout_w(j,nsrf,k) - END DO -! - ENDDO -!!! - ENDIF ! (iflag_split .eq.0) -!!! - -! print*,'Dans pbl OK43' -!IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique -!IM Ajoute dependance type surface - IF (thermcep) THEN -!!! jyg le 07/02/2012 - IF (iflag_split .eq.0) THEN - DO j = 1, knon - i=ni(j) - zdelta1 = MAX(0.,SIGN(1., rtt-yt2m(j) )) - zx_qs1 = r2es * FOEEW(yt2m(j),zdelta1)/paprs(i,1) - zx_qs1 = MIN(0.5,zx_qs1) - zcor1 = 1./(1.-RETV*zx_qs1) - zx_qs1 = zx_qs1*zcor1 - - rh2m(i) = rh2m(i) + yq2m(j)/zx_qs1 * pctsrf(i,nsrf) - qsat2m(i) = qsat2m(i) + zx_qs1 * pctsrf(i,nsrf) - ENDDO - ELSE ! (iflag_split .eq.0) - DO j = 1, knon - i=ni(j) - zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_x(j) )) - zx_qs1 = r2es * FOEEW(yt2m_x(j),zdelta1)/paprs(i,1) - zx_qs1 = MIN(0.5,zx_qs1) - zcor1 = 1./(1.-RETV*zx_qs1) - zx_qs1 = zx_qs1*zcor1 - - rh2m_x(i) = rh2m_x(i) + yq2m_x(j)/zx_qs1 * pctsrf(i,nsrf) - qsat2m_x(i) = qsat2m_x(i) + zx_qs1 * pctsrf(i,nsrf) - ENDDO - DO j = 1, knon - i=ni(j) - zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_w(j) )) - zx_qs1 = r2es * FOEEW(yt2m_w(j),zdelta1)/paprs(i,1) - zx_qs1 = MIN(0.5,zx_qs1) - zcor1 = 1./(1.-RETV*zx_qs1) - zx_qs1 = zx_qs1*zcor1 - - rh2m_w(i) = rh2m_w(i) + yq2m_w(j)/zx_qs1 * pctsrf(i,nsrf) - qsat2m_w(i) = qsat2m_w(i) + zx_qs1 * pctsrf(i,nsrf) - ENDDO -!!! - ENDIF ! (iflag_split .eq.0) -!!! - ENDIF -! - IF (prt_level >=10) THEN - print *, 'T2m, q2m, RH2m ', & - t2m(1:knon,:), q2m(1:knon,:), rh2m(1:knon) - ENDIF - -! print*,'OK pbl 5' -! -!!! jyg le 07/02/2012 - IF (iflag_split .eq.0) THEN - CALL hbtm(knon, ypaprs, ypplay, & - yt2m,yt10m,yq2m,yq10m,yustar,ywstar, & - y_flux_t,y_flux_q,yu,yv,yt,yq, & - ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, & - ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl) - IF (prt_level >=10) THEN - print *,' Arg. de HBTM: yt2m ',yt2m(1:knon) - print *,' Arg. de HBTM: yt10m ',yt10m(1:knon) - print *,' Arg. de HBTM: yq2m ',yq2m(1:knon) - print *,' Arg. de HBTM: yq10m ',yq10m(1:knon) - print *,' Arg. de HBTM: yustar ',yustar(1:knon) - print *,' Arg. de HBTM: y_flux_t ',y_flux_t(1:knon,:) - print *,' Arg. de HBTM: y_flux_q ',y_flux_q(1:knon,:) - print *,' Arg. de HBTM: yu ',yu(1:knon,:) - print *,' Arg. de HBTM: yv ',yv(1:knon,:) - print *,' Arg. de HBTM: yt ',yt(1:knon,:) - print *,' Arg. de HBTM: yq ',yq(1:knon,:) - ENDIF - ELSE ! (iflag_split .eq.0) - CALL HBTM(knon, ypaprs, ypplay, & - yt2m_x,yt10m_x,yq2m_x,yq10m_x,yustar_x,ywstar_x, & - y_flux_t_x,y_flux_q_x,yu_x,yv_x,yt_x,yq_x, & - ypblh_x,ycapCL_x,yoliqCL_x,ycteiCL_x,ypblT_x, & - ytherm_x,ytrmb1_x,ytrmb2_x,ytrmb3_x,ylcl_x) - IF (prt_level >=10) THEN - print *,' Arg. de HBTM: yt2m_x ',yt2m_x(1:knon) - print *,' Arg. de HBTM: yt10m_x ',yt10m_x(1:knon) - print *,' Arg. de HBTM: yq2m_x ',yq2m_x(1:knon) - print *,' Arg. de HBTM: yq10m_x ',yq10m_x(1:knon) - print *,' Arg. de HBTM: yustar_x ',yustar_x(1:knon) - print *,' Arg. de HBTM: y_flux_t_x ',y_flux_t_x(1:knon,:) - print *,' Arg. de HBTM: y_flux_q_x ',y_flux_q_x(1:knon,:) - print *,' Arg. de HBTM: yu_x ',yu_x(1:knon,:) - print *,' Arg. de HBTM: yv_x ',yv_x(1:knon,:) - print *,' Arg. de HBTM: yt_x ',yt_x(1:knon,:) - print *,' Arg. de HBTM: yq_x ',yq_x(1:knon,:) - ENDIF - CALL HBTM(knon, ypaprs, ypplay, & - yt2m_w,yt10m_w,yq2m_w,yq10m_w,yustar_w,ywstar_w, & - y_flux_t_w,y_flux_q_w,yu_w,yv_w,yt_w,yq_w, & - ypblh_w,ycapCL_w,yoliqCL_w,ycteiCL_w,ypblT_w, & - ytherm_w,ytrmb1_w,ytrmb2_w,ytrmb3_w,ylcl_w) -!!! - ENDIF ! (iflag_split .eq.0) -!!! - -!!! jyg le 07/02/2012 - IF (iflag_split .eq.0) THEN -!!! - DO j=1, knon - i = ni(j) - pblh(i,nsrf) = ypblh(j) - wstar(i,nsrf) = ywstar(j) - plcl(i,nsrf) = ylcl(j) - capCL(i,nsrf) = ycapCL(j) - oliqCL(i,nsrf) = yoliqCL(j) - cteiCL(i,nsrf) = ycteiCL(j) - pblT(i,nsrf) = ypblT(j) - therm(i,nsrf) = ytherm(j) - trmb1(i,nsrf) = ytrmb1(j) - trmb2(i,nsrf) = ytrmb2(j) - trmb3(i,nsrf) = ytrmb3(j) - ENDDO - IF (prt_level >=10) THEN - print *, 'After HBTM: pblh ', pblh(1:knon,:) - print *, 'After HBTM: plcl ', plcl(1:knon,:) - print *, 'After HBTM: cteiCL ', cteiCL(1:knon,:) - ENDIF - ELSE !(iflag_split .eq.0) - DO j=1, knon - i = ni(j) - pblh_x(i,nsrf) = ypblh_x(j) - wstar_x(i,nsrf) = ywstar_x(j) - plcl_x(i,nsrf) = ylcl_x(j) - capCL_x(i,nsrf) = ycapCL_x(j) - oliqCL_x(i,nsrf) = yoliqCL_x(j) - cteiCL_x(i,nsrf) = ycteiCL_x(j) - pblT_x(i,nsrf) = ypblT_x(j) - therm_x(i,nsrf) = ytherm_x(j) - trmb1_x(i,nsrf) = ytrmb1_x(j) - trmb2_x(i,nsrf) = ytrmb2_x(j) - trmb3_x(i,nsrf) = ytrmb3_x(j) - ENDDO - IF (prt_level >=10) THEN - print *, 'After HBTM: pblh_x ', pblh_x(1:knon,:) - print *, 'After HBTM: plcl_x ', plcl_x(1:knon,:) - print *, 'After HBTM: cteiCL_x ', cteiCL_x(1:knon,:) - ENDIF - DO j=1, knon - i = ni(j) - pblh_w(i,nsrf) = ypblh_w(j) - wstar_w(i,nsrf) = ywstar_w(j) - plcl_w(i,nsrf) = ylcl_w(j) - capCL_w(i,nsrf) = ycapCL_w(j) - oliqCL_w(i,nsrf) = yoliqCL_w(j) - cteiCL_w(i,nsrf) = ycteiCL_w(j) - pblT_w(i,nsrf) = ypblT_w(j) - therm_w(i,nsrf) = ytherm_w(j) - trmb1_w(i,nsrf) = ytrmb1_w(j) - trmb2_w(i,nsrf) = ytrmb2_w(j) - trmb3_w(i,nsrf) = ytrmb3_w(j) - ENDDO - IF (prt_level >=10) THEN - print *, 'After HBTM: pblh_w ', pblh_w(1:knon,:) - print *, 'After HBTM: plcl_w ', plcl_w(1:knon,:) - print *, 'After HBTM: cteiCL_w ', cteiCL_w(1:knon,:) - ENDIF -!!! - ENDIF ! (iflag_split .eq.0) -!!! - -! print*,'OK pbl 6' -#else -! T2m not defined -! No calculation - PRINT*,' Warning !!! No T2m calculation. Output is set to zero.' -#endif - -!**************************************************************************************** -! 15) End of loop over different surfaces -! -!**************************************************************************************** - ENDDO loop_nbsrf -! -!---------------------------------------------------------------------------------------- -! Reset iflag_split -! - iflag_split=iflag_split_ref - -#ifdef ISO -#ifdef ISOVERIF -! write(*,*) 'pbl_surface tmp 3249: d_q,d_xt(iso_eau,2,1)=',d_q(2,1),d_xt(iso_eau,2,1) - IF (iso_eau >= 0) THEN - call iso_verif_egalite_vect2D( & - d_xt,d_q, & - 'pbl_surface_mod 1276',ntraciso,klon,klev) - ENDIF !IF (iso_eau >= 0) THEN -#endif -#endif - -!**************************************************************************************** -! 16) Calculate the mean value over all sub-surfaces for some variables -! -!**************************************************************************************** - - z0m(:,nbsrf+1) = 0.0 - z0h(:,nbsrf+1) = 0.0 - DO nsrf = 1, nbsrf - DO i = 1, klon - z0m(i,nbsrf+1) = z0m(i,nbsrf+1) + z0m(i,nsrf)*pctsrf(i,nsrf) - z0h(i,nbsrf+1) = z0h(i,nbsrf+1) + z0h(i,nsrf)*pctsrf(i,nsrf) - ENDDO - ENDDO - -! print*,'OK pbl 7' - zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0 - zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0 - zxfluxt_x(:,:) = 0.0 ; zxfluxq_x(:,:) = 0.0 - zxfluxu_x(:,:) = 0.0 ; zxfluxv_x(:,:) = 0.0 - zxfluxt_w(:,:) = 0.0 ; zxfluxq_w(:,:) = 0.0 - zxfluxu_w(:,:) = 0.0 ; zxfluxv_w(:,:) = 0.0 -#ifdef ISO - zxfluxxt(:,:,:) = 0.0 - zxfluxxt_x(:,:,:) = 0.0 - zxfluxxt_w(:,:,:) = 0.0 -#endif - - -!!! jyg le 07/02/2012 - IF (iflag_split .ge.1) THEN -!!! -!!! nrlmd & jyg les 02/05/2011, 05/02/2012 - - DO nsrf = 1, nbsrf - DO k = 1, klev - DO i = 1, klon - zxfluxt_x(i,k) = zxfluxt_x(i,k) + flux_t_x(i,k,nsrf) * pctsrf(i,nsrf) - zxfluxq_x(i,k) = zxfluxq_x(i,k) + flux_q_x(i,k,nsrf) * pctsrf(i,nsrf) - zxfluxu_x(i,k) = zxfluxu_x(i,k) + flux_u_x(i,k,nsrf) * pctsrf(i,nsrf) - zxfluxv_x(i,k) = zxfluxv_x(i,k) + flux_v_x(i,k,nsrf) * pctsrf(i,nsrf) -! - zxfluxt_w(i,k) = zxfluxt_w(i,k) + flux_t_w(i,k,nsrf) * pctsrf(i,nsrf) - zxfluxq_w(i,k) = zxfluxq_w(i,k) + flux_q_w(i,k,nsrf) * pctsrf(i,nsrf) - zxfluxu_w(i,k) = zxfluxu_w(i,k) + flux_u_w(i,k,nsrf) * pctsrf(i,nsrf) - zxfluxv_w(i,k) = zxfluxv_w(i,k) + flux_v_w(i,k,nsrf) * pctsrf(i,nsrf) -#ifdef ISO - DO ixt=1,ntraciso - zxfluxxt_x(ixt,i,k) = zxfluxxt_x(ixt,i,k) + flux_xt_x(ixt,i,k,nsrf) * pctsrf(i,nsrf) - zxfluxxt_w(ixt,i,k) = zxfluxxt_w(ixt,i,k) + flux_xt_w(ixt,i,k,nsrf) * pctsrf(i,nsrf) - ENDDO ! DO ixt=1,ntraciso -#endif - ENDDO - ENDDO - ENDDO - - DO i = 1, klon - zxsens_x(i) = - zxfluxt_x(i,1) - zxsens_w(i) = - zxfluxt_w(i,1) - ENDDO -!!! - ENDIF ! (iflag_split .ge.1) -!!! - - DO nsrf = 1, nbsrf - DO k = 1, klev - DO i = 1, klon - zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf(i,nsrf) - zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf(i,nsrf) - zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf) - zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf) -#ifdef ISO - DO ixt=1,niso - zxfluxxt(ixt,i,k) = zxfluxxt(ixt,i,k) + flux_xt(ixt,i,k,nsrf) * pctsrf(i,nsrf) - ENDDO ! DO ixt=1,niso -#endif - ENDDO - ENDDO - ENDDO - - DO i = 1, klon - zxsens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol - zxevap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol - fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i) - ENDDO - - ! if blowing snow - if (ok_bs) then - DO nsrf = 1, nbsrf - DO k = 1, klev - DO i = 1, klon - zxfluxqbs(i,k) = zxfluxqbs(i,k) + flux_qbs(i,k,nsrf) * pctsrf(i,nsrf) - ENDDO - ENDDO - ENDDO - - DO i = 1, klon - zxsnowerosion(i) = zxfluxqbs(i,1) ! blowings snow flux at the surface - END DO - endif - -#ifdef ISO - DO i = 1, klon - DO ixt=1,ntraciso - zxxtevap(ixt,i) = - zxfluxxt(ixt,i,1) - ENDDO - ENDDO -#endif - -!!! - -! -! Incrementer la temperature du sol -! - zxtsol(:) = 0.0 ; zxfluxlat(:) = 0.0 - zt2m(:) = 0.0 ; zq2m(:) = 0.0 ; zn2mout(:,:) = 0 - zustar(:)=0.0 ; zu10m(:) = 0.0 ; zv10m(:) = 0.0 - s_pblh(:) = 0.0 ; s_plcl(:) = 0.0 -!!! jyg le 07/02/2012 - s_pblh_x(:) = 0.0 ; s_plcl_x(:) = 0.0 - s_pblh_w(:) = 0.0 ; s_plcl_w(:) = 0.0 -!!! - s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0 - s_cteiCL(:) = 0.0; s_pblT(:) = 0.0 - s_therm(:) = 0.0 ; s_trmb1(:) = 0.0 - s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0 - wstar(:,is_ave)=0. - -! print*,'OK pbl 9' - -!!! nrlmd le 02/05/2011 - zxfluxlat_x(:) = 0.0 ; zxfluxlat_w(:) = 0.0 -!!! - - DO nsrf = 1, nbsrf - DO i = 1, klon - ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf) - - wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) & - + flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf) - - wfevap(i,nsrf) = evap(i,nsrf)*pctsrf(i,nsrf) - - zxtsol(i) = zxtsol(i) + ts(i,nsrf) * pctsrf(i,nsrf) - zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf(i,nsrf) - ENDDO - ENDDO -! -!al1 - -!!! jyg le 07/02/2012 - IF (iflag_split .eq.0) THEN - DO nsrf = 1, nbsrf - DO i = 1, klon - zt2m(i) = zt2m(i) + t2m(i,nsrf) * pctsrf(i,nsrf) - zq2m(i) = zq2m(i) + q2m(i,nsrf) * pctsrf(i,nsrf) -! - DO k = 1, 6 - zn2mout(i,k) = zn2mout(i,k) + n2mout(i,nsrf,k) * pctsrf(i,nsrf) - ENDDO -! - zustar(i) = zustar(i) + ustar(i,nsrf) * pctsrf(i,nsrf) - wstar(i,is_ave)=wstar(i,is_ave)+wstar(i,nsrf)*pctsrf(i,nsrf) - zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf(i,nsrf) - zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf(i,nsrf) - - s_pblh(i) = s_pblh(i) + pblh(i,nsrf) * pctsrf(i,nsrf) - s_plcl(i) = s_plcl(i) + plcl(i,nsrf) * pctsrf(i,nsrf) - s_capCL(i) = s_capCL(i) + capCL(i,nsrf) * pctsrf(i,nsrf) - s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf(i,nsrf) - s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf(i,nsrf) - s_pblT(i) = s_pblT(i) + pblT(i,nsrf) * pctsrf(i,nsrf) - s_therm(i) = s_therm(i) + therm(i,nsrf) * pctsrf(i,nsrf) - s_trmb1(i) = s_trmb1(i) + trmb1(i,nsrf) * pctsrf(i,nsrf) - s_trmb2(i) = s_trmb2(i) + trmb2(i,nsrf) * pctsrf(i,nsrf) - s_trmb3(i) = s_trmb3(i) + trmb3(i,nsrf) * pctsrf(i,nsrf) - ENDDO - ENDDO - ELSE !(iflag_split .eq.0) - DO nsrf = 1, nbsrf - DO i = 1, klon -!!! nrlmd le 02/05/2011 - zxfluxlat_x(i) = zxfluxlat_x(i) + fluxlat_x(i,nsrf) * pctsrf(i,nsrf) - zxfluxlat_w(i) = zxfluxlat_w(i) + fluxlat_w(i,nsrf) * pctsrf(i,nsrf) -!!! -!!! jyg le 08/02/2012 -!! Pour le moment, on sort les valeurs dans (x) et (w) de pblh et de plcl ; -!! pour zt2m, on fait la moyenne surfacique sur les sous-surfaces ; -!! pour qsat2m, on fait la moyenne surfacique sur (x) et (w) ; -!! pour les autres variables, on sort les valeurs de la region (x). - zt2m(i) = zt2m(i) + (t2m_x(i,nsrf)+wake_s(i)*(t2m_w(i,nsrf)-t2m_x(i,nsrf))) * pctsrf(i,nsrf) - zq2m(i) = zq2m(i) + q2m_x(i,nsrf) * pctsrf(i,nsrf) -! - DO k = 1, 6 - zn2mout(i,k) = zn2mout(i,k) + n2mout_x(i,nsrf,k) * pctsrf(i,nsrf) - ENDDO -! - zustar(i) = zustar(i) + ustar_x(i,nsrf) * pctsrf(i,nsrf) - wstar(i,is_ave)=wstar(i,is_ave)+wstar_x(i,nsrf)*pctsrf(i,nsrf) - zu10m(i) = zu10m(i) + u10m_x(i,nsrf) * pctsrf(i,nsrf) - zv10m(i) = zv10m(i) + v10m_x(i,nsrf) * pctsrf(i,nsrf) -! - s_pblh(i) = s_pblh(i) + pblh_x(i,nsrf) * pctsrf(i,nsrf) - s_pblh_x(i) = s_pblh_x(i) + pblh_x(i,nsrf) * pctsrf(i,nsrf) - s_pblh_w(i) = s_pblh_w(i) + pblh_w(i,nsrf) * pctsrf(i,nsrf) -! - s_plcl(i) = s_plcl(i) + plcl_x(i,nsrf) * pctsrf(i,nsrf) - s_plcl_x(i) = s_plcl_x(i) + plcl_x(i,nsrf) * pctsrf(i,nsrf) - s_plcl_w(i) = s_plcl_w(i) + plcl_w(i,nsrf) * pctsrf(i,nsrf) -! - s_capCL(i) = s_capCL(i) + capCL_x(i,nsrf) * pctsrf(i,nsrf) - s_oliqCL(i) = s_oliqCL(i) + oliqCL_x(i,nsrf)* pctsrf(i,nsrf) - s_cteiCL(i) = s_cteiCL(i) + cteiCL_x(i,nsrf)* pctsrf(i,nsrf) - s_pblT(i) = s_pblT(i) + pblT_x(i,nsrf) * pctsrf(i,nsrf) - s_therm(i) = s_therm(i) + therm_x(i,nsrf) * pctsrf(i,nsrf) - s_trmb1(i) = s_trmb1(i) + trmb1_x(i,nsrf) * pctsrf(i,nsrf) - s_trmb2(i) = s_trmb2(i) + trmb2_x(i,nsrf) * pctsrf(i,nsrf) - s_trmb3(i) = s_trmb3(i) + trmb3_x(i,nsrf) * pctsrf(i,nsrf) - ENDDO - ENDDO - DO i = 1, klon - qsat2m(i)= qsat2m_x(i)+ wake_s(i)*(qsat2m_x(i)-qsat2m_w(i)) - ENDDO -!!! - ENDIF ! (iflag_split .eq.0) -!!! - - IF (check) THEN - amn=MIN(ts(1,is_ter),1000.) - amx=MAX(ts(1,is_ter),-1000.) - DO i=2, klon - amn=MIN(ts(i,is_ter),amn) - amx=MAX(ts(i,is_ter),amx) - ENDDO - PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx - ENDIF - -!jg ? -!!$! -!!$! If a sub-surface does not exsist for a grid point, the mean value for all -!!$! sub-surfaces is distributed. -!!$! -!!$ DO nsrf = 1, nbsrf -!!$ DO i = 1, klon -!!$ IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN -!!$ ts(i,nsrf) = zxtsol(i) -!!$ t2m(i,nsrf) = zt2m(i) -!!$ q2m(i,nsrf) = zq2m(i) -!!$ u10m(i,nsrf) = zu10m(i) -!!$ v10m(i,nsrf) = zv10m(i) -!!$ -!!$! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour -!!$ pblh(i,nsrf) = s_pblh(i) -!!$ plcl(i,nsrf) = s_plcl(i) -!!$ capCL(i,nsrf) = s_capCL(i) -!!$ oliqCL(i,nsrf) = s_oliqCL(i) -!!$ cteiCL(i,nsrf) = s_cteiCL(i) -!!$ pblT(i,nsrf) = s_pblT(i) -!!$ therm(i,nsrf) = s_therm(i) -!!$ trmb1(i,nsrf) = s_trmb1(i) -!!$ trmb2(i,nsrf) = s_trmb2(i) -!!$ trmb3(i,nsrf) = s_trmb3(i) -!!$ ENDIF -!!$ ENDDO -!!$ ENDDO - - - DO i = 1, klon - fder(i) = - 4.0*RSIGMA*zxtsol(i)**3 - ENDDO - - zxqsurf(:) = 0.0 - zxsnow(:) = 0.0 -#ifdef ISO - zxxtsnow(:,:) = 0.0 -#endif - - DO nsrf = 1, nbsrf - DO i = 1, klon - zxqsurf(i) = zxqsurf(i) + MAX(qsurf(i,nsrf),0.0) * pctsrf(i,nsrf) - zxsnow(i) = zxsnow(i) + snow(i,nsrf) * pctsrf(i,nsrf) -#ifdef ISO - DO ixt=1,niso - zxxtsnow(ixt,i) = zxxtsnow(ixt,i) + xtsnow(ixt,i,nsrf) * pctsrf(i,nsrf) - ENDDO ! DO ixt=1,niso -#endif - ENDDO - ENDDO - -! Premier niveau de vent sortie dans physiq.F - zu1(:) = u(:,1) - zv1(:) = v(:,1) - - END SUBROUTINE pbl_surface - -#endif - - - - !