! ! $Id: pbl_surface_mod.F90 1932 2014-01-17 18:03:16Z musat $ ! MODULE pbl_surface_mod ! ! Planetary Boundary Layer and Surface module ! ! This module manage the calculation of turbulent diffusion in the boundary layer ! and all interactions towards the differents sub-surfaces. ! ! USE dimphy USE mod_phys_lmdz_para, ONLY : mpi_size USE ioipsl USE surface_data, ONLY : type_ocean, ok_veget 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_wind_mod, ONLY : climb_wind_down, climb_wind_up USE coef_diff_turb_mod, ONLY : coef_diff_turb USE control_mod IMPLICIT NONE ! Declaration of variables saved in restart file REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: qsol ! water height in the soil (mm) !$OMP THREADPRIVATE(qsol) REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: fder ! flux drift !$OMP THREADPRIVATE(fder) REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: snow ! snow at surface !$OMP THREADPRIVATE(snow) REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: qsurf ! humidity at surface !$OMP THREADPRIVATE(qsurf) REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: evap ! evaporation at surface !$OMP THREADPRIVATE(evap) REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: rugos ! rugosity at surface (m) !$OMP THREADPRIVATE(rugos) REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: agesno ! age of snow at surface !$OMP THREADPRIVATE(agesno) ! Correction pour le cas AMMA (PRIVATE) REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: ftsoil ! soil temperature !$OMP THREADPRIVATE(ftsoil) CONTAINS ! !**************************************************************************************** ! SUBROUTINE pbl_surface_init(qsol_rst, fder_rst, snow_rst, qsurf_rst,& evap_rst, rugos_rst, agesno_rst, ftsoil_rst) ! This routine should be called after the restart file has been read. ! This routine initialize the restart variables and does some validation tests ! for the index of the different surfaces and tests the choice of type of ocean. USE indice_sol_mod INCLUDE "dimsoil.h" INCLUDE "iniprint.h" ! Input variables !**************************************************************************************** REAL, DIMENSION(klon), INTENT(IN) :: qsol_rst REAL, DIMENSION(klon), INTENT(IN) :: fder_rst REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: snow_rst REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: qsurf_rst REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: evap_rst REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: rugos_rst REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: agesno_rst REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(IN) :: ftsoil_rst ! Local variables !**************************************************************************************** INTEGER :: ierr CHARACTER(len=80) :: abort_message CHARACTER(len = 20) :: modname = 'pbl_surface_init' !**************************************************************************************** ! Allocate and initialize module variables with fields read from restart file. ! !**************************************************************************************** ALLOCATE(qsol(klon), stat=ierr) IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1) ALLOCATE(fder(klon), stat=ierr) IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1) ALLOCATE(snow(klon,nbsrf), stat=ierr) IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1) ALLOCATE(qsurf(klon,nbsrf), stat=ierr) IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1) ALLOCATE(evap(klon,nbsrf), stat=ierr) IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1) ALLOCATE(rugos(klon,nbsrf), stat=ierr) IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1) ALLOCATE(agesno(klon,nbsrf), stat=ierr) IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1) ALLOCATE(ftsoil(klon,nsoilmx,nbsrf), stat=ierr) IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1) qsol(:) = qsol_rst(:) fder(:) = fder_rst(:) snow(:,:) = snow_rst(:,:) qsurf(:,:) = qsurf_rst(:,:) evap(:,:) = evap_rst(:,:) rugos(:,:) = rugos_rst(:,:) agesno(:,:) = agesno_rst(:,:) ftsoil(:,:,:) = ftsoil_rst(:,:,:) !**************************************************************************************** ! Test for sub-surface indices ! !**************************************************************************************** IF (is_ter /= 1) THEN WRITE(lunout,*)" *** Warning ***" WRITE(lunout,*)" is_ter n'est pas le premier surface, is_ter = ",is_ter WRITE(lunout,*)"or on doit commencer par les surfaces continentales" abort_message="voir ci-dessus" CALL abort_gcm(modname,abort_message,1) ENDIF IF ( is_oce > is_sic ) THEN WRITE(lunout,*)' *** Warning ***' WRITE(lunout,*)' Pour des raisons de sequencement dans le code' WRITE(lunout,*)' l''ocean doit etre traite avant la banquise' WRITE(lunout,*)' or is_oce = ',is_oce, '> is_sic = ',is_sic abort_message='voir ci-dessus' CALL abort_gcm(modname,abort_message,1) ENDIF IF ( is_lic > is_sic ) THEN WRITE(lunout,*)' *** Warning ***' WRITE(lunout,*)' Pour des raisons de sequencement dans le code' WRITE(lunout,*)' la glace contineltalle doit etre traite avant la glace de mer' WRITE(lunout,*)' or is_lic = ',is_lic, '> is_sic = ',is_sic abort_message='voir ci-dessus' CALL abort_gcm(modname,abort_message,1) ENDIF !**************************************************************************************** ! Validation of ocean mode ! !**************************************************************************************** IF (type_ocean /= 'slab ' .AND. type_ocean /= 'force ' .AND. type_ocean /= 'couple') THEN WRITE(lunout,*)' *** Warning ***' WRITE(lunout,*)'Option couplage pour l''ocean = ', type_ocean abort_message='option pour l''ocean non valable' CALL abort_gcm(modname,abort_message,1) ENDIF END SUBROUTINE pbl_surface_init ! !**************************************************************************************** ! SUBROUTINE pbl_surface( & dtime, date0, itap, jour, & debut, lafin, & rlon, rlat, rugoro, rmu0, & zsig, sollwd_m, pphi, cldt, & rain_f, snow_f, solsw_m, sollw_m, & t, q, u, v, & pplay, paprs, pctsrf, & ts, alb1, alb2,ustar, u10m, v10m,wstar, & lwdown_m, cdragh, cdragm, zu1, zv1, & alb1_m, alb2_m, zxsens, zxevap, & alb3_lic, runoff, snowhgt, qsnow, to_ice, sissnow, & zxtsol, zxfluxlat, zt2m, qsat2m, & d_t, d_q, d_u, d_v, d_t_diss, & zcoefh, zcoefm, slab_wfbils, & qsol_d, zq2m, s_pblh, s_plcl, & s_capCL, s_oliqCL, s_cteiCL, s_pblT, & s_therm, s_trmb1, s_trmb2, s_trmb3, & zxrugs,zustar,zu10m, zv10m, fder_print, & zxqsurf, rh2m, zxfluxu, zxfluxv, & rugos_d, agesno_d, sollw, solsw, & d_ts, evap_d, fluxlat, t2m, & wfbils, wfbilo, flux_t, flux_u, flux_v,& dflux_t, dflux_q, zxsnow, & zxfluxt, zxfluxq, q2m, flux_q, tke ) !**************************************************************************************** ! 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 ! 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 ! rugos----input-R- longeur de rugosite (en m) ! Martin ! zsig-----input-R- slope ! cldt-----input-R- total cloud fraction ! pphi-----input-R- geopotentiel de chaque couche (g z) (reference sol) ! Martin ! ! 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---input/output-R- tke (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 ! USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, co2_send USE indice_sol_mod IMPLICIT NONE INCLUDE "dimsoil.h" INCLUDE "YOMCST.h" INCLUDE "iniprint.h" INCLUDE "YOETHF.h" INCLUDE "FCTTRE.h" INCLUDE "clesphys.h" INCLUDE "compbl.h" INCLUDE "dimensions.h" INCLUDE "temps.h" INCLUDE "flux_arp.h" !**************************************************************************************** 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) :: solsw_m ! net shortwave radiation at mean surface REAL, DIMENSION(klon), INTENT(IN) :: sollw_m ! net longwave radiation at mean surface REAL, DIMENSION(klon,klev), INTENT(IN) :: t ! temperature (K) REAL, DIMENSION(klon,klev), INTENT(IN) :: q ! water vapour (kg/kg) REAL, DIMENSION(klon,klev), INTENT(IN) :: u ! u speed REAL, DIMENSION(klon,klev), INTENT(IN) :: v ! v speed REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! mid-layer pression (Pa) REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! pression between layers (Pa) REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf ! sub-surface fraction ! Martin REAL, DIMENSION(klon), INTENT(IN) :: zsig ! slope REAL, DIMENSION(klon), INTENT(IN) :: sollwd_m ! net longwave radiation at mean s REAL, DIMENSION(klon), INTENT(IN) :: cldt ! total cloud fraction REAL, DIMENSION(klon,klev), INTENT(IN) :: pphi ! geopotential (m2/s2) ! Martin ! Input/Output variables !**************************************************************************************** REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: ts ! temperature at surface (K) REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: alb1 ! albedo in visible SW interval REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: alb2 ! albedo in near infra-red SW interval 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 ! Output variables !**************************************************************************************** REAL, DIMENSION(klon), INTENT(OUT) :: lwdown_m ! Downcoming longwave radiation REAL, DIMENSION(klon), INTENT(OUT) :: cdragh ! drag coefficient for T and Q REAL, DIMENSION(klon), INTENT(OUT) :: cdragm ! drag coefficient for wind REAL, DIMENSION(klon), INTENT(OUT) :: zu1 ! u wind speed in first layer REAL, DIMENSION(klon), INTENT(OUT) :: zv1 ! v wind speed in first layer REAL, DIMENSION(klon), INTENT(OUT) :: alb1_m ! mean albedo in visible SW interval REAL, DIMENSION(klon), INTENT(OUT) :: alb2_m ! mean albedo in near IR SW interval ! Martin REAL, DIMENSION(klon), INTENT(OUT) :: alb3_lic ! Martin REAL, DIMENSION(klon), INTENT(OUT) :: zxsens ! sensible heat flux at surface with inversed sign ! (=> positive sign upwards) REAL, DIMENSION(klon), INTENT(OUT) :: zxevap ! water vapour flux at surface, positiv upwards REAL, DIMENSION(klon), INTENT(OUT) :: zxtsol ! temperature at surface, mean for each grid point REAL, DIMENSION(klon), INTENT(OUT) :: zxfluxlat ! latent flux, mean for each grid point REAL, DIMENSION(klon), INTENT(OUT) :: zt2m ! temperature at 2m, mean for each grid point REAL, DIMENSION(klon), INTENT(OUT) :: qsat2m REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_t ! change in temperature REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_t_diss ! change in temperature REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_q ! change in water vapour REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_u ! change in u speed REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_v ! change in v speed REAL, INTENT(OUT):: zcoefh(:, :, :) ! (klon, klev, nbsrf + 1) ! coef for turbulent diffusion of T and Q, mean for each grid point REAL, INTENT(OUT):: zcoefm(:, :, :) ! (klon, klev, nbsrf + 1) ! coef for turbulent diffusion of U and V (?), mean for each grid point ! Output only for diagnostics REAL, DIMENSION(klon), INTENT(OUT) :: slab_wfbils! heat balance at surface only for slab at ocean points REAL, DIMENSION(klon), INTENT(OUT) :: qsol_d ! water height in the soil (mm) REAL, DIMENSION(klon), INTENT(OUT) :: zq2m ! water vapour at 2m, mean for each grid point REAL, DIMENSION(klon), INTENT(OUT) :: s_pblh ! height of the planetary boundary layer(HPBL) REAL, DIMENSION(klon), INTENT(OUT) :: s_plcl ! condensation level REAL, DIMENSION(klon), INTENT(OUT) :: s_capCL ! CAPE of PBL REAL, DIMENSION(klon), INTENT(OUT) :: s_oliqCL ! liquid water intergral of PBL REAL, DIMENSION(klon), INTENT(OUT) :: s_cteiCL ! cloud top instab. crit. of PBL REAL, DIMENSION(klon), INTENT(OUT) :: s_pblT ! temperature at PBLH REAL, DIMENSION(klon), INTENT(OUT) :: s_therm ! thermal virtual temperature excess REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb1 ! deep cape, mean for each grid point REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb2 ! inhibition, mean for each grid point REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb3 ! point Omega, mean for each grid point REAL, DIMENSION(klon), INTENT(OUT) :: zxrugs ! rugosity at surface (m), mean for each grid point REAL, DIMENSION(klon), INTENT(OUT) :: zustar ! u* REAL, DIMENSION(klon), INTENT(OUT) :: zu10m ! u speed at 10m, mean for each grid point REAL, DIMENSION(klon), INTENT(OUT) :: zv10m ! v speed at 10m, mean for each grid point REAL, DIMENSION(klon), INTENT(OUT) :: fder_print ! fder for printing (=fder(i) + dflux_t(i) + dflux_q(i)) REAL, DIMENSION(klon), INTENT(OUT) :: zxqsurf ! humidity at surface, mean for each grid point REAL, DIMENSION(klon), INTENT(OUT) :: rh2m ! relative humidity at 2m REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxu ! u wind tension, mean for each grid point REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxv ! v wind tension, mean for each grid point REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: rugos_d ! rugosity length (m) REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: agesno_d ! age of snow at surface REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: solsw ! net shortwave radiation at surface REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: sollw ! net longwave radiation at surface REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: d_ts ! change in temperature at surface REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: evap_d ! evaporation at surface REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: fluxlat ! latent flux REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: t2m ! temperature at 2 meter height REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfbils ! heat balance at surface REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfbilo ! water balance at surface REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_t ! sensible heat flux (CpT) J/m**2/s (W/m**2) ! positve orientation downwards REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_u ! u wind tension (kg m/s)/(m**2 s) or Pascal REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_v ! v wind tension (kg m/s)/(m**2 s) or Pascal ! Output not needed REAL, DIMENSION(klon), INTENT(OUT) :: dflux_t ! change of sensible heat flux REAL, DIMENSION(klon), INTENT(OUT) :: dflux_q ! change of water vapour flux REAL, DIMENSION(klon), INTENT(OUT) :: zxsnow ! snow at surface, mean for each grid point REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxt ! sensible heat flux, mean for each grid point REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxq ! water vapour flux, mean for each grid point REAL, DIMENSION(klon, nbsrf),INTENT(OUT) :: q2m ! water vapour at 2 meter height REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_q ! water vapour flux(latent flux) (kg/m**2/s) ! Martin ! sisvat REAL, DIMENSION(klon), INTENT(OUT) :: qsnow ! snow water content REAL, DIMENSION(klon), INTENT(OUT) :: snowhgt ! snow height REAL, DIMENSION(klon), INTENT(OUT) :: to_ice ! snow passed to ice REAL, DIMENSION(klon), INTENT(OUT) :: sissnow ! snow in snow model REAL, DIMENSION(klon), INTENT(OUT) :: runoff ! runoff on land ice ! Martin ! Local variables with attribute SAVE !**************************************************************************************** INTEGER, SAVE :: nhoridbg, nidbg ! variables for IOIPSL !$OMP THREADPRIVATE(nhoridbg, nidbg) LOGICAL, SAVE :: debugindex=.FALSE. !$OMP THREADPRIVATE(debugindex) LOGICAL, SAVE :: first_call=.TRUE. !$OMP THREADPRIVATE(first_call) CHARACTER(len=8), DIMENSION(nbsrf), SAVE :: cl_surf !$OMP THREADPRIVATE(cl_surf) ! Other local variables !**************************************************************************************** INTEGER :: i, k, nsrf INTEGER :: knon, j INTEGER :: idayref INTEGER , DIMENSION(klon) :: ni REAL :: zx_alf1, zx_alf2 !valeur ambiante par extrapola REAL :: amn, amx REAL :: f1 ! fraction de longeurs visibles parmi tout SW intervalle REAL, DIMENSION(klon) :: r_co2_ppm ! taux CO2 atmosphere REAL, DIMENSION(klon) :: yts, yrugos, ypct, yz0_new REAL, DIMENSION(klon) :: yalb, yalb1, yalb2 REAL, DIMENSION(klon) :: yu1, yv1,ytoto REAL, DIMENSION(klon) :: ysnow, yqsurf, yagesno, yqsol REAL, DIMENSION(klon) :: yrain_f, ysnow_f REAL, DIMENSION(klon) :: ysolsw, ysollw REAL, DIMENSION(klon) :: yfder REAL, DIMENSION(klon) :: yrugoro REAL, DIMENSION(klon) :: yfluxlat REAL, DIMENSION(klon) :: y_d_ts REAL, DIMENSION(klon) :: y_flux_t1, y_flux_q1 REAL, DIMENSION(klon) :: y_dflux_t, y_dflux_q REAL, DIMENSION(klon) :: y_flux_u1, y_flux_v1 REAL, DIMENSION(klon) :: yt2m, yq2m, yu10m REAL, DIMENSION(klon) :: yustar REAL, DIMENSION(klon) :: ywstar REAL, DIMENSION(klon) :: ywindsp REAL, DIMENSION(klon) :: yt10m, yq10m REAL, DIMENSION(klon) :: ypblh REAL, DIMENSION(klon) :: ylcl REAL, DIMENSION(klon) :: ycapCL REAL, DIMENSION(klon) :: yoliqCL REAL, DIMENSION(klon) :: ycteiCL REAL, DIMENSION(klon) :: ypblT REAL, DIMENSION(klon) :: ytherm REAL, DIMENSION(klon) :: ytrmb1 REAL, DIMENSION(klon) :: ytrmb2 REAL, DIMENSION(klon) :: ytrmb3 REAL, DIMENSION(klon) :: uzon, vmer REAL, DIMENSION(klon) :: tair1, qair1, tairsol REAL, DIMENSION(klon) :: psfce, patm REAL, DIMENSION(klon) :: qairsol, zgeo1 REAL, DIMENSION(klon) :: rugo1 REAL, DIMENSION(klon) :: yfluxsens REAL, DIMENSION(klon) :: AcoefH, AcoefQ, BcoefH, BcoefQ REAL, DIMENSION(klon) :: AcoefU, AcoefV, BcoefU, BcoefV REAL, DIMENSION(klon) :: ypsref REAL, DIMENSION(klon) :: yevap, ytsurf_new, yalb1_new, yalb2_new, yalb3_new REAL, DIMENSION(klon) :: ztsol REAL, DIMENSION(klon) :: alb_m ! mean albedo for whole SW interval REAL, DIMENSION(klon,klev) :: y_d_t, y_d_q, y_d_t_diss REAL, DIMENSION(klon,klev) :: y_d_u, y_d_v REAL, DIMENSION(klon,klev) :: y_flux_t, y_flux_q REAL, DIMENSION(klon,klev) :: y_flux_u, y_flux_v REAL, DIMENSION(klon,klev) :: ycoefh, ycoefm,ycoefq REAL, DIMENSION(klon) :: ycdragh, ycdragm REAL, DIMENSION(klon,klev) :: yu, yv REAL, DIMENSION(klon,klev) :: yt, yq REAL, DIMENSION(klon,klev) :: ypplay, ydelp REAL, DIMENSION(klon,klev) :: delp REAL, DIMENSION(klon,klev+1) :: ypaprs REAL, DIMENSION(klon,klev+1) :: ytke REAL, DIMENSION(klon,nsoilmx) :: ytsoil 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(klon) :: Kech_h ! Coefficient d'echange pour l'energie ! For debugging with IOIPSL INTEGER, DIMENSION(iim*(jjm+1)) :: ndexbg REAL :: zjulian REAL, DIMENSION(klon) :: tabindx REAL, DIMENSION(iim,jjm+1) :: zx_lon, zx_lat REAL, DIMENSION(iim,jjm+1) :: debugtab REAL, DIMENSION(klon,nbsrf) :: pblh ! height of the planetary boundary layer REAL, DIMENSION(klon,nbsrf) :: plcl ! condensation level REAL, DIMENSION(klon,nbsrf) :: capCL REAL, DIMENSION(klon,nbsrf) :: oliqCL REAL, DIMENSION(klon,nbsrf) :: cteiCL REAL, DIMENSION(klon,nbsrf) :: pblT REAL, DIMENSION(klon,nbsrf) :: therm REAL, DIMENSION(klon,nbsrf) :: trmb1 ! deep cape REAL, DIMENSION(klon,nbsrf) :: trmb2 ! inhibition REAL, DIMENSION(klon,nbsrf) :: trmb3 ! point Omega REAL, DIMENSION(klon,nbsrf) :: zx_rh2m, zx_qsat2m REAL, DIMENSION(klon,nbsrf) :: zx_t1 REAL, DIMENSION(klon, nbsrf) :: alb ! mean albedo for whole SW interval REAL, DIMENSION(klon) :: ylwdown ! jg : temporary (ysollwdown) REAL :: zx_qs1, zcor1, zdelta1 ! Martin REAL, DIMENSION(klon, nbsrf) :: sollwd ! net longwave radiation at surface REAL, DIMENSION(klon) :: ysollwd REAL, DIMENSION(klon) :: ytoice REAL, DIMENSION(klon) :: ysnowhgt, yqsnow, ysissnow, yrunoff REAL, DIMENSION(klon) :: yzsig REAL, DIMENSION(klon,klev) :: ypphi REAL, DIMENSION(klon) :: ycldt REAL, DIMENSION(klon) :: yrmu0 ! Martin !**************************************************************************************** ! End of declarations !**************************************************************************************** !**************************************************************************************** ! 1) Initialisation and validation tests ! Only done first time entering this subroutine ! !**************************************************************************************** IF (first_call) THEN first_call=.FALSE. ! Initialize ok_flux_surf (for 1D model) if (klon>1) ok_flux_surf=.FALSE. ! Initilize debug IO IF (debugindex .AND. mpi_size==1) THEN ! initialize IOIPSL output idayref = day_ini CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian) CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon) DO i = 1, iim zx_lon(i,1) = rlon(i+1) zx_lon(i,jjm+1) = rlon(i+1) ENDDO CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat) CALL histbeg("sous_index", iim,zx_lon(:,1),jjm+1,zx_lat(1,:), & 1,iim,1,jjm+1, & itau_phy,zjulian,dtime,nhoridbg,nidbg) ! no vertical axis cl_surf(1)='ter' cl_surf(2)='lic' cl_surf(3)='oce' cl_surf(4)='sic' DO nsrf=1,nbsrf CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim, & jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime) END DO CALL histend(nidbg) CALL histsync(nidbg) END IF ENDIF !**************************************************************************************** ! Force soil water content to qsol0 if qsol0>0 and VEGET=F (use bucket ! instead of ORCHIDEE) IF (qsol0>=0.) THEN PRINT*,'WARNING : On impose qsol=',qsol0 qsol(:)=qsol0 ENDIF !**************************************************************************************** !**************************************************************************************** ! 2) Initialization to zero ! Done for all local variables that will be compressed later ! and argument with INTENT(OUT) !**************************************************************************************** cdragh = 0.0 ; cdragm = 0.0 ; dflux_t = 0.0 ; dflux_q = 0.0 ypct = 0.0 ; yts = 0.0 ; ysnow = 0.0 zv1 = 0.0 ; yqsurf = 0.0 ; yalb1 = 0.0 ; yalb2 = 0.0 yrain_f = 0.0 ; ysnow_f = 0.0 ; yfder = 0.0 ; ysolsw = 0.0 ysollw = 0.0 ; yrugos = 0.0 ; yu1 = 0.0 yv1 = 0.0 ; ypaprs = 0.0 ; ypplay = 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 yrugoro = 0.0 ; ywindsp = 0.0 d_ts = 0.0 ; yfluxlat=0.0 ; flux_t = 0.0 ; flux_q = 0.0 flux_u = 0.0 ; flux_v = 0.0 ; d_t = 0.0 ; d_q = 0.0 d_t_diss= 0.0 ;d_u = 0.0 ; d_v = 0.0 ; yqsol = 0.0 ytherm = 0.0 ; ytke=0. ! Martin ysnowhgt = 0.0; yqsnow = 0.0 ; yrunoff = 0.0 ; ytoice =0.0 yalb3_new = 0.0 ; ysissnow = 0.0 ; ysollwd = 0.0 ypphi = 0.0 ; ycldt = 0.0 ; yrmu0 = 0.0 ! Martin tke(:,:,is_ave)=0. IF (iflag_pbl<20.or.iflag_pbl>=30) THEN zcoefh(:,:,:) = 0.0 zcoefh(:,1,:) = 999999. ! zcoefh(:,k=1) should never be used zcoefm(:,:,:) = 0.0 zcoefm(:,1,:) = 999999. ! ELSE zcoefm(:,:,is_ave)=0. zcoefh(:,:,is_ave)=0. ENDIF ytsoil = 999999. rh2m(:) = 0. qsat2m(:) = 0. !**************************************************************************************** ! 3) - Calculate pressure thickness of each layer ! - Calculate the wind at first layer ! - Mean calculations of albedo ! - Calculate net radiance at sub-surface !**************************************************************************************** DO k = 1, klev DO i = 1, klon delp(i,k) = paprs(i,k)-paprs(i,k+1) ENDDO ENDDO !**************************************************************************************** ! Test for rugos........ from physiq.. A la fin plutot??? ! !**************************************************************************************** zxrugs(:) = 0.0 DO nsrf = 1, nbsrf DO i = 1, klon rugos(i,nsrf) = MAX(rugos(i,nsrf),0.000015) zxrugs(i) = zxrugs(i) + rugos(i,nsrf)*pctsrf(i,nsrf) ENDDO ENDDO ! Mean calculations of albedo ! ! Albedo at sub-surface ! * alb1 : albedo in visible SW interval ! * alb2 : albedo in near infrared SW interval ! * alb : mean albedo for whole SW interval ! ! Mean albedo for grid point ! * alb1_m : albedo in visible SW interval ! * alb2_m : albedo in near infrared SW interval ! * alb_m : mean albedo at whole SW interval alb1_m(:) = 0.0 alb2_m(:) = 0.0 DO nsrf = 1, nbsrf DO i = 1, klon alb1_m(i) = alb1_m(i) + alb1(i,nsrf) * pctsrf(i,nsrf) alb2_m(i) = alb2_m(i) + alb2(i,nsrf) * pctsrf(i,nsrf) ENDDO ENDDO ! We here suppose the fraction f1 of incoming radiance of visible radiance ! as a fraction of all shortwave radiance f1 = 0.5 ! f1 = 1 ! put f1=1 to recreate old calculations DO nsrf = 1, nbsrf DO i = 1, klon alb(i,nsrf) = f1*alb1(i,nsrf) + (1-f1)*alb2(i,nsrf) ENDDO ENDDO DO i = 1, klon alb_m(i) = f1*alb1_m(i) + (1-f1)*alb2_m(i) END DO ! Calculation of mean temperature at surface grid points ztsol(:) = 0.0 DO nsrf = 1, nbsrf DO i = 1, klon ztsol(i) = ztsol(i) + ts(i,nsrf)*pctsrf(i,nsrf) ENDDO ENDDO ! Linear distrubution on sub-surface of long- and shortwave net radiance DO nsrf = 1, nbsrf DO i = 1, klon sollw(i,nsrf) = sollw_m(i) + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ts(i,nsrf)) ! Martin sollwd(i,nsrf)= sollwd_m(i) ! Martin solsw(i,nsrf) = solsw_m(i) * (1.-alb(i,nsrf)) / (1.-alb_m(i)) ENDDO ENDDO ! Downwelling longwave radiation at mean surface lwdown_m(:) = 0.0 DO i = 1, klon lwdown_m(i) = sollw_m(i) + RSIGMA*ztsol(i)**4 ENDDO !**************************************************************************************** ! 4) Loop over different surfaces ! ! Only points containing a fraction of the sub surface will be threated. ! !**************************************************************************************** loop_nbsrf: DO nsrf = 1, nbsrf ! 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 ! write index, with IOIPSL IF (debugindex .AND. mpi_size==1) THEN tabindx(:)=0. DO i=1,knon tabindx(i)=REAL(i) END DO debugtab(:,:) = 0. ndexbg(:) = 0 CALL gath2cpl(tabindx,debugtab,knon,ni) CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,iim*(jjm+1), ndexbg) ENDIF !**************************************************************************************** ! 5) Compress variables ! !**************************************************************************************** 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) yalb1(j) = alb1(i,nsrf) yalb2(j) = alb2(i,nsrf) yrain_f(j) = rain_f(i) ysnow_f(j) = snow_f(i) yagesno(j) = agesno(i,nsrf) yfder(j) = fder(i) ysolsw(j) = solsw(i,nsrf) ysollw(j) = sollw(i,nsrf) yrugos(j) = rugos(i,nsrf) yrugoro(j) = rugoro(i) yu1(j) = u(i,1) yv1(j) = v(i,1) ypaprs(j,klev+1) = paprs(i,klev+1) ywindsp(j) = SQRT(u10m(i,nsrf)**2 + v10m(i,nsrf)**2 ) ! Martin yzsig(j) = zsig(i) ycldt(j) = cldt(i) yrmu0(j) = rmu0(i) ! Martin END DO 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) ytke(j,k) = tke(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) ENDDO ENDDO DO k = 1, nsoilmx DO j = 1, knon i = ni(j) ytsoil(j,k) = ftsoil(i,k,nsrf) END DO END DO ! 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) END DO ENDIF !**************************************************************************************** ! 6a) Calculate coefficients for turbulent diffusion at surface, cdragh et cdragm. ! !**************************************************************************************** CALL clcdrag( knon, nsrf, ypaprs, ypplay, & yu(:,1), yv(:,1), yt(:,1), yq(:,1), & yts, yqsurf, yrugos, & ycdragm, ycdragh ) !**************************************************************************************** ! 6b) Calculate coefficients for turbulent diffusion in the atmosphere, ycoefm et ycoefm. ! !**************************************************************************************** CALL coef_diff_turb(dtime, nsrf, knon, ni, & ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf, ycdragm, & ycoefm, ycoefh, ytke) 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 !**************************************************************************************** ! ! 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 CALL climb_hq_down(knon, ycoefh, ypaprs, ypplay, & ydelp, yt, yq, dtime, & AcoefH, AcoefQ, BcoefH, BcoefQ) ! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V CALL climb_wind_down(knon, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv, & AcoefU, AcoefV, BcoefU, BcoefV) !**************************************************************************************** ! 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)) END DO ELSE r_co2_ppm(:) = co2_ppm ! Constant field END IF !**************************************************************************************** ! ! 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)) END DO ! Calculate the temperature et relative humidity at 2m and the wind at 10m CALL stdlevvar(klon, knon, is_ter, zxli, & yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, & yts, yqsurf, yrugos, ypaprs(:,1), ypplay(:,1), & yt2m, yq2m, yt10m, yq10m, yu10m, yustar) END IF !**************************************************************************************** ! ! 10) Switch selon current surface ! It is necessary to start with the continental surfaces because the ocean ! needs their run-off. ! !**************************************************************************************** SELECT CASE(nsrf) CASE(is_ter) ! ylwdown : to be removed, calculation is now done at land surface in surf_land ylwdown(:)=0.0 DO i=1,knon ylwdown(i)=lwdown_m(ni(i)) END DO CALL surf_land(itap, dtime, date0, jour, knon, ni,& rlon, rlat, & debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, & yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),& AcoefH, AcoefQ, BcoefH, BcoefQ, & AcoefU, AcoefV, BcoefU, BcoefV, & ypsref, yu1, yv1, yrugoro, pctsrf, & ylwdown, yq2m, yt2m, & ysnow, yqsol, yagesno, ytsoil, & yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, & yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, & y_flux_u1, y_flux_v1 ) CASE(is_lic) ! Martin CALL surf_landice(itap, dtime, knon, ni, & rlon, rlat, debut, lafin, & yrmu0, ysollwd, yalb, ypphi(:,1), & ysolsw, ysollw, yts, ypplay(:,1), & ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),& AcoefH, AcoefQ, BcoefH, BcoefQ, & AcoefU, AcoefV, BcoefU, BcoefV, & ypsref, yu1, yv1, yrugoro, pctsrf, & ysnow, yqsurf, yqsol, yagesno, & ytsoil, yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, & ytsurf_new, y_dflux_t, y_dflux_q, & yzsig, ycldt, & ysnowhgt, yqsnow, ytoice, ysissnow, & yalb3_new, yrunoff, & y_flux_u1, y_flux_v1) !CALL surf_landice(itap, dtime, knon, ni, & ! ysolsw, ysollw, yts, ypplay(:,1), & ! ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),& ! AcoefH, AcoefQ, BcoefH, BcoefQ, & ! AcoefU, AcoefV, BcoefU, BcoefV, & ! ypsref, yu1, yv1, yrugoro, pctsrf, & ! ysnow, yqsurf, yqsol, yagesno, & ! ytsoil, yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, & ! ytsurf_new, y_dflux_t, y_dflux_q, & ! y_flux_u1, y_flux_v1) alb3_lic(:)=0. 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) END DO ! Martin CASE(is_oce) CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb1, & yrugos, ywindsp, rmu0, yfder, yts, & itap, dtime, jour, knon, ni, & ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),& AcoefH, AcoefQ, BcoefH, BcoefQ, & AcoefU, AcoefV, BcoefU, BcoefV, & ypsref, yu1, yv1, yrugoro, pctsrf, & ysnow, yqsurf, yagesno, & yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, & ytsurf_new, y_dflux_t, y_dflux_q, slab_wfbils, & y_flux_u1, y_flux_v1) CASE(is_sic) CALL surf_seaice( & rlon, rlat, ysolsw, ysollw, yalb1, yfder, & itap, dtime, jour, knon, ni, & lafin, & yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),& AcoefH, AcoefQ, BcoefH, BcoefQ, & AcoefU, AcoefV, BcoefU, BcoefV, & ypsref, yu1, yv1, yrugoro, pctsrf, & ysnow, yqsurf, yqsol, yagesno, ytsoil, & yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, & ytsurf_new, y_dflux_t, y_dflux_q, & y_flux_u1, y_flux_v1) CASE DEFAULT WRITE(lunout,*) 'Surface index = ', nsrf abort_message = 'Surface index not valid' CALL abort_gcm(modname,abort_message,1) END SELECT !**************************************************************************************** ! 11) - Calcul the increment of surface temperature ! !**************************************************************************************** if (evap0>=0.) then yevap(:)=evap0 yevap(:)=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. ! !**************************************************************************************** ! H and Q IF (ok_flux_surf) THEN PRINT *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT y_flux_t1(:) = fsens y_flux_q1(:) = flat/RLVTT yfluxlat(:) = flat Kech_h(:) = ycdragh(:) * (1.0+SQRT(yu(:,1)**2+yv(:,1)**2)) * & ypplay(:,1)/(RD*yt(:,1)) ytoto(:)=(1./RCPD)*(AcoefH(:)+BcoefH(:)*y_flux_t1(:)*dtime) ytsurf_new(:)=ytoto(:)-y_flux_t1(:)/(Kech_h(:)*RCPD) y_d_ts(:) = ytsurf_new(:) - yts(:) ELSE y_flux_t1(:) = yfluxsens(:) y_flux_q1(:) = -yevap(:) ENDIF CALL climb_hq_up(knon, dtime, yt, yq, & y_flux_q1, y_flux_t1, ypaprs, ypplay, & y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:)) CALL climb_wind_up(knon, dtime, yu, yv, y_flux_u1, y_flux_v1, & 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,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,nsrf) ENDIF ! print*,'yamada_c OK' DO j = 1, knon y_dflux_t(j) = y_dflux_t(j) * ypct(j) y_dflux_q(j) = y_dflux_q(j) * ypct(j) ENDDO !**************************************************************************************** ! 13) Transform variables for output format : ! - Decompress ! - Multiply with pourcentage of current surface ! - Cumulate in global variable ! !**************************************************************************************** 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) 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) ENDDO ENDDO ! print*,'Dans pbl OK1' evap(:,nsrf) = - flux_q(:,1,nsrf) alb1(:, nsrf) = 0. alb2(:, nsrf) = 0. snow(:, nsrf) = 0. qsurf(:, nsrf) = 0. rugos(:, nsrf) = 0. fluxlat(:,nsrf) = 0. DO j = 1, knon i = ni(j) d_ts(i,nsrf) = y_d_ts(j) alb1(i,nsrf) = yalb1_new(j) alb2(i,nsrf) = yalb2_new(j) snow(i,nsrf) = ysnow(j) qsurf(i,nsrf) = yqsurf(j) rugos(i,nsrf) = yz0_new(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) dflux_q(i) = dflux_q(i) + y_dflux_q(j) END DO ! print*,'Dans pbl OK2' DO k = 2, klev DO j = 1, knon i = ni(j) tke(i,k,nsrf) = ytke(j,k) zcoefh(i,k,nsrf) = ycoefh(j,k) zcoefm(i,k,nsrf) = ycoefm(j,k) tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)*ypct(j) 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) END DO END DO ! print*,'Dans pbl OK3' IF ( nsrf .EQ. is_ter ) THEN DO j = 1, knon i = ni(j) qsol(i) = yqsol(j) END DO END IF ftsoil(:,:,nsrf) = 0. DO k = 1, nsoilmx DO j = 1, knon i = ni(j) ftsoil(i, k, nsrf) = ytsoil(j,k) END DO END DO 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) 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) END DO END DO ! print*,'Dans pbl OK4' !**************************************************************************************** ! 14) Calculate the temperature et relative humidity at 2m and the wind at 10m ! Call HBTM ! !**************************************************************************************** t2m(:,nsrf) = 0. q2m(:,nsrf) = 0. ustar(:,nsrf) = 0. wstar(:,nsrf) = 0. u10m(:,nsrf) = 0. v10m(:,nsrf) = 0. pblh(:,nsrf) = 0. ! Hauteur de couche limite plcl(:,nsrf) = 0. ! Niveau de condensation de la CLA capCL(:,nsrf) = 0. ! CAPE de couche limite oliqCL(:,nsrf) = 0. ! eau_liqu integree de couche limite cteiCL(:,nsrf) = 0. ! cloud top instab. crit. couche limite pblt(:,nsrf) = 0. ! T a la Hauteur de couche limite therm(:,nsrf) = 0. trmb1(:,nsrf) = 0. ! deep_cape trmb2(:,nsrf) = 0. ! inhibition trmb3(:,nsrf) = 0. ! Point Omega #undef T2m #define T2m #ifdef T2m ! Calculations of diagnostic t,q at 2m and u, v at 10m ! print*,'Dans pbl OK41' ! print*,'tair1,yt(:,1),y_d_t(:,1)' ! print*, tair1,yt(:,1),y_d_t(:,1) DO j=1, knon i = ni(j) 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) rugo1(j) = yrugos(j) IF(nsrf.EQ.is_oce) THEN rugo1(j) = rugos(i,nsrf) ENDIF psfce(j)=ypaprs(j,1) patm(j)=ypplay(j,1) qairsol(j) = yqsurf(j) END DO ! print*,'Dans pbl OK42A' ! print*,'tair1,yt(:,1),y_d_t(:,1)' ! print*, tair1,yt(:,1),y_d_t(:,1) ! Calculate the temperature et relative humidity at 2m and the wind at 10m CALL stdlevvar(klon, knon, nsrf, zxli, & uzon, vmer, tair1, qair1, zgeo1, & tairsol, qairsol, rugo1, psfce, patm, & yt2m, yq2m, yt10m, yq10m, yu10m, yustar) ! print*,'Dans pbl OK42B' 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))/SQRT(uzon(j)**2+vmer(j)**2) v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)**2+vmer(j)**2) END DO ! print*,'Dans pbl OK43' !IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique !IM Ajoute dependance type surface IF (thermcep) 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) END DO END IF ! print*,'OK pbl 5' 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) 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) END DO ! print*,'OK pbl 6' #else ! T2m not defined ! No calculation PRINT*,' Warning !!! No T2m calculation. Output is set to zero.' #endif !**************************************************************************************** ! 15) End of loop over different surfaces ! !**************************************************************************************** END DO loop_nbsrf !**************************************************************************************** ! 16) Calculate the mean value over all sub-surfaces for som variables ! !**************************************************************************************** ! print*,'OK pbl 7' zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0 zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0 DO nsrf = 1, nbsrf DO k = 1, klev DO i = 1, klon zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf(i,nsrf) zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf(i,nsrf) zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf) zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf) END DO END DO END DO ! print*,'OK pbl 8' DO i = 1, klon zxsens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol zxevap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i) ENDDO ! ! Incrementer la temperature du sol ! zxtsol(:) = 0.0 ; zxfluxlat(:) = 0.0 zt2m(:) = 0.0 ; zq2m(:) = 0.0 zustar(:)=0.0 ; zu10m(:) = 0.0 ; zv10m(:) = 0.0 s_pblh(:) = 0.0 ; s_plcl(:) = 0.0 s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0 s_cteiCL(:) = 0.0; s_pblT(:) = 0.0 s_therm(:) = 0.0 ; s_trmb1(:) = 0.0 s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0 wstar(:,is_ave)=0. ! print*,'OK pbl 9' DO nsrf = 1, nbsrf DO i = 1, klon ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf) wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) & + flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf) wfbilo(i,nsrf) = (evap(i,nsrf) - (rain_f(i) + snow_f(i))) * & pctsrf(i,nsrf) zxtsol(i) = zxtsol(i) + ts(i,nsrf) * pctsrf(i,nsrf) zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf(i,nsrf) zt2m(i) = zt2m(i) + t2m(i,nsrf) * pctsrf(i,nsrf) zq2m(i) = zq2m(i) + q2m(i,nsrf) * pctsrf(i,nsrf) zustar(i) = zustar(i) + ustar(i,nsrf) * pctsrf(i,nsrf) wstar(i,is_ave)=wstar(i,is_ave)+wstar(i,nsrf)*pctsrf(i,nsrf) zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf(i,nsrf) zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf(i,nsrf) s_pblh(i) = s_pblh(i) + pblh(i,nsrf) * pctsrf(i,nsrf) s_plcl(i) = s_plcl(i) + plcl(i,nsrf) * pctsrf(i,nsrf) s_capCL(i) = s_capCL(i) + capCL(i,nsrf) * pctsrf(i,nsrf) s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf(i,nsrf) s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf(i,nsrf) s_pblT(i) = s_pblT(i) + pblT(i,nsrf) * pctsrf(i,nsrf) s_therm(i) = s_therm(i) + therm(i,nsrf) * pctsrf(i,nsrf) s_trmb1(i) = s_trmb1(i) + trmb1(i,nsrf) * pctsrf(i,nsrf) s_trmb2(i) = s_trmb2(i) + trmb2(i,nsrf) * pctsrf(i,nsrf) s_trmb3(i) = s_trmb3(i) + trmb3(i,nsrf) * pctsrf(i,nsrf) END DO END DO ! print*,'OK pbl 10' IF (check) THEN amn=MIN(ts(1,is_ter),1000.) amx=MAX(ts(1,is_ter),-1000.) DO i=2, klon amn=MIN(ts(i,is_ter),amn) amx=MAX(ts(i,is_ter),amx) ENDDO PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx ENDIF !jg ? !!$! !!$! If a sub-surface does not exsist for a grid point, the mean value for all !!$! sub-surfaces is distributed. !!$! !!$ DO nsrf = 1, nbsrf !!$ DO i = 1, klon !!$ IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN !!$ ts(i,nsrf) = zxtsol(i) !!$ t2m(i,nsrf) = zt2m(i) !!$ q2m(i,nsrf) = zq2m(i) !!$ u10m(i,nsrf) = zu10m(i) !!$ v10m(i,nsrf) = zv10m(i) !!$ !!$! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour !!$ pblh(i,nsrf) = s_pblh(i) !!$ plcl(i,nsrf) = s_plcl(i) !!$ capCL(i,nsrf) = s_capCL(i) !!$ oliqCL(i,nsrf) = s_oliqCL(i) !!$ cteiCL(i,nsrf) = s_cteiCL(i) !!$ pblT(i,nsrf) = s_pblT(i) !!$ therm(i,nsrf) = s_therm(i) !!$ trmb1(i,nsrf) = s_trmb1(i) !!$ trmb2(i,nsrf) = s_trmb2(i) !!$ trmb3(i,nsrf) = s_trmb3(i) !!$ ENDIF !!$ ENDDO !!$ ENDDO DO i = 1, klon fder(i) = - 4.0*RSIGMA*zxtsol(i)**3 ENDDO zxqsurf(:) = 0.0 zxsnow(:) = 0.0 DO nsrf = 1, nbsrf DO i = 1, klon zxqsurf(i) = zxqsurf(i) + qsurf(i,nsrf) * pctsrf(i,nsrf) zxsnow(i) = zxsnow(i) + snow(i,nsrf) * pctsrf(i,nsrf) END DO END DO ! Premier niveau de vent sortie dans physiq.F zu1(:) = u(:,1) zv1(:) = v(:,1) ! Some of the module declared variables are returned for printing in physiq.F qsol_d(:) = qsol(:) evap_d(:,:) = evap(:,:) rugos_d(:,:) = rugos(:,:) agesno_d(:,:) = agesno(:,:) END SUBROUTINE pbl_surface ! !**************************************************************************************** ! SUBROUTINE pbl_surface_final(qsol_rst, fder_rst, snow_rst, qsurf_rst, & evap_rst, rugos_rst, agesno_rst, ftsoil_rst) USE indice_sol_mod INCLUDE "dimsoil.h" ! Ouput variables !**************************************************************************************** REAL, DIMENSION(klon), INTENT(OUT) :: qsol_rst REAL, DIMENSION(klon), INTENT(OUT) :: fder_rst REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: snow_rst REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: qsurf_rst REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: evap_rst REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: rugos_rst REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: agesno_rst REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst !**************************************************************************************** ! Return module variables for writing to restart file ! !**************************************************************************************** qsol_rst(:) = qsol(:) fder_rst(:) = fder(:) snow_rst(:,:) = snow(:,:) qsurf_rst(:,:) = qsurf(:,:) evap_rst(:,:) = evap(:,:) rugos_rst(:,:) = rugos(:,:) agesno_rst(:,:) = agesno(:,:) ftsoil_rst(:,:,:) = ftsoil(:,:,:) !**************************************************************************************** ! Deallocate module variables ! !**************************************************************************************** ! DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil) IF (ALLOCATED(qsol)) DEALLOCATE(qsol) IF (ALLOCATED(fder)) DEALLOCATE(fder) IF (ALLOCATED(snow)) DEALLOCATE(snow) IF (ALLOCATED(qsurf)) DEALLOCATE(qsurf) IF (ALLOCATED(evap)) DEALLOCATE(evap) IF (ALLOCATED(rugos)) DEALLOCATE(rugos) IF (ALLOCATED(agesno)) DEALLOCATE(agesno) IF (ALLOCATED(ftsoil)) DEALLOCATE(ftsoil) END SUBROUTINE pbl_surface_final ! !**************************************************************************************** ! SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, tsurf, alb1, alb2, ustar, u10m, v10m, tke) ! Give default values where new fraction has appread USE indice_sol_mod INCLUDE "dimsoil.h" INCLUDE "clesphys.h" INCLUDE "compbl.h" ! Input variables !**************************************************************************************** INTEGER, INTENT(IN) :: itime REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf_new, pctsrf_old ! InOutput variables !**************************************************************************************** REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: tsurf REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: alb1, alb2 REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: ustar,u10m, v10m REAL, DIMENSION(klon,klev+1,nbsrf), INTENT(INOUT) :: tke ! Local variables !**************************************************************************************** INTEGER :: nsrf, nsrf_comp1, nsrf_comp2, nsrf_comp3, i CHARACTER(len=80) :: abort_message CHARACTER(len=20) :: modname = 'pbl_surface_newfrac' INTEGER, DIMENSION(nbsrf) :: nfois=0, mfois=0, pfois=0 ! ! All at once !! !**************************************************************************************** DO nsrf = 1, nbsrf ! First decide complement sub-surfaces SELECT CASE (nsrf) CASE(is_oce) nsrf_comp1=is_sic nsrf_comp2=is_ter nsrf_comp3=is_lic CASE(is_sic) nsrf_comp1=is_oce nsrf_comp2=is_ter nsrf_comp3=is_lic CASE(is_ter) nsrf_comp1=is_lic nsrf_comp2=is_oce nsrf_comp3=is_sic CASE(is_lic) nsrf_comp1=is_ter nsrf_comp2=is_oce nsrf_comp3=is_sic END SELECT ! Initialize all new fractions DO i=1, klon IF (pctsrf_new(i,nsrf) > 0. .AND. pctsrf_old(i,nsrf) == 0.) THEN IF (pctsrf_old(i,nsrf_comp1) > 0.) THEN ! Use the complement sub-surface, keeping the continents unchanged qsurf(i,nsrf) = qsurf(i,nsrf_comp1) evap(i,nsrf) = evap(i,nsrf_comp1) rugos(i,nsrf) = rugos(i,nsrf_comp1) tsurf(i,nsrf) = tsurf(i,nsrf_comp1) alb1(i,nsrf) = alb1(i,nsrf_comp1) alb2(i,nsrf) = alb2(i,nsrf_comp1) ustar(i,nsrf) = ustar(i,nsrf_comp1) u10m(i,nsrf) = u10m(i,nsrf_comp1) v10m(i,nsrf) = v10m(i,nsrf_comp1) if (iflag_pbl > 1) then tke(i,:,nsrf) = tke(i,:,nsrf_comp1) endif mfois(nsrf) = mfois(nsrf) + 1 ELSE ! The continents have changed. The new fraction receives the mean sum of the existent fractions qsurf(i,nsrf) = qsurf(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + qsurf(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3) evap(i,nsrf) = evap(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + evap(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3) rugos(i,nsrf) = rugos(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + rugos(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3) tsurf(i,nsrf) = tsurf(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + tsurf(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3) alb1(i,nsrf) = alb1(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + alb1(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3) alb2(i,nsrf) = alb2(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + alb2(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3) ustar(i,nsrf) = ustar(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + ustar(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3) u10m(i,nsrf) = u10m(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + u10m(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3) v10m(i,nsrf) = v10m(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + v10m(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3) if (iflag_pbl > 1) then tke(i,:,nsrf) = tke(i,:,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + tke(i,:,nsrf_comp3)*pctsrf_old(i,nsrf_comp3) endif ! Security abort. This option has never been tested. To test, comment the following line. ! abort_message='The fraction of the continents have changed!' ! CALL abort_gcm(modname,abort_message,1) nfois(nsrf) = nfois(nsrf) + 1 END IF snow(i,nsrf) = 0. agesno(i,nsrf) = 0. ftsoil(i,:,nsrf) = tsurf(i,nsrf) ELSE pfois(nsrf) = pfois(nsrf)+ 1 END IF END DO END DO END SUBROUTINE pbl_surface_newfrac ! !**************************************************************************************** ! END MODULE pbl_surface_mod