MODULE pbl_surface_subsrf_mod ! ! Planetary Boundary Layer and Surface module ! ! This module manages the calculation of turbulent diffusion in the boundary layer ! and all interactions towards the differents sub-surfaces. ! ! USE pbl_surface_data USE dimphy USE mod_phys_lmdz_para, ONLY : mpi_size USE mod_grid_phy_lmdz, ONLY : klon_glo USE ioipsl USE surface_data, ONLY : type_ocean, ok_veget, landice_opt, iflag_leads USE surf_land_mod, ONLY : surf_land USE surf_landice_mod, ONLY : surf_landice USE surf_ocean_mod, ONLY : surf_ocean USE surf_seaice_mod, ONLY : surf_seaice USE cpl_mod, ONLY : gath2cpl USE climb_hq_mod, ONLY : climb_hq_down, climb_hq_up USE climb_qbs_mod, ONLY : climb_qbs_down, climb_qbs_up USE climb_wind_mod, ONLY : climb_wind_down, climb_wind_up USE coef_diff_turb_mod, ONLY : coef_diff_turb USE lmdz_call_atke, ONLY : call_atke USE ioipsl_getin_p_mod, ONLY : getin_p USE cdrag_mod USE stdlevvar_mod USE wx_pbl_var_mod, ONLY : wx_pbl_init, wx_pbl_final, & wx_pbl_prelim_0, wx_pbl_prelim_beta USE wx_pbl_mod, ONLY : wx_pbl0_merge, wx_pbl_split, wx_pbl_dts_merge, & wx_pbl_check, wx_pbl_dts_check, wx_evappot use config_ocean_skin_m, only: activate_ocean_skin #ifdef ISO USE infotrac_phy, ONLY: niso,ntraciso=>ntiso #endif IMPLICIT NONE PRIVATE PUBLIC :: pbl_surface_subsrf CONTAINS SUBROUTINE pbl_surface_subsrf( nsrf, knon, ni, & dtime, date0, itap, jour, & debut, lafin, & rlon, rlat, rugoro, rmu0, & lwdown_m, pphi, cldt, & rain_f, snow_f, bs_f, & gustiness, & t, q, qbs, u, v, & wake_dlt, wake_dlq, & wake_cstar, wake_s, & pplay, paprs, pctsrf, & ts,SFRWL, alb_dir, alb_dif,ustar, u10m, v10m,wstar, & cdragh, cdragm, & beta, & icesub_ice, icemelt_ice, alb3_lic, runoff, snowhgt, qsnow, to_ice, sissnow, & qsat2m, & d_t, d_q, d_qbs, d_u, d_v, d_t_diss, & d_t_w, d_q_w, & d_t_x, d_q_x, & delta_tsurf,wake_dens,cdragh_x,cdragh_w, & cdragm_x,cdragm_w,kh,kh_x,kh_w, & zcoefh, zcoefm, slab_wfbils, & qsol, s_pblh, & s_pblh_x, s_pblh_w, & delta_qsurf, & rh2m, & z0m, z0h, agesno, sollw, solsw, & d_ts, evap, fluxlat, t2m, & flux_t, flux_u, flux_v, & dflux_t, dflux_q, & q2m, flux_q, flux_qbs, tke_x, eps_x, & wake_dltke, & treedrg,hice ,tice, bilg_cumul, & fcds, fcdi, dh_basal_growth, dh_basal_melt, & dh_top_melt, dh_snow2sic, & dtice_melt, dtice_snow2sic , & 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 , n2mout, n2mout_x, n2mout_w, d_u_x, d_u_w, d_v_x, d_v_w, windsp, t2m_x, & q2m_x, rh2m_x, qsat2m_x, u10m_x, v10m_x, ustar_x, wstar_x, pblh_x, plcl_x, capCL_x, & oliqCL_x, cteiCL_x, pblt_x, therm_x, trmb1_x, trmb2_x, trmb3_x, t2m_w, q2m_w, rh2m_w, & qsat2m_w, u10m_w, v10m_w, ustar_w, wstar_w, pblh_w, plcl_w, capCL_w, oliqCL_w, cteiCL_w,& pblt_w, therm_w, trmb1_w, trmb2_w, trmb3_w, pblh, plcl, capCL, oliqCL, cteiCL, pblT, & therm, trmb1, trmb2, trmb3, alb, snowerosion, iflag_split_ref, & delp, d_t_diss_x, d_t_diss_w, flux_t_x, flux_q_x, flux_t_w, flux_q_w,& flux_u_x, flux_v_x, flux_u_w, flux_v_w, fluxlat_x, fluxlat_w) !$gpum horizontal knon klon !**************************************************************************************** ! 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 USE carbon_cycle_mod, ONLY : co2_send, nbcf_out, fields_out, yfields_out use hbtm_mod, only: hbtm USE indice_sol_mod USE mod_grid_phy_lmdz, ONLY : 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 !AM use phys_output_var_mod, only: tkt, tks, taur, sss use lmdz_blosno_ini, only : zeta_bs USE dimsoil_mod_h, ONLY: nsoilmx USE surf_param_mod, ONLY: eff_surf_param !AM use wxios_mod, ONLY: missing_val_xios => missing_val, using_xios USE netcdf, only: missing_val_netcdf => nf90_fill_real USE yomcst_mod_h USE mod_phys_lmdz_para, ONLY : is_master USE cdrag_mod, ONLY : cdrag USE freinage_mod, ONLY : freinage USE yamada_c_mod, ONLY : yamada_c IMPLICIT NONE INCLUDE "FCTTRE.h" !**************************************************************************************** INTEGER, INTENT(IN) :: nsrf ! indice current subsurface INTEGER, INTENT(IN) :: knon ! number of compressed points for current subsurface INTEGER, INTENT(IN) :: ni(knon)! index for compressed current sub-surface 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,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 REAL, DIMENSION(klon), INTENT(IN) :: lwdown_m ! downward longwave radiation at mean s REAL, DIMENSION(klon), INTENT(IN) :: gustiness ! gustiness REAL, DIMENSION(klon,klev), INTENT(IN) :: pphi ! geopotential (m2/s2) 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 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 !**************************************************************************************** REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: beta ! Aridity factor REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: ts ! temperature at surface (K) REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: delta_tsurf !surface temperature difference between REAL, DIMENSIOn(6),intent(in) :: SFRWL REAL, DIMENSION(klon, nsw, nbsrf), INTENT(INOUT) :: alb_dir,alb_dif 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 REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke_x 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(INOUT) :: eps_x ! TKE dissipation rate REAL, DIMENSION(klon), INTENT(INOUT) :: cdragh ! drag coefficient for T and Q REAL, DIMENSION(klon), INTENT(INOUT) :: cdragm ! drag coefficient for wind REAL, DIMENSION(klon), INTENT(INOUT) :: alb3_lic REAL, DIMENSION(klon), INTENT(INOUT) :: icesub_ice ! ice (no snow!) sublimation flux over ice sheet and sea ice REAL, DIMENSION(klon), INTENT(INOUT) :: icemelt_ice ! ice (no snow!) melting flux over ice sheet and sea ice REAL, DIMENSION(klon,klev), INTENT(INOUT) :: d_t_w ! ! REAL, DIMENSION(klon,klev), INTENT(INOUT) :: d_q_w ! ! Tendances dans les poches REAL, DIMENSION(klon,klev), INTENT(INOUT) :: d_t_x ! ! REAL, DIMENSION(klon,klev), INTENT(INOUT) :: d_q_x ! ! Tendances hors des poches REAL, DIMENSION(klon), INTENT(INOUT) :: qsat2m REAL, DIMENSION(klon, klev), INTENT(INOUT) :: d_t ! change in temperature REAL, DIMENSION(klon, klev), INTENT(INOUT) :: d_t_diss ! change in temperature REAL, DIMENSION(klon, klev), INTENT(INOUT) :: d_q ! change in water vapour REAL, DIMENSION(klon, klev), INTENT(INOUT) :: d_u ! change in u speed REAL, DIMENSION(klon, klev), INTENT(INOUT) :: d_v ! change in v speed REAL, DIMENSION(klon, klev), INTENT(INOUT) :: d_qbs ! change in blowing snow specific content REAL, INTENT(INOUT):: 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(INOUT):: 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(INOUT) :: zxxtevap ! water vapour flux at surface, positiv upwards REAL, DIMENSION(ntraciso,klon, klev), INTENT(INOUT) :: d_xt ! change in water vapour REAL, DIMENSION(klon), INTENT(INOUT) :: runoff_diag REAL, DIMENSION(niso,klon), INTENT(INOUT) :: xtrunoff_diag REAL, DIMENSION(ntraciso,klon,klev), INTENT(INOUT) :: d_xt_w REAL, DIMENSION(ntraciso,klon,klev), INTENT(INOUT) :: d_xt_x #endif ! Output only for diagnostics REAL, DIMENSION(klon), INTENT(INOUT) :: cdragh_x REAL, DIMENSION(klon), INTENT(INOUT) :: cdragh_w REAL, DIMENSION(klon), INTENT(INOUT) :: cdragm_x REAL, DIMENSION(klon), INTENT(INOUT) :: cdragm_w REAL, DIMENSION(klon), INTENT(INOUT) :: kh REAL, DIMENSION(klon), INTENT(INOUT) :: kh_x REAL, DIMENSION(klon), INTENT(INOUT) :: kh_w REAL, DIMENSION(klon), INTENT(INOUT) :: slab_wfbils! heat balance at surface only for slab at ocean points REAL, DIMENSION(klon), INTENT(INOUT) :: qsol ! water height in the soil (mm) REAL, DIMENSION(klon), INTENT(INOUT) :: s_pblh ! height of the planetary boundary layer(HPBL) REAL, DIMENSION(klon), INTENT(INOUT) :: s_pblh_x ! height of the PBL in the off-wake region REAL, DIMENSION(klon), INTENT(INOUT) :: s_pblh_w ! height of the PBL in the wake region REAL, DIMENSION(klon), INTENT(INOUT) :: delta_qsurf! humidity difference at surface, mean for each grid point REAL, DIMENSION(klon), INTENT(INOUT) :: rh2m ! relative humidity at 2m 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(IN) :: solsw ! net shortwave radiation at surface REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: sollw ! net longwave radiation at surface REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: d_ts ! change in temperature at surface REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: evap ! evaporation at surface REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: fluxlat ! latent flux REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: t2m ! temperature at 2 meter height REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: flux_t ! sensible heat flux (CpT) J/m**2/s (W/m**2) ! positve orientation downwards REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: flux_u ! u wind tension (kg m/s)/(m**2 s) or Pascal REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: flux_v ! v wind tension (kg m/s)/(m**2 s) or Pascal 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(INOUT) :: xtsol ! water height in the soil (mm) REAL, DIMENSION(ntraciso,klon, nbsrf) :: xtevap ! evaporation at surface REAL, DIMENSION(klon), INTENT(INOUT) :: h1_diag ! just diagnostic, not useful #endif ! Output not needed REAL, DIMENSION(klon), INTENT(INOUT) :: dflux_t ! change of sensible heat flux REAL, DIMENSION(klon), INTENT(INOUT) :: dflux_q ! change of water vapour flux REAL, DIMENSION(klon, nbsrf),INTENT(INOUT) :: q2m ! water vapour at 2 meter height REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: flux_q ! water vapour flux(latent flux) (kg/m**2/s) REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: flux_qbs ! blowind snow vertical flux (kg/m**2 #ifdef ISO REAL, DIMENSION(ntraciso,klon), INTENT(INOUT) :: dflux_xt ! change of water vapour flux REAL, DIMENSION(niso,klon), INTENT(INOUT) :: zxxtsnow ! snow at surface, mean for each grid point REAL, DIMENSION(ntraciso,klon, klev), INTENT(INOUT) :: zxfluxxt ! water vapour flux, mean for each grid point REAL, DIMENSION(ntraciso,klon, klev, nbsrf), INTENT(INOUT) :: flux_xt ! water vapour flux(latent flux) (kg/m**2/s) #endif REAL, DIMENSION(klon), INTENT(INOUT) :: qsnow ! snow water content REAL, DIMENSION(klon), INTENT(INOUT) :: snowhgt ! snow height REAL, DIMENSION(klon), INTENT(INOUT) :: to_ice ! snow passed to ice REAL, DIMENSION(klon), INTENT(INOUT) :: sissnow ! snow in snow model REAL, DIMENSION(klon), INTENT(INOUT) :: runoff ! runoff on land ice 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 ! variables temporaires en "klon" (nom compressée) passée en argument pour les sous-surface INTEGER, DIMENSION(klon, nbsrf, 6), INTENT(INOUT) :: n2mout INTEGER, DIMENSION(klon, nbsrf, 6), INTENT(INOUT) :: n2mout_x INTEGER, DIMENSION(klon, nbsrf, 6), INTENT(INOUT) :: n2mout_w REAL, DIMENSION(klon, klev), INTENT(INOUT) :: d_u_x REAL, DIMENSION(klon, klev), INTENT(INOUT) :: d_u_w REAL, DIMENSION(klon, klev), INTENT(INOUT) :: d_v_x REAL, DIMENSION(klon, klev), INTENT(INOUT) :: d_v_w REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: windsp REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: t2m_x REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: q2m_x REAL, DIMENSION(klon), INTENT(INOUT) :: rh2m_x REAL, DIMENSION(klon), INTENT(INOUT) :: qsat2m_x REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: u10m_x REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: v10m_x REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: ustar_x REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: wstar_x REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: pblh_x REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: plcl_x REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: capCL_x REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: oliqCL_x REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: cteiCL_x REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: pblt_x REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: therm_x REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: trmb1_x REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: trmb2_x REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: trmb3_x REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: t2m_w REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: q2m_w REAL, DIMENSION(klon), INTENT(INOUT) :: rh2m_w REAL, DIMENSION(klon), INTENT(INOUT) :: qsat2m_w REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: u10m_w REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: v10m_w REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: ustar_w REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: wstar_w ! REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: pblh_w REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: plcl_w REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: capCL_w REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: oliqCL_w REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: cteiCL_w REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: pblt_w REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: therm_w REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: trmb1_w REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: trmb2_w REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: trmb3_w ! REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pblh ! height of the planetary boundary layer REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: plcl ! condensation level REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: capCL REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: oliqCL REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: cteiCL REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pblT REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: therm REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: trmb1 ! deep cape REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: trmb2 ! inhibition REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: trmb3 ! point Omega REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: alb ! mean albedo for whole SW interval REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: snowerosion INTEGER, INTENT(INOUT) :: iflag_split_ref ! ! Other local variables !**************************************************************************************** INTEGER :: n INTEGER :: iflag_split INTEGER :: i, k INTEGER :: j REAL, DIMENSION(knon) :: r_co2_ppm ! taux CO2 atmosphere REAL :: yt1_new REAL, DIMENSION(knon) :: yts, yz0m, yz0h, ypct REAL, DIMENSION(knon) :: yz0h_old REAL, DIMENSION(knon) :: yalb,yalb_vis REAL, DIMENSION(knon) :: yt1, yq1, yu1, yv1, yqbs1 REAL, DIMENSION(knon) :: yqa REAL, DIMENSION(knon) :: ysnow, yqsurf, yagesno, yqsol REAL, DIMENSION(knon) :: yrain_f, ysnow_f, ybs_f #ifdef ISO REAL, DIMENSION(ntraciso,knon) :: yxt1 REAL, DIMENSION(niso,knon) :: yxtsnow, yxtsol REAL, DIMENSION(ntraciso,knon) :: yxtrain_f, yxtsnow_f REAL, DIMENSION(knon) :: yrunoff_diag REAL, DIMENSION(niso,knon) :: yxtrunoff_diag REAL, DIMENSION(niso,knon) :: yRland_ice #endif REAL, DIMENSION(knon) :: ysolsw, ysollw REAL, DIMENSION(knon) :: yfder REAL, DIMENSION(knon) :: yrugoro REAL, DIMENSION(knon) :: yfluxlat REAL, DIMENSION(knon) :: yfluxbs REAL, DIMENSION(knon) :: y_d_ts REAL, DIMENSION(knon) :: y_flux_t1, y_flux_q1 REAL, DIMENSION(knon) :: y_dflux_t, y_dflux_q #ifdef ISO REAL, DIMENSION(ntraciso,knon) :: y_flux_xt1 REAL, DIMENSION(ntraciso,knon) :: y_dflux_xt #endif REAL, DIMENSION(knon) :: y_flux_u1, y_flux_v1 REAL, DIMENSION(knon) :: y_flux_bs, y_flux0 REAL, DIMENSION(knon) :: yt2m, yq2m, yu10m INTEGER, DIMENSION(knon, nbsrf, 6) :: yn2mout, yn2mout_x, yn2mout_w REAL, DIMENSION(knon) :: yustar REAL, DIMENSION(knon) :: ywstar REAL, DIMENSION(knon) :: ywindsp REAL, DIMENSION(knon) :: yt10m, yq10m REAL, DIMENSION(knon) :: ypblh REAL, DIMENSION(knon) :: ylcl REAL, DIMENSION(knon) :: ycapCL REAL, DIMENSION(knon) :: yoliqCL REAL, DIMENSION(knon) :: ycteiCL REAL, DIMENSION(knon) :: ypblT REAL, DIMENSION(knon) :: ytherm REAL, DIMENSION(knon) :: ytrmb1 REAL, DIMENSION(knon) :: ytrmb2 REAL, DIMENSION(knon) :: ytrmb3 REAL, DIMENSION(knon) :: yt2m_x REAL, DIMENSION(knon) :: yq2m_x REAL, DIMENSION(knon) :: yt10m_x REAL, DIMENSION(knon) :: yq10m_x REAL, DIMENSION(knon) :: yu10m_x REAL, DIMENSION(knon) :: yustar_x REAL, DIMENSION(knon) :: ywstar_x ! REAL, DIMENSION(knon) :: ypblh_x REAL, DIMENSION(knon) :: ylcl_x REAL, DIMENSION(knon) :: ycapCL_x REAL, DIMENSION(knon) :: yoliqCL_x REAL, DIMENSION(knon) :: ycteiCL_x REAL, DIMENSION(knon) :: ypblt_x REAL, DIMENSION(knon) :: ytherm_x REAL, DIMENSION(knon) :: ytrmb1_x REAL, DIMENSION(knon) :: ytrmb2_x REAL, DIMENSION(knon) :: ytrmb3_x REAL, DIMENSION(knon) :: uzon, vmer REAL, DIMENSION(knon) :: tair1, qair1, tairsol REAL, DIMENSION(knon) :: psfce, patm REAL, DIMENSION(knon) :: qairsol, zgeo1, speed, zri1, pref !speed, zri1, pref, added by Fuxing WANG, 04/03/2015 REAL, DIMENSION(knon) :: yz0h_oupas REAL, DIMENSION(knon) :: yfluxsens REAL, DIMENSION(knon) :: AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0 REAL, DIMENSION(knon) :: AcoefH, AcoefQ, BcoefH, BcoefQ #ifdef ISO REAL, DIMENSION(ntraciso,knon) :: AcoefXT, BcoefXT #endif REAL, DIMENSION(knon) :: AcoefU, AcoefV, BcoefU, BcoefV REAL, DIMENSION(knon) :: AcoefQBS, BcoefQBS REAL, DIMENSION(knon) :: ypsref REAL, DIMENSION(knon) :: yevap, yevap_pot, ytsurf_new, yalb3_new, yicesub, yicemelt REAL, DIMENSION(knon,nsw) :: yalb_dir_new, yalb_dif_new REAL, DIMENSION(knon,klev) :: y_d_t, y_d_q, y_d_t_diss, y_d_qbs REAL, DIMENSION(knon,klev) :: y_d_u, y_d_v REAL, DIMENSION(knon,klev) :: y_flux_t, y_flux_q, y_flux_qbs REAL, DIMENSION(knon,klev) :: y_flux_u, y_flux_v REAL, DIMENSION(knon,klev) :: ycoefh,ycoefm,ycoefq,ycoefqbs REAL, DIMENSION(knon) :: ycdragh, ycdragq, ycdragm REAL, DIMENSION(knon,klev) :: yu, yv REAL, DIMENSION(knon,klev) :: yt, yq, yqbs #ifdef ISO REAL, DIMENSION(ntraciso,knon) :: yxtevap REAL, DIMENSION(ntraciso,knon,klev) :: y_d_xt REAL, DIMENSION(ntraciso,knon,klev) :: y_flux_xt REAL, DIMENSION(ntraciso,knon,klev) :: yxt #endif REAL, DIMENSION(knon,klev) :: ypplay, ydelp REAL, DIMENSION(klon,klev),INTENT(IN) :: delp REAL, DIMENSION(knon,klev+1) :: ypaprs REAL, DIMENSION(knon,klev+1) :: ytke, yeps REAL, DIMENSION(knon,nsoilmx) :: ytsoil REAL, DIMENSION(knon,nvm_lmdz) :: yveget REAL, DIMENSION(knon,nvm_lmdz) :: ylai REAL, DIMENSION(knon,nvm_lmdz) :: yheight REAL, DIMENSION(knon,klev) :: y_d_u_frein REAL, DIMENSION(knon,klev) :: y_d_v_frein REAL, DIMENSION(knon,klev) :: y_treedrg 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. REAL, DIMENSION(knon) :: ywake_s, ywake_cstar, ywake_dens REAL, DIMENSION(knon,klev+1) :: ytke_x, ytke_w, yeps_x, yeps_w REAL, DIMENSION(knon,klev+1) :: ywake_dltke REAL, DIMENSION(knon,klev) :: yu_x, yv_x, yu_w, yv_w REAL, DIMENSION(knon,klev) :: yt_x, yq_x, yt_w, yq_w REAL, DIMENSION(knon,klev) :: ycoefh_x, ycoefm_x, ycoefh_w, ycoefm_w REAL, DIMENSION(knon,klev) :: ycoefq_x, ycoefq_w REAL, DIMENSION(knon) :: ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w REAL, DIMENSION(knon) :: ycdragm_x, ycdragm_w REAL, DIMENSION(knon) :: AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x REAL, DIMENSION(knon) :: AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w REAL, DIMENSION(knon) :: AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x REAL, DIMENSION(knon) :: AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w REAL, DIMENSION(knon) :: y_flux_t1_x, y_flux_q1_x, y_flux_t1_w, y_flux_q1_w REAL, DIMENSION(knon) :: y_flux_u1_x, y_flux_v1_x, y_flux_u1_w, y_flux_v1_w REAL, DIMENSION(knon,klev) :: y_flux_t_x, y_flux_q_x, y_flux_t_w, y_flux_q_w REAL, DIMENSION(knon,klev) :: y_flux_u_x, y_flux_v_x, y_flux_u_w, y_flux_v_w REAL, DIMENSION(knon) :: yfluxlat_x, yfluxlat_w REAL, DIMENSION(knon,klev) :: y_d_t_x, y_d_q_x, y_d_t_w, y_d_q_w REAL, DIMENSION(knon,klev) :: y_d_t_diss_x, y_d_t_diss_w REAL, DIMENSION(klon,klev), INTENT(INOUT) :: d_t_diss_x, d_t_diss_w REAL, DIMENSION(knon,klev) :: y_d_u_x, y_d_v_x, y_d_u_w, y_d_v_w REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: flux_t_x, flux_q_x, flux_t_w, flux_q_w REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: flux_u_x, flux_v_x, flux_u_w, flux_v_w REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: fluxlat_x, fluxlat_w REAL, DIMENSION(knon) :: ybeta REAL, DIMENSION(knon) :: ybeta_prev REAL, DIMENSION(knon,klev) :: CcoefH, CcoefQ, DcoefH, DcoefQ REAL, DIMENSION(knon,klev) :: CcoefU, CcoefV, DcoefU, DcoefV REAL, DIMENSION(knon,klev) :: CcoefQBS, DcoefQBS REAL, DIMENSION(knon,klev) :: CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x REAL, DIMENSION(knon,klev) :: CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w REAL, DIMENSION(knon,klev) :: CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x REAL, DIMENSION(knon,klev) :: CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w REAL, DIMENSION(knon,klev) :: Kcoef_hq, Kcoef_m, gama_h, gama_q REAL, DIMENSION(knon,klev) :: gama_qbs, Kcoef_qbs REAL, DIMENSION(knon,klev) :: Kcoef_hq_x, Kcoef_m_x, gama_h_x, gama_q_x REAL, DIMENSION(knon,klev) :: Kcoef_hq_w, Kcoef_m_w, gama_h_w, gama_q_w REAL, DIMENSION(knon) :: alf_1, alf_2, alf_1_x, alf_2_x, alf_1_w, alf_2_w #ifdef ISO REAL, DIMENSION(ntraciso,knon,klev) :: yxt_x, yxt_w REAL, DIMENSION(ntraciso,knon) :: y_flux_xt1_x , y_flux_xt1_w REAL, DIMENSION(ntraciso,knon,klev) :: y_flux_xt_x,y_d_xt_x REAL, DIMENSION(ntraciso,knon,klev) :: y_flux_xt_w,y_d_xt_w REAL, DIMENSION(ntraciso,klon,klev),INTENT(INOUT) :: zxfluxxt_w, zxfluxxt_x REAL, DIMENSION(ntraciso,klon,klev,nbsrf), INTENT(INOUT) :: flux_xt_x, flux_xt_w REAL, DIMENSION(ntraciso,knon) :: AcoefXT_x, BcoefXT_x REAL, DIMENSION(ntraciso,knon) :: AcoefXT_w, BcoefXT_w REAL, DIMENSION(ntraciso,knon,klev) :: CcoefXT, DcoefXT REAL, DIMENSION(ntraciso,knon,klev) :: CcoefXT_x, DcoefXT_x REAL, DIMENSION(ntraciso,knon,klev) :: CcoefXT_w, DcoefXT_w REAL, DIMENSION(ntraciso,knon,klev) :: gama_xt,gama_xt_x,gama_xt_w #endif REAL, DIMENSION(knon) :: yt2m_w REAL, DIMENSION(knon) :: yq2m_w REAL, DIMENSION(knon) :: yt10m_w REAL, DIMENSION(knon) :: yq10m_w REAL, DIMENSION(knon) :: yu10m_w REAL, DIMENSION(knon) :: yustar_w REAL, DIMENSION(knon) :: ywstar_w ! REAL, DIMENSION(knon) :: ypblh_w REAL, DIMENSION(knon) :: ylcl_w REAL, DIMENSION(knon) :: ycapCL_w REAL, DIMENSION(knon) :: yoliqCL_w REAL, DIMENSION(knon) :: ycteiCL_w REAL, DIMENSION(knon) :: ypblt_w REAL, DIMENSION(knon) :: ytherm_w REAL, DIMENSION(knon) :: ytrmb1_w REAL, DIMENSION(knon) :: ytrmb2_w REAL, DIMENSION(knon) :: ytrmb3_w ! REAL, DIMENSION(knon) :: 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(knon) :: zgeo1_x, tair1_x, qair1_x, tairsol_x ! REAL, DIMENSION(knon) :: 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(knon) :: zgeo1_w, tair1_w, qair1_w, tairsol_w REAL, DIMENSION(knon) :: yus0, yvs0 ! REAL, DIMENSION(knon) :: y_delta_flux_t1 REAL, DIMENSION(knon) :: y_delta_tsurf, y_delta_tsurf_new REAL, DIMENSION(knon) :: delta_coef, tau_eq REAL, DIMENSION(knon) :: HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn REAL, DIMENSION(knon) :: phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0 REAL, DIMENSION(knon) :: y_delta_qsurf REAL, DIMENSION(knon) :: y_delta_qsats REAL, DIMENSION(knon) :: yg_T, yg_Q REAL, DIMENSION(knon) :: yGamma_dTs_phiT, yGamma_dQs_phiQ REAL, DIMENSION(knon) :: ydTs_ins, ydqs_ins ! REAL, PARAMETER :: facteur = 2. / 1.772 ! ( == 2. / SQRT(3.14)) REAL, PARAMETER :: inertia=2000. REAL, DIMENSION(knon) :: ydtsurf_th REAL, DIMENSION(knon) :: Kech_h ! Coefficient d'echange pour l'energie REAL, DIMENSION(knon) :: Kech_h_x, Kech_h_w REAL, DIMENSION(knon) :: yts_x, yts_w REAL, DIMENSION(knon) :: yqsurf_x, yqsurf_w REAL :: fact_cdrag REAL :: z1lay REAL :: vent REAL, DIMENSION(knon) :: ylwdown ! jg : temporary (ysollwdown) REAL, DIMENSION(knon) :: ygustiness ! jg : temporary (ysollwdown) REAL :: zx_qs1, zcor1, zdelta1 REAL, DIMENSION(knon) :: ytoice REAL, DIMENSION(knon) :: ysnowhgt, yqsnow, ysissnow, yrunoff REAL, DIMENSION(knon) :: yzmea REAL, DIMENSION(knon) :: yzsig REAL, DIMENSION(knon) :: ycldt REAL, DIMENSION(knon) :: yrmu0 REAL, DIMENSION(knon) :: yri0 REAL, DIMENSION(knon):: 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 REAL, DIMENSION(knon,klev) :: ytheta REAL, DIMENSION(knon,klev) :: ypphii REAL, DIMENSION(knon,klev) :: ypphi REAL, DIMENSION(knon,klev) :: ydthetadz REAL, DIMENSION(knon) :: ydthetadz300 REAL, DIMENSION(knon) :: Ampl REAL, DIMENSION(knon, nbtersrf) :: yfrac_tersrf, yz0m_tersrf, yz0h_tersrf REAL, DIMENSION(knon) :: yzxtsol ! temperature at surface REAL, DIMENSION(knon) :: ypblh_tmp ! temporaire pblh compressed REAL, DIMENSION(knon) :: zgeo1_over_RG #ifdef ISO REAL, DIMENSION(knon) :: h1 INTEGER :: ixt #endif IF (using_xios) THEN missing_val=missing_val_xios ELSE missing_val=missing_val_netcdf ENDIF yus0(:)=0. ; yvs0(:)=0. ! 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 !!! ! 2b) Initialization of all local variables that will be compressed later !**************************************************************************************** ypct = 0.0 ; yts = 0.0 ; ysnow = 0.0 yqsurf = 0.0 ; yalb = 0.0 ; yalb_vis = 0.0 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 yfluxlat=0.0 ; y_flux0(:)=0.0 yqsol = 0.0 ; yzxtsol = 0.0 ytke=0. yeps=0. yri0(:)=0. y_treedrg=0. 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 y_d_qbs(:,:)=0.0 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. yfluxlat_x=0. ; yfluxlat_w=0. ywake_s=0. ; ywake_cstar=0. ;ywake_dens=0. 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. y_d_u_frein(:,:)=0. y_d_v_frein(:,:)=0. #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 y_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. ypphi = 0.0 !**************************************************************************************** ! 5) Compress variables ! !**************************************************************************************** ! 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) ! 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) 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 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) ywindsp(j) = windsp(i,nsrf) yzmea(j) = zmea(i) yzsig(j) = zsig(i) ycldt(j) = cldt(i) yrmu0(j) = rmu0(i) 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 DO ixt=1,niso yRland_ice(ixt,j)= Rland_ice(ixt,i) ENDDO #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 !--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 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 DO k = 1, klev+1 DO j = 1, knon i = ni(j) ytke(j,k) = tke_x(i,k,nsrf) ENDDO ENDDO 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) 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 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 DO k = 1, klev+1 DO j = 1, knon i = ni(j) 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) ENDDO ENDDO 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 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. ! !**************************************************************************************** IF (iflag_split .eq.0) 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)) 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 ! zgeo1_over_RG(1:knon) = zgeo1(1:knon)/RG CALL eff_surf_param(knon, nbtersrf, yz0m_tersrf, yfrac_tersrf, 'CDN', yz0m, zgeo1_over_RG) CALL eff_surf_param(knon, nbtersrf, yz0h_tersrf, yfrac_tersrf, 'CDN', yz0h, zgeo1_over_RG) ! ENDIF ENDIF ! ypblh_tmp(:)=s_pblh(ni(:)) CALL cdrag(knon, nsrf, & speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1), ypblh_tmp, & yts, yqsurf, yz0m, yz0h, yri0, 0, & ycdragm, ycdragh, zri1, pref, rain_f, yzxtsol, ypplay(:,1)) s_pblh(ni(:)) = ypblh_tmp(:) ! --- 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.+vent)/vent ENDDO ENDIF IF (prt_level >=10) print *,'cdrag -> ycdragh ', ycdragh(1:knon) ELSE !(iflag_split .eq.0) 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 ypblh_tmp(:)=s_pblh_x(ni(:)) CALL cdrag(knon, nsrf, & speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),ypblh,& yts_x, yqsurf_x, yz0m, yz0h, yri0, 0, & ycdragm_x, ycdragh_x, zri1_x, pref_x, rain_f, yzxtsol, ypplay(:,1) ) s_pblh_x(ni(:)) = ypblh_tmp(:) ! --- special Dice. JYG+MPL 25112013 IF (ok_prescr_ust) THEN DO i = 1, knon 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 ENDDO ENDIF IF (prt_level >=10) print *,'clcdrag -> ycdragh_x ', ycdragh_x(1:knon) 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 ypblh_tmp(:)=s_pblh_w(ni(:)) 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, yzxtsol, ypplay(:,1) ) s_pblh_w(ni(:)) = ypblh_tmp(:) ! 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 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 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. ! !**************************************************************************************** IF (iflag_split .eq.0) THEN 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, ycdragh,yus0,yvs0,yts, & yu, yv,yt,yq,ypplay,ypaprs, & ytke,yeps, ycoefm, ycoefh) ELSE CALL coef_diff_turb(dtime, nsrf, knon, ni, & ypaprs, ypplay, yu, yv, yq, yt, yts, yqsurf, ycdragm, & ycoefm, ycoefh, ytke, yeps, y_treedrg) 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) !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,ycdragh_w,yus0,yvs0,yts_w, & yu_w, yv_w, yt_w, yq_w, ypplay, ypaprs, & ytke_w, yeps_w, ycoefm_w, ycoefh_w) 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) 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 IF (iflag_split .eq.0) THEN CALL climb_hq_down(knon, ni, ycoefh, ypaprs, ypplay, & ydelp, yt, yq, dtime, & 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, & 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, & 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 IF (iflag_split .eq.0) THEN CALL climb_wind_down(knon, ni, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv, & 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, & 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, & 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 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 ENDDO ! j = 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(knon, 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) ELSE CALL stdlevvar(knon, 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, yzxtsol) 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) CALL surf_land(itap, dtime, date0, jour, knon, ni,& rlon, rlat, yrmu0, & debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, & 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 & #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 IF (ifl_pbltree .ge. 1) THEN CALL freinage(knon, knon, yu, yv, yt, & 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 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 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) 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, yicemelt, 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 & ) 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_ice(i) = icesub_ice(i) + yicesub(j)*ypct(j) icemelt_ice(i) = icemelt_ice(i) + yicemelt(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 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) ! 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 !ym minloc does't work on GPU (nvfortran + openacc) ! DO j = 1, knon ! k= minloc(abs(ypphii(j,:)-300),1) ! ydthetadz300(j)=ydthetadz(j,k) ! ENDDO ydthetadz300(1:knon) = ydthetadz(1:knon,1) DO k=2, klev DO j=1,knon IF (abs(ypphii(j,k)-300) < ydthetadz300(j)) ydthetadz300(j) = abs(ypphii(j,k)-300) ENDDO ENDDO end if zgeo1_over_RG(1:knon) = zgeo1(1:knon)/RG CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb_vis, & ywindsp, yrmu0, yfder, yts, & itap, dtime, jour, knon, ni, & ypplay(:,1), zgeo1_over_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), & ydt_ds(:knon), ytkt(:knon), ytks(:knon), ytaur(:knon), ysss, & ydthetadz300,Ampl & #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( & rlon, rlat, ysolsw, ysollw, yalb_vis, yfder, & itap, dtime, jour, knon, ni, & lafin, & 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, & yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yicesub, yicemelt, yfluxsens,yfluxlat,& ytsurf_new, y_dflux_t, y_dflux_q, & 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 & #ifdef ISO & ,yxtrain_f, yxtsnow_f,yxt1,Roce, & & yxtsnow,yxtsol,yxtevap,Rland_ice & #endif & ) DO j = 1, knon i = ni(j) icesub_ice(i) = icesub_ice(i) + yicesub(j)*ypct(j) icemelt_ice(i) = icemelt_ice(i) + yicemelt(j)*ypct(j) ENDDO ! 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 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' !ym 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. ! !**************************************************************************************** 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, & 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) 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) ENDIF ! (iflag_split .ge.1) IF (iflag_split .eq.0) THEN CALL climb_hq_up(knon, ni, dtime, yt, yq, & y_flux_q1, y_flux_t1, ypaprs, ypplay, & 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, & 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, & 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) IF (iflag_split .eq.0) THEN CALL climb_wind_up(knon, ni, dtime, yu, yv, y_flux_u1, y_flux_v1, & 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 ELSE !(iflag_split .eq.0) CALL climb_wind_up( knon, ni, dtime, yu_x, yv_x, y_flux_u1_x, y_flux_v1_x, & 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 CALL climb_wind_up(knon, ni, dtime, yu_w, yv_w, y_flux_u1_w, y_flux_v1_w, & 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 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 !**************************************************************************************** ! 13) Transform variables for output format : ! - Decompress ! - Multiply with pourcentage of current surface ! - Cumulate in global variable ! !**************************************************************************************** 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) IF (nsrf .EQ. is_ter .and. ifl_pbltree .GE. 1) THEN 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 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 IF (iflag_split .ge.1) THEN DO j = 1, knon i = ni(j) fluxlat_x(i,nsrf) = yfluxlat_x(j) fluxlat_w(i,nsrf) = yfluxlat_w(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) IF (iflag_split .eq.0) THEN wake_dltke(:,:,nsrf) = 0. DO k = 1, klev+1 DO j = 1, knon i = ni(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) 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) 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) 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 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 DO k = 1, nsoilmx DO j = 1, knon i = ni(j) ftsoil(i, k, nsrf) = ytsoil(j,k) ENDDO ENDDO #ifdef ISO #ifdef ISOVERIF 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 IF (iflag_split .ge.1) THEN 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 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 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 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 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_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 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 ! Calculate the temperature and relative humidity at 2m and the wind at 10m IF (iflag_split .eq.0) THEN IF (iflag_new_t2mq2m==1) THEN CALL stdlevvarn(knon, knon, nsrf, zxli, & uzon, vmer, tair1, qair1, zgeo1, & tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, & yt2m, yq2m, yt10m, yq10m, yu10m, yustar, & yn2mout(:, :, :)) ELSE CALL stdlevvar(knon, knon, nsrf, zxli, & uzon, vmer, tair1, qair1, zgeo1, & tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, & yt2m, yq2m, yt10m, yq10m, yu10m, yustar, ypblh, rain_f, yzxtsol) ENDIF ELSE !(iflag_split .eq.0) IF (iflag_new_t2mq2m==1) THEN CALL stdlevvarn(knon, 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(:, :, :)) CALL stdlevvarn(knon, 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(:, :, :)) ELSE CALL stdlevvar(knon, 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, yzxtsol) CALL stdlevvar(knon, 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, yzxtsol) ENDIF ENDIF ! (iflag_split .eq.0) 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) !IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique !IM Ajoute dependance type surface IF (thermcep) THEN 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 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) 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) #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 END SUBROUTINE pbl_surface_subsrf END MODULE pbl_surface_subsrf_mod