! ============================================================================================================================== ! MODULE : enerbil ! ! CONTACT : orchidee-help _at_ listes.ipsl.fr ! ! LICENCE : IPSL (2006) ! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC ! !>\BRIEF This module computes the energy balance on !! continental surfaces. The module contains the following subroutines: enerbil_initialize, enerbil_main, !! enerbil_clear, enerbil_begin, enerbil_surftemp, enerbil_flux and enerbil_evapveg !! !!\n DESCRIPTION : !! \n !! \latexonly !! \input{enerbil_intro2.tex} !! \endlatexonly !! !! IMPORTANT NOTE: The coefficients A and B are defined differently than those in the referenced !! literature and from those in the code and documentation of the atmospheric model LMDZ. For the !! avoidance of doubt, the coefficients as described here always refer to the ORCHIDEE coefficients, !! and are denoted as such (with the marker: ORC). The re-definition of the coefficients takes place !! within LMDZ before they are passed to ORCHIDEE. The following sequence of expressions is to be !! found within the LMDZ module 'surf_land_orchidee':\n !! !! \latexonly !! \input{surflandLMDZ1.tex} !! \input{surflandLMDZ2.tex} !! \input{surflandLMDZ3.tex} !! \input{surflandLMDZ4.tex} !! \endlatexonly !! \n !! !! \latexonly !! \input{enerbil_symbols.tex} !! \endlatexonly !! !! RECENT CHANGE(S) : None !! !! REFERENCE(S) : None !! !! SVN : !! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/enerbil.f90 $ !! $Date: 2021-07-30 18:36:29 +0200 (Fri, 30 Jul 2021) $ !! $Revision: 7265 $ !! \n !_ ================================================================================================================================ MODULE enerbil ! routines called : restput, restget ! ! modules used USE ioipsl USE xios_orchidee USE ioipsl_para USE constantes USE time, ONLY : one_day, dt_sechiba USE pft_parameters USE qsat_moisture USE sechiba_io_p USE constantes_soil USE explicitsnow IMPLICIT NONE PRIVATE PUBLIC :: enerbil_main, enerbil_initialize, enerbil_finalize, enerbil_clear ! variables used inside enerbil module : declaration and initialisation ! one dimension array allocated, computed and used in enerbil module exclusively REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: psold !! Old surface dry static energy (J kg^{-1}) !$OMP THREADPRIVATE(psold) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: qsol_sat !! Saturated specific humudity for old temperature (kg kg^{-1}) !$OMP THREADPRIVATE(qsol_sat) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: pdqsold !! Deriv. of saturated specific humidity at old temp !! (kg (kg s)^{-1}) !$OMP THREADPRIVATE(pdqsold) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: psnew !! New surface static energy (J kg^{-1}) !$OMP THREADPRIVATE(psnew) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: qsol_sat_new !! New saturated surface air moisture (kg kg^{-1}) !$OMP THREADPRIVATE(qsol_sat_new) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: netrad !! Net radiation (W m^{-2}) !$OMP THREADPRIVATE(netrad) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: lwabs !! LW radiation absorbed by the surface (W m^{-2}) !$OMP THREADPRIVATE(lwabs) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: lwup !! Long-wave up radiation (W m^{-2}) !$OMP THREADPRIVATE(lwup) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: lwnet !! Net Long-wave radiation (W m^{-2}) !$OMP THREADPRIVATE(lwnet) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: fluxsubli !! Energy of sublimation (mm day^{-1}) !$OMP THREADPRIVATE(fluxsubli) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: qsat_air !! Air saturated specific humidity (kg kg^{-1}) !$OMP THREADPRIVATE(qsat_air) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tair !! Air temperature (K) !$OMP THREADPRIVATE(tair) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: q_sol_pot !! Potential surface humidity !$OMP THREADPRIVATE(q_sol_pot) CONTAINS !! ============================================================================================================================= !! SUBROUTINE: enerbil_initialize !! !>\BRIEF Allocate module variables, read from restart file or initialize with default values !! !! DESCRIPTION: The purpose of this module is, firstly, to allocate space !! in memory for key variables within the 'enerbil' module. The second task is to retrieve previous data !! from the restart file. If the variables are not in restart file, default initialization is done. !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): !! !! REFERENCE(S): None !! !! FLOWCHART: None !! \n !_ ============================================================================================================================== SUBROUTINE enerbil_initialize (kjit, kjpindex, index, rest_id, & qair, & temp_sol, temp_sol_new, tsol_rad, & evapot, evapot_corr, qsurf, fluxsens, & fluxlat, vevapp ) !! 0 Variable and parameter description !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-) INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-) INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (-) INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier (-) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity (kg kg^{-1}) !! 0.2 Output variables REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol !! Soil temperature (K) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_new !! New soil temperature (K) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tsol_rad !! Tsol_rad (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: evapot !! Soil Potential Evaporation (mm/tstep) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: evapot_corr !! Soil Potential Evaporation Correction (mm/tstep) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsurf !! Surface specific humidity (kg kg^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsens !! Sensible heat flux (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxlat !! Latent heat flux (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapp !! Total of evaporation (mm day^{-1}) !! 0.4 Local variables INTEGER(i_std) :: ier !_ ================================================================================================================================ IF (printlev>=3) WRITE (numout,*) 'enerbil_initialize : Start initillization' !! 1. Allocate module variables ALLOCATE (psold(kjpindex),stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable ','','') ALLOCATE (qsol_sat(kjpindex),stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable ','','') ALLOCATE (pdqsold(kjpindex),stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable ','','') ALLOCATE (psnew(kjpindex),stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable ','','') ALLOCATE (qsol_sat_new(kjpindex),stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable qsol_sat_new','','') ALLOCATE (netrad(kjpindex),stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable netrad','','') ALLOCATE (lwabs(kjpindex),stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable lwabs','','') ALLOCATE (lwup(kjpindex),stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable lwup','','') ALLOCATE (lwnet(kjpindex),stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable lwnet','','') ALLOCATE (fluxsubli(kjpindex),stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable fluxsubli','','') ALLOCATE (qsat_air(kjpindex),stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable qsat_air','','') ALLOCATE (tair(kjpindex),stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable tair','','') ALLOCATE (q_sol_pot(kjpindex),stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable q_sol_pot','','') !! 2. Initialize variables from restart file or by default values !! The variables read are: temp_sol (surface temperature), qsurf (near surface specific humidity), !! evapot (soil potential evaporation), evapot_corr (corrected soil potential evaporation), tsolrad !! (radiative surface temperature), evapora (evaporation), fluxlat (latent heat flux), fluxsens !! (sensible heat flux) and temp_sol_new (new surface temperature). IF (printlev>=3) WRITE (numout,*) 'Read a restart file for ENERBIL variables' CALL ioconf_setatt_p('UNITS', 'K') CALL ioconf_setatt_p('LONG_NAME','Surface temperature') CALL restget_p (rest_id, 'temp_sol', nbp_glo, 1, 1, kjit, .TRUE., temp_sol, "gather", nbp_glo, index_g) !Config Key = ENERBIL_TSURF !Config Desc = Initial temperature if not found in restart !Config If = OK_SECHIBA !Config Def = 280. !Config Help = The initial value of surface temperature if its value is not found !Config in the restart file. This should only be used if the model is !Config started without a restart file. !Config Units = Kelvin [K] CALL setvar_p (temp_sol, val_exp,'ENERBIL_TSURF', 280._r_std) ! Initialize temp_sol_new with temp_sol. These variables are always equal in the beginning of a new time step. temp_sol_new(:) = temp_sol(:) CALL ioconf_setatt_p('UNITS', 'g/g') CALL ioconf_setatt_p('LONG_NAME','near surface specific humidity') CALL restget_p (rest_id, 'qsurf', nbp_glo, 1, 1, kjit, .TRUE., qsurf, "gather", nbp_glo, index_g) IF ( ALL( qsurf(:) .EQ. val_exp ) ) THEN qsurf(:) = qair(:) ENDIF CALL ioconf_setatt_p('UNITS', 'mm day^{-1}') CALL ioconf_setatt_p('LONG_NAME','Soil Potential Evaporation') CALL restget_p (rest_id, 'evapot', nbp_glo, 1, 1, kjit, .TRUE., evapot, "gather", nbp_glo, index_g) CALL setvar_p (evapot, val_exp, 'ENERBIL_EVAPOT', zero) CALL ioconf_setatt_p('UNITS', 'mm day^{-1}') CALL ioconf_setatt_p('LONG_NAME','Corrected Soil Potential Evaporation') CALL restget_p (rest_id, 'evapot_corr', nbp_glo, 1, 1, kjit, .TRUE., evapot_corr, "gather", nbp_glo, index_g) !Config Key = ENERBIL_EVAPOT !Config Desc = Initial Soil Potential Evaporation !Config If = OK_SECHIBA !Config Def = 0.0 !Config Help = The initial value of soil potential evaporation if its value !Config is not found in the restart file. This should only be used if !Config the model is started without a restart file. !Config Units = CALL setvar_p (evapot_corr, val_exp, 'ENERBIL_EVAPOT', zero) CALL ioconf_setatt_p('UNITS', 'K') CALL ioconf_setatt_p('LONG_NAME','Radiative surface temperature') CALL restget_p (rest_id, 'tsolrad', nbp_glo, 1, 1, kjit, .TRUE., tsol_rad, "gather", nbp_glo, index_g) IF ( ALL( tsol_rad(:) .EQ. val_exp ) ) THEN tsol_rad(:) = temp_sol(:) ENDIF !! 1.3 Set the fluxes so that we have something reasonable and not NaN on some machines CALL ioconf_setatt_p('UNITS', 'Kg/m^2/dt') CALL ioconf_setatt_p('LONG_NAME','Evaporation') CALL restget_p (rest_id, 'evapora', nbp_glo, 1, 1, kjit, .TRUE., vevapp, "gather", nbp_glo, index_g) IF ( ALL( vevapp(:) .EQ. val_exp ) ) THEN vevapp(:) = zero ENDIF CALL ioconf_setatt_p('UNITS', 'W/m^2') CALL ioconf_setatt_p('LONG_NAME','Latent heat flux') CALL restget_p (rest_id, 'fluxlat', nbp_glo, 1, 1, kjit, .TRUE., fluxlat, "gather", nbp_glo, index_g) IF ( ALL( fluxlat(:) .EQ. val_exp ) ) THEN fluxlat(:) = zero ENDIF CALL ioconf_setatt_p('UNITS', 'W/m^2') CALL ioconf_setatt_p('LONG_NAME','Sensible heat flux') CALL restget_p (rest_id, 'fluxsens', nbp_glo, 1, 1, kjit, .TRUE., fluxsens, "gather", nbp_glo, index_g) IF ( ALL( fluxsens(:) .EQ. val_exp ) ) THEN fluxsens(:) = zero ENDIF CALL ioconf_setatt_p('UNITS', 'kg/m^2') CALL ioconf_setatt_p('LONG_NAME','Potential saturated surface humidity') CALL restget_p (rest_id, 'qsolpot', nbp_glo, 1, 1, kjit, .TRUE., q_sol_pot, "gather", nbp_glo, index_g) IF ( ALL( q_sol_pot(:) .EQ. val_exp ) ) THEN q_sol_pot = qsurf ENDIF END SUBROUTINE enerbil_initialize !! =========================================================================================================================== !! SUBROUTINE : enerbil_main !! !! !>\BRIEF Calls each part of the energy budget calculation in sequence !! !! DESCRIPTION : !! This is the main routine for the 'enerbil' module. It is called !! once during the initialisation of ORCHIDEE, and then once for each time step. It is called a final !! time at the culmination of the last time step, for the writing of a restart file.\n !! !! The algorithm first calls 'enerbil_begin' for initialisation, followed by 'enerbil_surftemp' to !! calculate the surface static energy and the saturated surface humidity for the new time step. !! Next is the module 'enerbil_flux' which calculates the new surface temperature, the net radiation in !! the surface layer, the total evaporation and the latent and sensible heat fluxes. Finally comes !! 'enerbil_evapveg', which calculates the evaporation and transpiration from the vegetation.\n !! \n !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S) : evapot, evapot_corr, temp_sol, qsurf, fluxsens, fluxlat, tsol_rad, !! vevapp, temp_sol_new, vevapnu, vevapsno, transpir, vevapwet !! !! REFERENCE(S) : !! - Best, MJ, Beljaars, A, Polcher, J & Viterbo, P, 2004. A proposed structure for coupling tiled !! surfaces with the planetary boundary layer. Journal of Hydrometeorology, 5, pp.1271-1278 !! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations !! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General !! Circulation Model. Journal of Climate, 6, pp.248-273 !! - Dufresne, J-L & Ghattas, J, 2009. Description du schéma de la couche limite turbulente et la !! interface avec la surface planetaire dans LMDZ, Technical note, available (22/12/11): !! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques/ressources/pbl_surface.pdf !! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influences de l'irrigation !! sur le cycle de l'eau, PhD Thesis, available (22/12/11): !! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf !! - Polcher, J. McAvaney, B, Viterbo, P, Gaertner, MA, Hahmann, A, Mahfouf, J-F, Noilhan, J !! Phillips, TJ, Pitman, AJ, Schlosser, CA, Schulz, J-P, Timbal, B, Verseghy, D & !! Xue, Y, 1998. A proposal for a general interface between land surface schemes and !! general circulation models. Global and Planetary Change, 19, pp.261-276 !! - Richtmeyer, RD, Morton, KW, 1967. Difference Methods for Initial-Value Problems. !! Interscience Publishers\n !! - Schulz, Jan-Peter, Lydia Dümenil, Jan Polcher, 2001: On the Land Surface–Atmosphere !! Coupling and Its Impact in a Single-Column Atmospheric Model. J. Appl. Meteor., 40, 642–663. !! !! FLOWCHART : !! \latexonly !! \includegraphics[scale=0.5]{enerbil_main_flowchart.png} !! \endlatexonly !! \n !_ ============================================================================================================================== SUBROUTINE enerbil_main (kjit, kjpindex, & & index, indexveg, zlev, lwdown, swnet, epot_air, temp_air, u, v, petAcoef, petBcoef, & & qair, peqAcoef, peqBcoef, pb, rau, vbeta, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, & & emis, soilflx, soilcap, q_cdrag, humrel, fluxsens, fluxlat, & & vevapp, transpir, transpot, vevapnu, vevapwet, vevapsno, vevapflo, temp_sol, tsol_rad, & & temp_sol_new, qsurf, evapot, evapot_corr, rest_id, hist_id, hist2_id, & & precip_rain, pgflux, snowdz, temp_sol_add) !! 0 Variable and parameter description !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-) INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-) INTEGER(i_std),INTENT (in) :: rest_id,hist_id !! _Restart_ file and _history_ file identifier (-) INTEGER(i_std),INTENT (in) :: hist2_id !! _history_ file 2 identifier (-) INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (-) INTEGER(i_std),DIMENSION(kjpindex*nvm), INTENT(in) :: indexveg !! Indeces of the points on the 3D map REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Down-welling long-wave flux (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: epot_air !! Air potential energy (?? J) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature (K) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petAcoef !! PetAcoef (see note) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petBcoef !! PetBcoef (see note) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity (kg kg^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqAcoef !! PeqAcoef (see note) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqBcoef !! PeqBcoef (see note) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure (hPa) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air density (kg m^{-3}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta !! Resistance coefficient (-) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta1 !! Snow resistance (-) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta4 !! Bare soil resistance (-) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta5 !! Floodplains resistance REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: emis !! Emissivity (-) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilflx !! Effective ground heat flux including both snow and soil (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil calorific capacity including both snow and soil (J K^{-1]) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-) REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (in) :: humrel !! Soil moisture stress (within range 0 to 1) REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta2 !! Interception resistance (-) REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta3 !! Vegetation resistance (-) REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta3pot !! Vegetation resistance for potential transpiration REAL(r_std),DIMENSION (kjpindex),INTENT (in) :: precip_rain !! Rainfall REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowdz !! Snow depth at each snow layer !! 0.2 Output variables REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapnu !! Bare soil evaporation (mm day^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapsno !! Snow evaporation (mm day^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapflo !! Floodplains evaporation REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: transpir !! Transpiration (mm day^{-1}) REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: transpot !! Potential transpiration REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vevapwet !! Interception (mm day^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_new !! New soil temperature (K) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K) !! 0.3 Modified variables REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evapot !! Soil Potential Evaporation (mm/tstep) REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evapot_corr !! Soil Potential Evaporation Correction (mm/tstep) REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: temp_sol !! Soil temperature (K) REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: qsurf !! Surface specific humidity (kg kg^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsens !! Sensible heat flux (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxlat !! Latent heat flux (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: tsol_rad !! Tsol_rad (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapp !! Total of evaporation (mm day^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: pgflux !! Net energy into snowpack(W/m^2) !! 0.4 Local variables REAL(r_std),DIMENSION (kjpindex) :: epot_air_new, qair_new REAL(r_std),DIMENSION (kjpindex) :: diffevap !! Difference betwence vevapp and composing fluxes (Kg/m^2/s) INTEGER(i_std) :: ji,ii,iv !_ ================================================================================================================================ !! 1. Computes some initialisation variables ! Computes some initialisation variables psold (the old surface static energy), qsol_sat ! (the saturated surface humidity) and pdqsold (the derivative of the saturated surface humidity, ! calculated with respect to the surface temperature at the 'old' timestep) CALL enerbil_begin (kjpindex, temp_sol, lwdown, swnet, pb, psold, qsol_sat, pdqsold, netrad, emis) !! 2. Computes psnew (the surface static energy at the 'new' timestep) ! Computes psnew (the surface static energy at the 'new' timestep), qsol_sat_new (the surface ! saturated humidity at the 'new' timestep), temp_sol_new (the surface temperature at the 'new' ! timestep), qair_new (the lowest atmospheric humidity at the 'new' timestep) and epot_air_new ! (the lowest atmospheric evaporation potential at the 'new' timestep) CALL enerbil_surftemp (kjpindex, zlev, emis, & & epot_air, petAcoef, petBcoef, qair, peqAcoef, peqBcoef, soilflx, rau, u, v, q_cdrag, vbeta,& & vbeta1, vbeta5, soilcap, lwdown, swnet, psnew, qsol_sat_new, temp_sol_new, & & qair_new, epot_air_new) !! 3. Diagnose components of the energy budget ! Diagnoses lwup (upwards long wave radiation), lwnet (net long wave radiation), tsol_rad (radiative ! ground temperature), netrad (net radiation), qsurf (surface humidity), vevapp (total evaporation), ! evapot (evaporation potential), evapot_corr (evaporation potential correction factor), ! fluxlat (latent heat flux), fluxsubli (sublimination heat flux) and fluxsens (sensible heat flux). CALL enerbil_flux (kjpindex, emis, temp_sol, rau, u, v, q_cdrag, vbeta, vbeta1, vbeta5, & & qair_new, epot_air_new, psnew, qsurf, & & fluxsens , fluxlat , fluxsubli, vevapp, temp_sol_new, lwdown, swnet, & & lwup, lwnet, pb, tsol_rad, netrad, evapot, evapot_corr, & & precip_rain,snowdz,temp_air,pgflux, soilcap, temp_sol_add) !! 4. Diagnoses the values for evaporation and transpiration ! Diagnoses the values for evaporation and transpiration: vevapsno (snow evaporation), vevapnu ! (bare soil evaporation), transpir (transpiration) and vevapwet (interception) CALL enerbil_evapveg (kjpindex, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, & & rau, u, v, q_cdrag, qair_new, humrel, vevapsno, vevapnu , vevapflo, vevapwet, & & transpir, transpot, evapot) !! 5. Write diagnosics CALL xios_orchidee_send_field("netrad",netrad) CALL xios_orchidee_send_field("evapot",evapot/dt_sechiba) CALL xios_orchidee_send_field("evapot_corr",evapot_corr/dt_sechiba) CALL xios_orchidee_send_field("lwdown",lwabs) CALL xios_orchidee_send_field("lwnet",lwnet) CALL xios_orchidee_send_field("Qv",fluxsubli) DO ji=1,kjpindex diffevap(ji) = vevapp(ji) - ( SUM(vevapwet(ji,:)) + & SUM(transpir(ji,:)) + vevapnu(ji) + vevapsno(ji) + vevapflo(ji) ) ENDDO CALL xios_orchidee_send_field("diffevap",diffevap/dt_sechiba) ! mm/s IF ( .NOT. almaoutput ) THEN CALL histwrite_p(hist_id, 'netrad', kjit, netrad, kjpindex, index) CALL histwrite_p(hist_id, 'evapot', kjit, evapot, kjpindex, index) CALL histwrite_p(hist_id, 'evapot_corr', kjit, evapot_corr, kjpindex, index) CALL histwrite_p(hist_id, 'lwdown', kjit, lwabs, kjpindex, index) CALL histwrite_p(hist_id, 'lwnet', kjit, lwnet, kjpindex, index) IF ( hist2_id > 0 ) THEN CALL histwrite_p(hist2_id, 'netrad', kjit, netrad, kjpindex, index) CALL histwrite_p(hist2_id, 'evapot', kjit, evapot, kjpindex, index) CALL histwrite_p(hist2_id, 'evapot_corr', kjit, evapot_corr, kjpindex, index) CALL histwrite_p(hist2_id, 'lwdown', kjit, lwabs, kjpindex, index) CALL histwrite_p(hist2_id, 'lwnet', kjit, lwnet, kjpindex, index) ENDIF ELSE CALL histwrite_p(hist_id, 'LWnet', kjit, lwnet, kjpindex, index) CALL histwrite_p(hist_id, 'Qv', kjit, fluxsubli, kjpindex, index) CALL histwrite_p(hist_id, 'PotEvap', kjit, evapot_corr, kjpindex, index) CALL histwrite_p(hist_id, 'PotEvapOld', kjit, evapot, kjpindex, index) IF ( hist2_id > 0 ) THEN CALL histwrite_p(hist2_id, 'LWnet', kjit, lwnet, kjpindex, index) CALL histwrite_p(hist2_id, 'Qv', kjit, fluxsubli, kjpindex, index) CALL histwrite_p(hist2_id, 'PotEvap', kjit, evapot_corr, kjpindex, index) ENDIF ENDIF IF (printlev>=3) WRITE (numout,*) ' enerbil_main Done ' END SUBROUTINE enerbil_main !! ============================================================================================================================= !! SUBROUTINE: enerbil_finalize !! !>\BRIEF Write to restart file !! !! DESCRIPTION: This subroutine writes the module variables and variables calculated in enerbil !! to restart file !! !! RECENT CHANGE(S): None !! !! REFERENCE(S): None !! !! FLOWCHART: None !! \n !_ ============================================================================================================================== SUBROUTINE enerbil_finalize (kjit, kjpindex, rest_id, & evapot, evapot_corr, temp_sol, tsol_rad, & qsurf, fluxsens, fluxlat, vevapp ) !! 0 Variable and parameter description !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-) INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-) INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier (-) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation (mm/tstep) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot_corr !! Soil Potential Evaporation Correction (mm/tstep) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Soil temperature (K) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: tsol_rad !! Tsol_rad (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsurf !! Surface specific humidity (kg kg^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: fluxsens !! Sensible heat flux (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: fluxlat !! Latent heat flux (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vevapp !! Total of evaporation (mm day^{-1}) !_ ================================================================================================================================ !! 1. Write variables to restart file to be used for the next simulation IF (printlev>=3) WRITE (numout,*) 'Write restart file with ENERBIL variables' CALL restput_p(rest_id, 'temp_sol', nbp_glo, 1, 1, kjit, temp_sol, 'scatter', nbp_glo, index_g) CALL restput_p(rest_id, 'qsurf', nbp_glo, 1, 1, kjit, qsurf, 'scatter', nbp_glo, index_g) CALL restput_p(rest_id, 'evapot', nbp_glo, 1, 1, kjit, evapot, 'scatter', nbp_glo, index_g) CALL restput_p(rest_id, 'evapot_corr', nbp_glo, 1, 1, kjit, evapot_corr, 'scatter', nbp_glo, index_g) CALL restput_p(rest_id, 'tsolrad', nbp_glo, 1, 1, kjit, tsol_rad, 'scatter', nbp_glo, index_g) CALL restput_p(rest_id, 'evapora', nbp_glo, 1, 1, kjit, vevapp, 'scatter', nbp_glo, index_g) CALL restput_p(rest_id, 'fluxlat', nbp_glo, 1, 1, kjit, fluxlat, 'scatter', nbp_glo, index_g) CALL restput_p(rest_id, 'fluxsens', nbp_glo, 1, 1, kjit, fluxsens, 'scatter', nbp_glo, index_g) CALL restput_p(rest_id, 'qsolpot', nbp_glo, 1, 1, kjit, q_sol_pot, 'scatter', nbp_glo, index_g) END SUBROUTINE enerbil_finalize !! ============================================================================================================================= !! SUBROUTINE : enerbil_clear !! !>\BRIEF Routine deallocates clear output variables if already allocated. !! !! DESCRIPTION : This is a 'housekeeping' routine that deallocates the key output !! variables, if they have already been allocated. The variables that are deallocated are psold, !! qsol_sat, pdqsold, psnew, qsol_sat_new, netrad, lwabs, lwup, lwnet, fluxsubli, qsat_air, tair !! !! RECENT CHANGE(S) : None !! !! MAIN OUTPUT VARIABLE(S) : None !! !! REFERENCES : None !! !! FLOWCHART : None !! \n !_ ============================================================================================================================== SUBROUTINE enerbil_clear () IF ( ALLOCATED (psold)) DEALLOCATE (psold) IF ( ALLOCATED (qsol_sat)) DEALLOCATE (qsol_sat) IF ( ALLOCATED (pdqsold)) DEALLOCATE (pdqsold) IF ( ALLOCATED (psnew)) DEALLOCATE (psnew) IF ( ALLOCATED (qsol_sat_new)) DEALLOCATE (qsol_sat_new) IF ( ALLOCATED (netrad)) DEALLOCATE (netrad) IF ( ALLOCATED (lwabs)) DEALLOCATE (lwabs) IF ( ALLOCATED (lwup)) DEALLOCATE (lwup) IF ( ALLOCATED (lwnet)) DEALLOCATE (lwnet) IF ( ALLOCATED (fluxsubli)) DEALLOCATE (fluxsubli) IF ( ALLOCATED (qsat_air)) DEALLOCATE (qsat_air) IF ( ALLOCATED (tair)) DEALLOCATE (tair) IF ( ALLOCATED (q_sol_pot)) DEALLOCATE (q_sol_pot) END SUBROUTINE enerbil_clear !! ============================================================================================================================= !! SUBROUTINE : enerbil_begin !! !>\BRIEF Preliminary variables required for the calculation of !! the energy budget are derived. !! !! DESCRIPTION : This routines computes preliminary variables required !! for the calculation of the energy budget: the old surface static energy (psold), the surface saturation !! humidity (qsol_sat), the derivative of satured specific humidity at the old temperature (pdqsold) !! and the net radiation (netrad). !! !! MAIN OUTPUT VARIABLE(S) : psold, qsol_sat, pdqsold, netrad !! !! REFERENCE(S) : !! - Best, MJ, Beljaars, A, Polcher, J & Viterbo, P, 2004. A proposed structure for coupling tiled !! surfaces with the planetary boundary layer. Journal of Hydrometeorology, 5, pp.1271-1278 !! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations !! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General !! Circulation Model. Journal of Climate, 6, pp.248-273 !! - Dufresne, J-L & Ghattas, J, 2009. Description du schéma de la couche limite turbulente et la !! interface avec la surface planetaire dans LMDZ, Technical note, available (22/12/11): !! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques/ressources/pbl_surface.pdf !! - Polcher, J. McAvaney, B, Viterbo, P, Gaertner, MA, Hahmann, A, Mahfouf, J-F, Noilhan, J !! Phillips, TJ, Pitman, AJ, Schlosser, CA, Schulz, J-P, Timbal, B, Verseghy, D & !! Xue, Y, 1998. A proposal for a general interface between land surface schemes and !! general circulation models. Global and Planetary Change, 19, pp.261-276 !! - Richtmeyer, RD, Morton, KW, 1967. Difference Methods for Initial-Value Problems. !! Interscience Publishers\n !! - Schulz, Jan-Peter, Lydia Dümenil, Jan Polcher, 2001: On the Land Surface–Atmosphere !! Coupling and Its Impact in a Single-Column Atmospheric Model. J. Appl. Meteor., 40, 642–663. !! !! FLOWCHART : None !! \n !_ ============================================================================================================================== SUBROUTINE enerbil_begin (kjpindex, temp_sol, lwdown, swnet, pb, psold, qsol_sat, pdqsold, netrad, emis) !! 0. Variable and parameter declaration !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Soil temperature (K) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Down-welling long-wave flux (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure (hPa) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: emis !! Emissivity (-) !! 0.2 Output variables REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: psold !! Old surface dry static energy (J kg^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsol_sat !! Saturated specific humudity for old temperature !! (kg kg^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: pdqsold !! Derivative of satured specific humidity at the old !! temperature (kg (kg s)^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: netrad !! Net radiation (W m^{-2}) !! 0.3 Modified variables !! 0.4 Local variables INTEGER(i_std) :: ji REAL(r_std), DIMENSION(kjpindex) :: dev_qsol REAL(r_std), PARAMETER :: missing = 999998. !_ ================================================================================================================================ !! 1. Computes psold (the surface static energy for the old timestep) !! We here define the surface static energy for the 'old' timestep, in terms of the surface !! temperature and heat capacity. !! \latexonly !! \input{enerbilbegin1.tex} !! \endlatexonly psold(:) = temp_sol(:)*cp_air !! 2. Computes qsol_sat (the surface saturated humidity). !! We call the routine 'qsatcalc' from within the module 'src_parameters/constantes_veg'. CALL qsatcalc (kjpindex, temp_sol, pb, qsol_sat) IF ( diag_qsat ) THEN IF ( ANY(ABS(qsol_sat(:)) .GT. missing) ) THEN DO ji = 1, kjpindex IF ( ABS(qsol_sat(ji)) .GT. missing) THEN WRITE(numout,*) 'ERROR on ji = ', ji WRITE(numout,*) 'temp_sol(ji), pb(ji) :', temp_sol(ji), pb(ji) CALL ipslerr_p (3,'enerbil_begin', & & 'qsol too high ','','') ENDIF ENDDO ENDIF ENDIF !! 3. Computes pdqsold !! Computes pdqsold (the derivative of the saturated humidity with respect to temperature !! analysed at the surface temperature at the 'old' timestep. !! We call the routine 'dev_qsatcalc' from qsat_moisture module. CALL dev_qsatcalc (kjpindex, temp_sol, pb, dev_qsol) !! \latexonly !! \input{enerbilbegin2.tex} !! \endlatexonly pdqsold(:) = dev_qsol(:) IF ( diag_qsat ) THEN IF ( ANY(ABS( pdqsold(:)) .GT. missing) ) THEN DO ji = 1, kjpindex IF ( ABS( pdqsold(ji)) .GT. missing ) THEN WRITE(numout,*) 'ERROR on ji = ', ji WRITE(numout,*) 'temp_sol(ji), pb(ji) :', temp_sol(ji), pb(ji) CALL ipslerr_p (3,'enerbil_begin', & & 'pdqsold too high ','','') ENDIF ENDDO ENDIF ENDIF !! 4. Computes the net radiation and the absorbed LW radiation absorbed at the surface. !! Long wave radiation absorbed by the surface is the product of emissivity and downwelling LW radiation !! \latexonly !! \input{enerbilbegin3.tex} !! \endlatexonly lwabs(:) = emis(:) * lwdown(:) !! Net radiation is calculated as: !! \latexonly !! \input{enerbilbegin4.tex} !! \endlatexonly netrad(:) = lwdown(:) + swnet (:) - (emis(:) * c_stefan * temp_sol(:)**4 + (un - emis(:)) * lwdown(:)) IF (printlev>=3) WRITE (numout,*) ' enerbil_begin done ' END SUBROUTINE enerbil_begin !! ============================================================================================================================= !! SUBROUTINE : enerbil_surftemp !! !>\BRIEF This routine computes the new surface static energy !! (psnew) and the saturated humidity at the surface (qsol_sat_new). !! !! DESCRIPTION : This is the key part of the enerbil module, for which !! the energy budget in the surface layer is solved and changes over the model timestep 'dt_sechiba' are !! quantified for the surface static energy, surface temperature, surface humidity, the 'air' (or lowest !! level atmospheric model) temperature, the 'air' (or lowest level atmospheric model) humidity, !! according to the method that is laid out by Dufresne \& Ghattas (2009) and Shultz et al. (2001). !! !! It computes the energy balance at the surface with an implicit scheme, that is connected to the !! Richtmyer and Morton algorithm of the PBL. By computing the surface temperature and surface humidity !! the routine also implicitly estimates the various fluxes which balance in order to give the new !! temperature. Thus once the surface temperature has been obtained all the fluxes need to be diagnosed. !! !! Since the explicit snow scheme is adapted, the notion of skin layer in ORCHIDEE is abandoned and !! the first thin snow layer is adopted to solve the surface energy budget. !! Implicit method is still used for coupling to the atmosphere. !! In this new scheme, the snow temperature profile is simulatenously solved !! along with the surface snow temperature, and the details are referenced to Boone et al. (2010) !! !! MAIN OUTPUT VARIABLE(S) : psnew, qsol_sat_new, temp_sol_new, qair_new, epot_air_new !! !! REFERENCE(S) : !! - Best, MJ, Beljaars, A, Polcher, J & Viterbo, P, 2004. A proposed structure for coupling tiled !! surfaces with the planetary boundary layer. Journal of Hydrometeorology, 5, pp.1271-1278 !! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations !! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General !! Circulation Model. Journal of Climate, 6, pp.248-273 !! - Dufresne, J-L & Ghattas, J, 2009. Description du schéma de la couche limite turbulente et la !! interface avec la surface planetaire dans LMDZ, Technical note, available (22/12/11): !! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques/ressources/pbl_surface.pdf !! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influences de l'irrigation !! sur le cycle de l'eau, PhD Thesis, available (22/12/11): !! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pd !! - Polcher, J. McAvaney, B, Viterbo, P, Gaertner, MA, Hahmann, A, Mahfouf, J-F, Noilhan, J !! Phillips, TJ, Pitman, AJ, Schlosser, CA, Schulz, J-P, Timbal, B, Verseghy, D & !! Xue, Y, 1998. A proposal for a general interface between land surface schemes and !! general circulation models. Global and Planetary Change, 19, pp.261-276 !! - Richtmeyer, RD, Morton, KW, 1967. Difference Methods for Initial-Value Problems. !! Interscience Publishers !! - Schulz, Jan-Peter, Lydia Dümenil, Jan Polcher, 2001: On the Land Surface–Atmosphere !! Coupling and Its Impact in a Single-Column Atmospheric Model. J. Appl. Meteor., 40, 642–663. !! !! FLOWCHART : None !! !_ ============================================================================================================================== SUBROUTINE enerbil_surftemp (kjpindex, zlev, emis, epot_air, & & petAcoef, petBcoef, qair, peqAcoef, peqBcoef, soilflx, rau, u, v, q_cdrag, vbeta,& & vbeta1, vbeta5, soilcap, lwdown, swnet, psnew, qsol_sat_new, temp_sol_new, & & qair_new, epot_air_new) !! 0. Variable and parameter declaration !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: emis !! Emissivity (-) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: epot_air !! Air potential energy (?? J) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petAcoef !! PetAcoef (see note) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petBcoef !! PetBcoef (see note) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity (kg kg^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqAcoef !! PeqAcoef (see note) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqBcoef !! PeqBcoef (see note) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilflx !! Soil flux (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air density (kg m^{-3}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u, v !! Wind velocity by directional components !! u (Eastwards) and v (Northwards) (m s^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta !! Resistance coefficient (-) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta1 !! Snow resistance (-) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta5 !! Floodplains resistance REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil calorific capacity (J K^{-1]) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Down-welling long-wave flux (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux (W m^{-2}) !! 0.2 Output variables REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: psnew !! New surface static energy (J kg^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsol_sat_new !! New saturated surface air moisture (kg kg^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_new !! New soil temperature (K) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qair_new !! New air moisture (kg kg^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: epot_air_new !! New air temperature (K) !! 0.4 Local variables INTEGER(i_std) :: ji REAL(r_std),DIMENSION (kjpindex) :: zicp REAL(r_std) :: fevap REAL(r_std) :: sensfl_old, larsub_old, lareva_old, dtheta, sum_old, sum_sns REAL(r_std) :: zikt, zikq, netrad_sns, sensfl_sns, larsub_sns, lareva_sns REAL(r_std) :: speed REAL(r_std),DIMENSION (kjpindex) :: dtes !! Diagnostic output variable for the change of the thermal energy content of surface for which the energy balance is calculated (J m-2) !_ ================================================================================================================================ zicp = un / cp_air DO ji=1,kjpindex !! 1. Derivation of auxiliary variables !! \latexonly !! \input{enerbilsurftemp1.tex} !! \endlatexonly speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji))) ! !! \latexonly !! \input{enerbilsurftemp2.tex} !! \endlatexonly zikt = 1/(rau(ji) * speed * q_cdrag(ji)) zikq = 1/(rau(ji) * speed * q_cdrag(ji)) !! 2. Computation of fluxes for the old surface conditions !! As a reference, we first calculate the sensible and latent heat for the 'old' timestep. !! 2.1 Sensible heat for the 'old' timestep !! This is equation (64) of (Dufresne & Ghattas, 2009), rewritten in terms of H_old. We make the !! approximation that (P_r/P)^{kappa}=1 and convert the A and B coefficients to the ORCHIDEE !! format (see introductory note, above). Also, the equation here is in terms of surface static !! energy, which is defined as per section 3.1, above. !! \latexonly !! \input{enerbilsurftemp3.tex} !! \endlatexonly sensfl_old = (petBcoef(ji) - psold(ji)) / (zikt - petAcoef(ji)) !! 2.2 Latent heat for the 'old' timestep !! This is equation (70) of (Dufresne & Ghattas, 2009), rewritten in terms of {\lambda E}_{old}. !! Again we convert the A and B coefficients from the LMDZ format to the ORCHIDEE format.\n !! There are two forms of the equation - first for the latent heat as a result of sublimination !! processes and then as a result of evaporative processes: !! \latexonly !! \input{enerbilsurftemp4.tex} !! \endlatexonly larsub_old = chalsu0 * vbeta1(ji) * (un - vbeta5(ji)) * & (peqBcoef(ji) - qsol_sat(ji)) / (zikq - vbeta1(ji) * (un - vbeta5(ji))* peqAcoef(ji)) !! \latexonly !! \input{enerbilsurftemp5.tex} !! \endlatexonly lareva_old = chalev0 * (un - vbeta1(ji)) * (un - vbeta5(ji)) * vbeta(ji) * & (peqBcoef(ji) - qsol_sat(ji)) / & (zikq - (un - vbeta1(ji)) * (un - vbeta5(ji)) * vbeta(ji) * peqAcoef(ji)) & + chalev0 * vbeta5(ji) * (peqBcoef(ji) - qsol_sat(ji)) / (zikq - vbeta5(ji) * peqAcoef(ji)) !! 3. Calculation of sensitivity terms !! We here calculate the sensitivity for the different fluxes to changes in dtheta, which is the !! change in the surface static energy over the model time step (dt_sechiba). !! 3.1 Net radiation sensitivity !! This is an explicit-implicit representation of the Stefan-Boltzmann law - the explicit terms !! are ${ps}_{old}^3$, and it is completed when multiplied by dtheta (defined above). !! \latexonly !! \input{enerbilsurftemp6.tex} !! \endlatexonly netrad_sns = zicp(ji) * quatre * emis(ji) * c_stefan * ((zicp(ji) * psold(ji))**3) !! 3.2 Sensible heat flux sensitivity !! This is the temperature sensitivity term N_1^h derived in equation (66) of (Dufresne & Ghattas, 2009), !! where we again assume that (P_r/P)^{kappa}=1, convert the A and B coefficients to the ORCHIDEE format, !! and rewrite in terms of surface static energy, rather than temperature (see section 3.1, above). !! \latexonly !! \input{enerbilsurftemp7.tex} !! \endlatexonly sensfl_sns = un / (zikt - petAcoef(ji)) !! 3.3 Latent heat flux sensitivity !! This is the humidity sensitivity term N_1^q derived in equation (72) of (Dufresne & Ghattas, 2009), where !! the A and B coefficients are written in the ORCHIDEE format. !! larsub_sns is the latent heat sensitivity for sublimination. The coefficient vbeta1 is the evaporation !! coefficient. It represents the relationship between the real and potential evaporation. It is derived !! in the module 'src_sechiba/diffuco_snow'. !! \latexonly !! \input{enerbilsurftemp8.tex} !! \endlatexonly larsub_sns = chalsu0 * vbeta1(ji) * (un - vbeta5(ji)) * zicp(ji) * & pdqsold(ji) / (zikq - vbeta1(ji) * (un - vbeta5(ji)) * peqAcoef(ji)) !! lareva_sns is the latent heat sensitivity for evaporation. vbeta1 (src_sechiba/diffuco_snow), !! and vbeta (src_sechiba/diffuco_comb) are the evaporation !! coefficients. !! \latexonly !! \input{enerbilsurftemp9.tex} !! \endlatexonly lareva_sns = chalev0 * ((un - vbeta1(ji))*(un - vbeta5(ji)) * vbeta(ji) + vbeta5(ji)) * & & zicp(ji) * pdqsold(ji) / & (zikq - ((un - vbeta1(ji))*(un - vbeta5(ji)) * vbeta(ji) + vbeta5(ji))* peqAcoef(ji)) !! 4. Solution of the energy balance !! 4.1 We calculate the total flux for the 'old' timestep. !! \latexonly !! \input{enerbilsurftemp10.tex} !! \endlatexonly sum_old = netrad(ji) + sensfl_old + larsub_old + lareva_old + soilflx(ji) !! 4.2 We calculate the total sensitivity dtheta (the change in the !! surface static energy over the timestep. !! \latexonly !! \input{enerbilsurftemp11.tex} !! \endlatexonly sum_sns = netrad_sns + sensfl_sns + larsub_sns + lareva_sns !! 4.3 Calculation of dtheta (change in surface static energy over the !! timestep. !! \latexonly !! \input{enerbilsurftemp12.tex} !! \endlatexonly dtheta = dt_sechiba * sum_old / (zicp(ji) * soilcap(ji) + dt_sechiba * sum_sns) !! 4.4 Determination of state variables for the 'new' timestep !! No we have dtheta, we can determine the surface static energy that !! corresponds to the 'new' timestep. !! \latexonly !! \input{enerbilsurftemp13.tex} !! \endlatexonly psnew(ji) = psold(ji) + dtheta !! The new surface saturated humidity can be calculated by equation (69) !! of (Dufresne & Ghattas, 2009), in which we substitute dtheta for the !! change between old and new temperature using the relationship from 3.1, !! above. !! \latexonly !! \input{enerbilsurftemp14.tex} !! \endlatexonly qsol_sat_new(ji) = qsol_sat(ji) + zicp(ji) * pdqsold(ji) * dtheta !! The new surface temperature is determined from the new surface static temperature, !! again by using the relationship in section 3.1. !! \latexonly !! \input{enerbilsurftemp15.tex} !! \endlatexonly temp_sol_new(ji) = psnew(ji) / cp_air !! 4.5 Calculation of new evaporation potential and new evaporation latent heat !! flux (???) !! \latexonly !! \input{enerbilsurftemp16.tex} !! \endlatexonly epot_air_new(ji) = zikt * (sensfl_old - sensfl_sns * dtheta) + psnew(ji) !! \latexonly !! \input{enerbilsurftemp17.tex} !! \endlatexonly fevap = (lareva_old - lareva_sns * dtheta) + (larsub_old - larsub_sns * dtheta) IF ( ABS(fevap) < EPSILON(un) ) THEN !! \latexonly !! \input{enerbilsurftemp18.tex} !! \endlatexonly qair_new(ji) = qair(ji) ELSE !! \latexonly !! \input{enerbilsurftemp19.tex} !! \endlatexonly qair_new(ji) = zikq * un / ( chalsu0 * vbeta1(ji) * (un - vbeta5(ji)) + & & chalev0 * ((un - vbeta1(ji))*(un - vbeta5(ji)) * vbeta(ji) + vbeta5(ji)) ) & & * fevap + qsol_sat_new(ji) ENDIF ! Calculate dtes only for diagnostic output and send it to xios dtes(ji) = (soilcap(ji)/cp_air)*dtheta/dt_sechiba ENDDO ! Send diagnostic variable to XIOS CALL xios_orchidee_send_field("dtes", dtes) IF (printlev>=3) WRITE (numout,*) ' enerbil_surftemp done ' END SUBROUTINE enerbil_surftemp !! ============================================================================================================================= !! SUBROUTINE : enerbil_pottemp !! !>\BRIEF This subroutine computes the surface temperature and humidity should the surface been unstressed !! !! DESCRIPTION : !! !! MAIN OUTPUT VARIABLE(S) : !! !! REFERENCE(S) : !! !! FLOWCHART : None !! !_ ============================================================================================================================== SUBROUTINE enerbil_pottemp (kjpindex, zlev, emis, epot_air, & & petAcoef, petBcoef, qair, peqAcoef, peqBcoef, soilflx, rau, u, v, q_cdrag, vbeta,& & vbeta1, vbeta5, soilcap, lwdown, swnet, q_sol_pot, temp_sol_pot) ! interface ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: emis !! Emissivity REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: epot_air !! Air potential energy REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petAcoef !! PetAcoef REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petBcoef !! PetBcoef REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqAcoef !! PeqAcoef REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqBcoef !! PeqBcoef REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilflx !! Soil flux REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u, v !! Wind REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta !! Resistance coefficient REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta1 !! Snow resistance REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta5 !! Floodplains resistance REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil calorific capacity REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Down-welling long-wave flux REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux ! output fields REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: q_sol_pot !! Potential surface air moisture REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_pot !! Potential soil temperature ! local declaration INTEGER(i_std) :: ji REAL(r_std),DIMENSION (kjpindex) :: zicp REAL(r_std) :: fevap REAL(r_std) :: sensfl_old, larsub_old, lareva_old, dtheta, sum_old, sum_sns REAL(r_std) :: zikt, zikq, netrad_sns, sensfl_sns, larsub_sns, lareva_sns REAL(r_std) :: speed !_ ============================================================================================================================== zicp = un / cp_air ! DO ji=1,kjpindex dtheta = zero fevap = zero temp_sol_pot(ji) = temp_sol_pot(ji) + dtheta q_sol_pot(ji) = q_sol_pot(ji) + fevap ENDDO END SUBROUTINE enerbil_pottemp !! ============================================================================================================================= !! SUBROUTINE : enerbil_flux !! !>\BRIEF Computes the new soil temperature, net radiation and !! latent and sensible heat flux for the new time step. !! !! DESCRIPTION : This routine diagnoses, based on the new soil temperature, !! the net radiation, the total evaporation, the latent heat flux, the sensible heat flux and the !! sublimination flux. It also diagnoses the potential evaporation used for the fluxes (evapot) and the potential !! as defined by Penman & Monteith (Monteith, 1965) based on the correction term developed by Chris !! Milly (1992). This Penman-Monteith formulation is required for the estimation of bare soil evaporation !! when the 11 layer CWRR moisture scheme is used for the soil hydrology. !! !! MAIN OUTPUT VARIABLE(S) : qsurf, fluxsens, fluxlat, fluxsubli, vevapp, lwup, lwnet, !! tsol_rad, netrad, evapot, evapot_corr !! !! REFERENCE(S) : !! - Best, MJ, Beljaars, A, Polcher, J & Viterbo, P, 2004. A proposed structure for coupling tiled !! surfaces with the planetary boundary layer. Journal of Hydrometeorology, 5, pp.1271-1278 !! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations !! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General !! Circulation Model. Journal of Climate, 6, pp.248-273 !! - Dufresne, J-L & Ghattas, J, 2009. Description du schéma de la couche limite turbulente et la !! interface avec la surface planetaire dans LMDZ, Technical note, available (22/12/11): !! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques/ressources/pbl_surface.pdf !! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influences de l'irrigation !! sur le cycle de l'eau, PhD Thesis, available (22/12/11): !! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf !! - Monteith, JL, 1965. Evaporation and Environment, paper presented at Symposium of the Society !! for Experimental Biology !! - Monteith & Unsworth, 2008. Principles of Environmental Physics (third edition), published Elsevier !! ISBN 978-0-12-505103-3 !! - Milly, P. C. D., 1992: Potential Evaporation and Soil Moisture in General Circulation Models. !! Journal of Climate, 5, pp. 209–226. !! - Polcher, J. McAvaney, B, Viterbo, P, Gaertner, MA, Hahmann, A, Mahfouf, J-F, Noilhan, J !! Phillips, TJ, Pitman, AJ, Schlosser, CA, Schulz, J-P, Timbal, B, Verseghy, D & !! Xue, Y, 1998. A proposal for a general interface between land surface schemes and !! general circulation models. Global and Planetary Change, 19, pp.261-276 !! - Richtmeyer, RD, Morton, KW, 1967. Difference Methods for Initial-Value Problems. !! Interscience Publishers !! - Schulz, Jan-Peter, Lydia Dümenil, Jan Polcher, 2001: On the Land Surface–Atmosphere !! Coupling and Its Impact in a Single-Column Atmospheric Model. J. Appl. Meteor., 40, 642–663. !! !! FLOWCHART : None !! \n !_ ============================================================================================================================== SUBROUTINE enerbil_flux (kjpindex, emis, temp_sol, rau, u, v, q_cdrag, vbeta, vbeta1, vbeta5, & & qair, epot_air, psnew, qsurf, fluxsens, fluxlat, fluxsubli, vevapp, temp_sol_new, & & lwdown, swnet, lwup, lwnet, pb, tsol_rad, netrad, evapot, evapot_corr,& & precip_rain,snowdz,temp_air,pgflux, soilcap, temp_sol_add) !! 0. Variable and parameter declaration !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: emis !! Emissivity (-) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Surface temperature (K) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density (kg m^{-3}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u,v !! Wind velocity in components u and v (m s^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta !! Resistance coefficient (-) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta1 !! Snow resistance (-) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta5 !! Flood resistance REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity (kg kg^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: epot_air !! Air potential energy (J) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: psnew !! New surface static energy (J kg^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature (K) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure (hPa) REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowdz !! Snow depth REAL(r_std), DIMENSION (kjpindex),INTENT(in) :: precip_rain !! Rainfall REAL(r_std), DIMENSION (kjpindex),INTENT(in) :: temp_air !! Air temperature REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Downward long wave radiation (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net short wave radiation (W m^{-2}) REAL(r_std), DIMENSION (kjpindex),INTENT(in) :: soilcap !! Soil calorific capacity including snow and soil (J K^{-1}) !! 0.2 Output variables REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsurf !! Surface specific humidity (kg kg^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsens !! Sensible heat flux (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxlat !! Latent heat flux (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsubli !! Energy of sublimation (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapp !! Total of evaporation (mm day^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: lwup !! Long-wave up radiation (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: lwnet !! Long-wave net radiation (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tsol_rad !! Radiative surface temperature (W m^{-2}) REAL(r_std), DIMENSION (kjpindex),INTENT(out) :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K) !! 0.3 Modified variables REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: netrad !! Net radiation (W m^{-2}) REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evapot !! Soil Potential Evaporation (mm/tstep) REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evapot_corr !! Soil Potential Evaporation Correction (mm/tstep) REAL(r_std), DIMENSION (kjpindex),INTENT(inout) :: pgflux !! Net energy into snowpack(W/m^2) !! 0.4 Local variables INTEGER(i_std) :: ji REAL(r_std),DIMENSION (kjpindex) :: grad_qsat REAL(r_std) :: correction REAL(r_std) :: speed !! Speed (m s^{-1}), REAL(r_std) :: qc !! Surface drag coefficient (??) LOGICAL,DIMENSION (kjpindex) :: warning_correction REAL(r_std), DIMENSION (kjpindex) :: zgflux REAL(r_std), DIMENSION (kjpindex) :: qsol_sat_tmp REAL(r_std), DIMENSION (kjpindex) :: zerocelcius REAL(r_std), DIMENSION (kjpindex) :: PHPSNOW REAL(r_std), DIMENSION (kjpindex) :: lwup_tmp REAL(r_std), DIMENSION (kjpindex) :: netrad_tmp REAL(r_std), DIMENSION (kjpindex) :: fluxsens_tmp REAL(r_std), DIMENSION (kjpindex) :: fluxlat_tmp REAL(r_std), DIMENSION (kjpindex) :: betatot !! beta used in the computation of qsurf (-) !_ ================================================================================================================================ zerocelcius(:) = tp_00 CALL qsatcalc (kjpindex, zerocelcius, pb, qsol_sat_tmp) DO ji=1,kjpindex !! 1. Determination of 'housekeeping' variables !! The horizontal wind speed is calculated from the velocities along each axis using Pythagorus' theorem. !! \latexonly !! \input{enerbilflux1.tex} !! \endlatexonly speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji))) !! From definition, the surface drag coefficient is defined: !! \latexonly !! \input{enerbilflux2.tex} !! \endlatexonly qc = speed * q_cdrag(ji) !! 2. Calculation of the net upward flux of longwave radiation !! We first of all calculate the radiation as a result of the Stefan-Boltzmann equation, !! which is the sum of the calculated values at the surface temperature at the 'old' !! temperature and the value that corresponds to the difference between the old temperature !! and the temperature at the 'new' timestep. !! \latexonly !! \input{enerbilflux3.tex} !! \endlatexonly lwup(ji) = emis(ji) * c_stefan * temp_sol(ji)**4 + & & quatre * emis(ji) * c_stefan * temp_sol(ji)**3 * & & (temp_sol_new(ji) - temp_sol(ji)) !! We then add to this value that from the reflected longwave radiation: !! \latexonly !! \input{enerbilflux4.tex} !! \endlatexonly lwup(ji) = lwup(ji) + (un - emis(ji)) * lwdown(ji) !! The radiative surface temperature is calculated by inverting the Stefan-Boltzmann relation: !! \latexonly !! \input{enerbilflux5.tex} !! \endlatexonly !! The implicit solution computes an emitted long-wave flux which is a limited Taylor expansion !! of the future surface temperature around the current values. Thus the long-wave flux does not !! correspond to the new surface temperature but some intermediate value. So we need to deduce the !! radiative surface temperature to which this flux corresponds. !! tsol_rad(ji) = (lwup(ji)/ (emis(ji) * c_stefan)) **(1./quatre) !! qsurf (the surface specific humidity) is a simple diagnostic which will be used by the !! GCM to compute the dependence of of the surface layer stability on moisture. !! \latexonly !! \input{enerbilflux6.tex} !! \endlatexonly betatot(ji) = vbeta1(ji) * (un - vbeta5(ji)) + vbeta5(ji) + & & (un - vbeta1(ji))*(un - vbeta5(ji)) * vbeta(ji) qsurf(ji) = betatot(ji) * qsol_sat_new(ji) + (un - betatot(ji)) * qair(ji) ! When no turbulence and no wind (qc small) the near surface air moisture is equal to qair. ! This is related to problem with interpolation in stdlevvar in LMDZ. IF ( qc .LE. min_qc ) qsurf(ji) = qair(ji) !! \latexonly !! \input{enerbilflux7.tex} !! \endlatexonly qsurf(ji) = MAX(qsurf(ji), qair(ji)) !! Net downward radiation is the sum of the down-welling less the up-welling long wave flux, plus !! the short wave radiation. !! \latexonly !! \emissivity absorbed lw radiationinput{enerbilflux8.tex} !! \endlatexonly netrad(ji) = lwdown(ji) + swnet(ji) - lwup(ji) !! 'vevapp' is the sum of the total evaporative processes (snow plus non-snow processes). !! \latexonly !! \input{enerbilflux9.tex} !! \endlatexonly vevapp(ji) = dt_sechiba * rau(ji) * qc * (vbeta1(ji) * (un - vbeta5(ji)) + vbeta5(ji)) * & & (qsol_sat_new(ji) - qair(ji)) + & & dt_sechiba * rau(ji) * qc * (un - vbeta1(ji))*(un-vbeta5(ji)) * vbeta(ji) * & & (qsol_sat_new(ji) - qair(ji)) !! The total latent heat flux is the sum of the snow plus non-snow processes. !! \latexonly !! \input{enerbilflux10.tex} !! \endlatexonly fluxlat(ji) = chalsu0 * rau(ji) * qc * vbeta1(ji) * (un - vbeta5(ji)) * & & (qsol_sat_new(ji) - qair(ji)) + & & chalev0 * rau(ji) * qc * vbeta5(ji) *& & (qsol_sat_new(ji) - qair(ji)) + & & chalev0 * rau(ji) * qc * (un - vbeta1(ji)) * (un - vbeta5(ji)) * vbeta(ji) * & & (qsol_sat_new(ji) - qair(ji)) !! The sublimination flux concerns is calculated using vbeta1, the snow resistance. !! \latexonly !! \input{enerbilflux11.tex} !! \endlatexonly fluxsubli(ji) = chalsu0 * rau(ji) * qc * vbeta1(ji) * (un - vbeta5(ji)) * & & (qsol_sat_new(ji) - qair(ji)) !! The sensible heat flux is a factor of the difference between the new surface static energy !! and the potential energy of air. !! \latexonly !! \input{enerbilflux12.tex} !! \endlatexonly fluxsens(ji) = rau(ji) * qc * (psnew(ji) - epot_air(ji)) !! This is the net longwave downwards radiation. !! \latexonly !! \input{enerbilflux13.tex} !! \endlatexonly lwnet(ji) = lwdown(ji) - lwup(ji) !! Diagnoses the potential evaporation used for the fluxes (evapot) !! \latexonly !! \input{enerbilflux14.tex} !! \endlatexonly evapot(ji) = MAX(zero, dt_sechiba * rau(ji) * qc * (qsol_sat_new(ji) - qair(ji))) !! From definition we can say: !! \latexonly !! \input{enerbilflux15.tex} !! \endlatexonly tair(ji) = epot_air(ji) / cp_air !! To calculate net energy flux into the snowpack PHPSNOW(ji) = precip_rain(ji)*(4.218E+3)*(MAX(tp_00,temp_air(ji))-tp_00)/dt_sechiba ! (w/m2) pgflux(ji) = netrad(ji) - fluxsens(ji) - fluxlat(ji) + PHPSNOW(ji) !! To get the extra energy used to melt the snowpack IF (temp_sol_new (ji) > tp_00 .AND. & SUM(snowdz(ji,:)) .GT. zero .AND. soilcap(ji) .GT. min_sechiba) THEN lwup_tmp(ji) = emis(ji) * c_stefan * temp_sol(ji)**4 + & quatre * emis(ji) * c_stefan * temp_sol(ji)**3 * (tp_00 - temp_sol(ji)) lwup_tmp(ji) = lwup_tmp(ji) + (un - emis(ji)) * lwdown(ji) netrad_tmp(ji) = lwdown(ji) + swnet(ji) - lwup_tmp(ji) fluxsens_tmp(ji) = rau(ji) * qc * cp_air * (tp_00 - epot_air(ji)/cp_air) fluxlat_tmp(ji) = chalsu0 * rau(ji) * qc * vbeta1(ji) * (un-vbeta5(ji)) * & & (qsol_sat_tmp(ji) - qair(ji)) + & & chalev0 * rau(ji) * qc * vbeta5(ji) *& & (qsol_sat_tmp(ji) - qair(ji)) + & & chalev0 * rau(ji) * qc * (un - vbeta1(ji)) * (un-vbeta5(ji))* vbeta(ji) * & & (qsol_sat_tmp(ji) - qair(ji)) zgflux(ji) = netrad_tmp(ji) - fluxsens_tmp(ji) - fluxlat_tmp(ji)+PHPSNOW(ji) temp_sol_add(ji) = -(pgflux(ji) - zgflux(ji))*dt_sechiba/soilcap(ji) pgflux(ji) = zgflux(ji) ELSE temp_sol_add(ji) = zero ENDIF ENDDO !! 3. Define qsat_air with the subroutine src_parameter: CALL qsatcalc(kjpindex, tair, pb, qsat_air) CALL dev_qsatcalc(kjpindex, tair, pb, grad_qsat) ! grad_qsat(:)= (qsol_sat_new(:)- qsat_air(:)) / ((psnew(:) - epot_air(:)) / cp_air) ! * dt_sechiba warning_correction(:)=.FALSE. DO ji=1,kjpindex !! \latexonly !! \input{enerbilflux16.tex} !! \endlatexonly speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji))) !! \latexonly !! \input{enerbilflux17.tex} !! \endlatexonly qc = speed * q_cdrag(ji) !! Derive the potential as defined by Penman & Monteith (Monteith, 1965) based on the correction !! term developed by Chris Milly (1992). IF ((evapot(ji) .GT. min_sechiba) .AND. ((psnew(ji) - epot_air(ji)) .NE. zero )) THEN ! !! \latexonly !! \input{enerbilflux18.tex} !! \endlatexonly correction = (quatre * emis(ji) * c_stefan * tair(ji)**3 + rau(ji) * qc * cp_air + & & chalev0 * rau(ji) * qc * grad_qsat(ji) * vevapp(ji) / evapot(ji) ) !! \latexonly !! \input{enerbilflux19.tex} !! \endlatexonly IF (ABS(correction) .GT. min_sechiba) THEN correction = chalev0 * rau(ji) * qc * grad_qsat(ji) * (un - vevapp(ji)/evapot(ji)) / correction ELSE warning_correction(ji)=.TRUE. ENDIF ELSE correction = zero ENDIF correction = MAX (zero, correction) !! \latexonly !! \input{enerbilflux20.tex} !! \endlatexonly evapot_corr(ji) = evapot(ji) / (un + correction) ENDDO IF ( ANY(warning_correction) ) THEN DO ji=1,kjpindex IF ( warning_correction(ji) ) THEN WRITE(numout,*) ji,"Denominateur de la correction de milly nul! Aucune correction appliquee" ENDIF ENDDO ENDIF IF (printlev>=3) WRITE (numout,*) ' enerbil_flux done ' END SUBROUTINE enerbil_flux !! =========================================================================================================================== !! SUBROUTINE : enerbil_evapveg !! !>\BRIEF Splits total evaporation loss into individual process !! components and calculates the assimilation. !! !! DESCRIPTION : Based on the estimation of the fluxes in enerbil_flux, !! this routine splits the total evaporation into transpiration interception loss from vegetation, !! bare soil evaporation and snow sublimation. It then calculates the assimilation. !! !! MAIN OUTPUT VARIABLE(S) : vevapsno, vevapnu, transpir, vevapwet !! !! REFERENCE(S) : !! - Best, MJ, Beljaars, A, Polcher, J & Viterbo, P, 2004. A proposed structure for coupling tiled !! surfaces with the planetary boundary layer. Journal of Hydrometeorology, 5, pp.1271-1278 !! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations !! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General !! Circulation Model. Journal of Climate, 6, pp.248-273 !! - Dufresne, J-L & Ghattas, J, 2009. Description du schéma de la couche limite turbulente et la !! interface avec la surface planetaire dans LMDZ, Technical note, available (22/12/11): !! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques/ressources/pbl_surface.pdf !! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influences de l'irrigation !! sur le cycle de l'eau, PhD Thesis, available (22/12/11): !! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf !! - Polcher, J. McAvaney, B, Viterbo, P, Gaertner, MA, Hahmann, A, Mahfouf, J-F, Noilhan, J !! Phillips, TJ, Pitman, AJ, Schlosser, CA, Schulz, J-P, Timbal, B, Verseghy, D & !! Xue, Y, 1998. A proposal for a general interface between land surface schemes and !! general circulation models. Global and Planetary Change, 19, pp.261-276 !! - Richtmeyer, RD, Morton, KW, 1967. Difference Methods for Initial-Value Problems. !! Interscience Publishers !! - Schulz, Jan-Peter, Lydia Dümenil, Jan Polcher, 2001: On the Land Surface–Atmosphere !! Coupling and Its Impact in a Single-Column Atmospheric Model. J. Appl. Meteor., 40, 642–663. !! !! FLOWCHART : None !! \n !_ ============================================================================================================================== SUBROUTINE enerbil_evapveg (kjpindex, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, & & rau, u, v, q_cdrag, qair, humrel, vevapsno, vevapnu , vevapflo, vevapwet, & & transpir, transpot, evapot) !! 0. Variable and parameter declaration !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta1 !! Snow resistance (-) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta4 !! Bare soil resistance (-) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta5 !! Floodplains resistance REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density (kg m^{-3}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u, v !! Wind velocity in directions u and v (m s^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity (kg kg^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation (mm/tstep) REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (in) :: humrel !! Soil moisture stress (within range 0 to 1) !!$ DS 15022011 humrel was used in a previous version of Orchidee, developped by Nathalie. Need to be discussed if it should be introduce again REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta2 !! Interception resistance (-) REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta3 !! Vegetation resistance (-) REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta3pot !! Vegetation resistance for potential transpiration !! 0.2 Output variables REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapsno !! Snow evaporation (mm day^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapnu !! Bare soil evaporation (mm day^{-1}) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapflo !! Floodplains evaporation REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: transpir !! Transpiration (mm day^{-1}) REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: transpot !! Potential Transpiration REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vevapwet !! Interception (mm day^{-1}) !! 0.3 Modified variables !! 0.4 Local variables INTEGER(i_std) :: ji, jv REAL(r_std), DIMENSION(kjpindex) :: xx REAL(r_std), DIMENSION(kjpindex) :: vbeta2sum, vbeta3sum REAL(r_std) :: speed !! Speed (m s^{-1}) !_ ============================================================================================================================== ! initialisation: utile pour calculer l'evaporation des floodplains dans lesquelles il y a de la vegetation vbeta2sum(:) = 0. vbeta3sum(:) = 0. DO jv = 1, nvm vbeta2sum(:) = vbeta2sum(:) + vbeta2(:,jv) vbeta3sum(:) = vbeta3sum(:) + vbeta3(:,jv) ENDDO !! 1. Compute vevapsno (evaporation from snow) and vevapnu (bare soil evaporation) DO ji=1,kjpindex !! \latexonly !! \input{enerbilevapveg1.tex} !! \endlatexonly speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji))) !! 1.1 Snow sublimation !! \latexonly !! \input{enerbilevapveg2.tex} !! \endlatexonly vevapsno(ji) = (un - vbeta5(ji)) * vbeta1(ji) * dt_sechiba * rau(ji) * speed * q_cdrag(ji) * (qsol_sat_new(ji) - qair(ji)) !! 1.2 Bare soil evaporation !! \latexonly !! \input{enerbilevapveg3.tex} !! \endlatexonly vevapnu(ji) = (un - vbeta1(ji)) * (un-vbeta5(ji)) * vbeta4(ji) * dt_sechiba * rau(ji) * speed * q_cdrag(ji) & & * (qsol_sat_new(ji) - qair(ji)) ! ! 1.3 floodplains evaporation - transpiration et interception prioritaires dans les floodplains ! vevapflo(ji) = vbeta5(ji) & & * dt_sechiba * rau(ji) * speed * q_cdrag(ji) * (qsol_sat_new(ji) - qair(ji)) END DO !! 2. Compute transpir (transpiration) and vevapwet (interception) !! Preliminaries DO ji = 1, kjpindex !! \latexonly !! \input{enerbilevapveg4.tex} !! \endlatexonly speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji))) !! \latexonly !! \input{enerbilevapveg5.tex} !! \endlatexonly xx(ji) = dt_sechiba * (un-vbeta1(ji)) * (un-vbeta5(ji)) * (qsol_sat_new(ji)-qair(ji)) * rau(ji) * speed * q_cdrag(ji) ENDDO DO jv=1,nvm DO ji=1,kjpindex !! 2.1 Calculate interception loss !! \latexonly !! \input{enerbilevapveg6.tex} !! \endlatexonly vevapwet(ji,jv) = xx(ji) * vbeta2(ji,jv) ! !! 2.2 Calculate transpiration !! \latexonly !! \input{enerbilevapveg7.tex} !! \endlatexonly transpir (ji,jv) = xx(ji) * vbeta3(ji,jv) transpot(ji,jv) = xx(ji) * vbeta3pot(ji,jv) END DO END DO IF (printlev>=3) WRITE (numout,*) ' enerbil_evapveg done ' END SUBROUTINE enerbil_evapveg END MODULE enerbil